Main

Enzymes rely on intricate intramolecular interactions between amino acids to organize the active site and achieve efficient catalysis. Rewiring these interactions through evolution often leads to unexpected and non-additive effects on protein fitness in a phenomenon known as epistasis1,2,3,4. Positive epistasis, wherein the effect of mutations is more beneficial when combined than in isolation, is highly desirable in protein engineering and often drives evolutionary trajectories. In contrast, negative epistasis is caused by mutations that act antagonistically, which adversely affects protein fitness and can be detrimental to protein engineering campaigns. Consequently, epistasis often dictates enzyme evolution by providing access to, or restricting, mutational paths1,2,3,4,5,6,7,8. Such non-additive interactions can be introduced by various factors, including changes in structural interactions and protein conformational dynamics. For example, epistasis can arise by rewiring the interactions of active site residues, thereby establishing new interactions with the substrate. Also, mutations that change conformational dynamics can improve enzymatic activity in a highly synergistic fashion9,10,11,12,13,14,15, for instance, by epistatically altering the dynamics of solvent-exposed loops to aid substrate entry or enhance active-site organization16,17,18.

While the link between structural changes and epistasis has been intensively studied, the mechanistic relationship between epistatic mutations and the overall catalytic cycle is often overlooked. Concerning epistasis, mechanistic studies may be particularly important because evolution typically enhances the slowest steps in the catalytic cycle19,20,21,22, which might change the rate-limiting step resulting in non-additive effects. Here we hypothesize that studying how mutations impact each step in the reaction could reveal novel mechanisms of epistasis, which would improve the predictability of evolution and provide a better understanding of the overall permissiveness of adaptive landscapes.

The β-lactamase OXA-48 is an excellent model system for studying epistasis in the context of antimicrobial resistance development23,24,25. OXA-48 uses a catalytic serine (S70) to cleave β-lactams in a three-step mechanism comprising enzyme–substrate complex formation (ES), formation of an acyl–enzyme intermediate (EI), and hydrolytic product release (E + P) (ref. 26). OXA-48 confers resistance to many carbapenem and penicillin β-lactams, but only slowly hydrolyses oxyimino-cephalosporins such as ceftazidime (CAZ)27. Low catalytic activity for CAZ hydrolysis has been attributed to the substrate’s bulkiness, potentially requiring sampling of alternative loop conformations to promote cephalosporin binding and hydrolysis16,28. While we have recently demonstrated that single mutations in OXA-48 can result in low-level CAZ resistance by increasing the flexibility of active site loops23,24,25, the mutational effects on the overall catalytic cycle and their potential for epistasis remain elusive.

In this Article, we used directed evolution to study the mechanistic drivers of epistasis during the adaptation of an antimicrobial resistance gene. After subjecting the β-lactamase OXA-48 to iterative rounds of mutagenesis and selection, we constructed an adaptive fitness landscape of the introduced mutations that identifies positive epistasis as a key driver for enzyme evolution. By combining biochemical, structural and computational methods, we reveal how evolution alters the protein’s conformational dynamics to accelerate substrate binding and hydrolysis, thereby changing the rate-limiting step and introducing epistasis. This detailed understanding of mutational effects on the whole catalytic cycle is crucial for predicting epistatic interactions1,2,3,4, which is relevant to the evolution of enzymatic activity and the design and engineering of novel enzymes29.

Results

Evolution of OXA-48 is driven by positive epistasis

To investigate how epistasis drives the evolution of OXA-48, we performed five cycles of directed evolution starting from the wild-type OXA-48 (wtOXA-48), using error-prone polymerase chain reaction (PCR) mutagenesis followed by selection on increasing CAZ concentrations (Fig. 1a). Variants arising along the evolutionary trajectory were characterized by their half-maximal inhibitory concentrations (IC50) from antibiotic dose–response growth curves. Five mutations accumulated along the evolutionary trajectory, F72L → S212A → T213A → A33V → K51E, resulting in Q5 conferring 43-fold increased CAZ resistance in Escherichia coli (Fig. 1b and Supplementary Table 1). Although minimum inhibitory concentrations (MICs) are typically used in clinical contexts, we evaluate functional activity on the basis of IC50 values due to their higher resolution (twofold resolution; see Supplementary Table 1 for MIC values)23. The mutations acquired during evolution cluster around the active site or structural elements known to alter substrate specificity, such as the nucleophilic S70 or the Ω (D143–I164) and β5–β6 (T213–K218) loops (Fig. 1c)27. Interestingly, F72L and T213A have been described in environmental and clinical isolates30,31, suggesting that laboratory evolution reflects the natural adaptation of pathogens to some degree. K51E increased resistance development by only 1.1-fold compared with Q4. Thus, we focused analysis on the molecular origins of epistasis in Q4 (A33V/F72L/S212A/T213A), which conferred a 40-fold higher CAZ resistance over wtOXA-48 in E. coli (Supplementary Table 1).

Fig. 1: Positive epistasis drives the evolution of OXA-48.
figure 1

a, During directed evolution of OXA-48, selection for resistance against the oxyimino-cephalosporin CAZ was performed at increasing CAZ concentrations from 0.5 to 14 µM. b, CAZ resistance conferred by OXA-48 was improved 43-fold over five rounds of evolution (see Supplementary Table 1 for all IC50 values). c, Mutations acquired during evolution, such as F72L, S212 and T213A, cluster around the active site serine (S70) and the Ω and β5–β6 loops that affect substrate specificity (purple). d, The adaptive landscape of the mutations found during evolution shows high epistasis. Each node represents a unique variant indicated by single-letter amino acid codes. Values within each node reflect the CAZ IC50 fold change relative to wtOXA-48. Purple arrows indicate the trajectory followed during evolution (see Supplementary Table 2 for all IC50 values). e, Comparison of the effects of single mutations (grey, F72L; dark grey, S212A; black, A33V; no expected effect for T213A) on the IC50 fold changes along the evolutionary trajectory reveals a high degree of epistasis (purple). f, Comparison of the fold-change improvements relative to the previous variants reveals diminishing returns (purple area) in CAZ resistance.

Source data

Full size image

To study the interplay between mutations acquired along the evolutionary trajectory, we constructed an adaptive landscape of all 16 mutational combinations based on the four mutations in Q4 and determined their IC50 values when produced in E. coli (Fig. 1d and Supplementary Table 2). Overall, positive epistasis greatly shaped the evolution of Q4 and resulted in a resistance increase of 40-fold in contrast to the 3.4-fold predicted increase for strictly additive gains from single-point mutations (Fig. 1e and Supplementary Fig. 1a). To quantify the apparent positive epistasis, we analysed the contribution of each mutation to the IC50 fold change across every possible genetic background (Supplementary Fig. 1b). Strikingly, epistasis is primarily driven by interactions with F72L, the first mutation acquired during evolution and the only single-point mutation that significantly (twofold) increased resistance in the wild-type background (analysis of variance (ANOVA), degrees of freedom (d.f.) 4, P < 0.001; Supplementary Table 3). For example, the combination of F72L with either S212A or T213A (referred to as alanine mutations hereafter) confers 8.2-fold and 11.7-fold higher resistance than expected (Supplementary Fig. 1c). We note that changes in thermostability (Tm) did not drive epistatic adaptation (Supplementary Fig. 1d and Supplementary Table 4). Except for F72L, which reduced the Tm by 6 °C, the selected mutations barely affected the thermostability. Thus, the observed epistasis does not stem from changes in protein stability, but from specific intramolecular interactions that affect the enzyme’s function. Interestingly, we observed diminishing returns epistasis in resistance development for higher-order variants along the evolutionary trajectory (Fig. 1f). As described above, F72L increased the IC50 by twofold, followed by a fivefold gain in IC50 through S212A. All later mutations, however, provided lower improvements than expected. For instance, addition of T213A to F72L/S212A resulting in Q3 conferred only a 3-fold increase in resistance instead of the 6.3-fold effect of T213A in F72L. Given that the alanine mutations are located in neighbouring positions, they probably have redundant effects on the CAZ IC50. Thus, their effects do not combine additively and lead to lower-than-expected improvements.

Evolution selects for an epistatic burst phase

To understand the molecular origins of epistasis, we studied the reaction kinetics of the OXA-48 variants (Fig. 2a). We monitored the conversion of CAZ by wtOXA-48, F72L and Q4 using a stopped flow and discovered that all enzymes catalysed the conversion of CAZ with an initial activity burst (Fig. 2b and Supplementary Fig. 2). At 400 µM CAZ, wtOXA-48 possessed a 1.5-fold higher rate in the burst phase than in its subsequent steady-state phase. Interestingly, the burst phase became more pronounced during evolution. F72L and Q4 show 4.6-fold and 48-fold higher burst-phase rates than their respective steady states. Evolution thus selectively improved the burst-phase over the steady-state activity.

Fig. 2: Kinetic changes drive the evolution of OXA-48.
figure 2

Kinetics are shown for wtOXA-48 (red) and the corresponding variants F72L (purple) and Q4 (blue, A33V/F72L/S212A/T213A). a, Enzymatic hydrolysis of β-lactams proceeds via an enzyme–substrate complex (ES), formation of an acyl–enzyme intermediate (EI), and hydrolytic deacylation (E + P). b, Evolution amplified the super-stoichiometric burst behaviour of OXA-48 (CAZ concentration: 50–400 µM, light to dark colours). c, In vitro burst-phase activities correlate well with the in vivo IC50 fold changes. The line represents the Pearson correlation, and the error bands display the 95% confidence interval. d, Michaelis–Menten kinetics of the burst phase determined at 4 °C. e, Substrate binding was measured by W-fluorescence and was substantially accelerated during evolution (CAZ concentration: 100–1,200 μM, light to dark colours). f, Comparison of k1 and kcat/KM between wtOXA-48 (red), F72L (purple) and Q4 (blue) reveals that binding is no longer rate-limiting in the burst phase of Q4 (determined at 25 °C; point above the diagonal line at k1 = kcat/KM indicates that binding is not rate-limiting).

Source data

Full size image

To obtain kinetic parameters for the burst phase for a range of variants along the evolutionary trajectory, we assayed substrate conversion with a microtitre plate reader at 4 °C (Table 1 and Supplementary Fig. 3). As expected, F72L improved the kcat/KM of the burst phase by 8-fold compared with wtOXA-48, while S212A and T213A provided only marginal improvements (~1.5-fold). Similar to the IC50 effects, F72L displayed strong pairwise epistasis with either S212A or T213A at the kinetic level and improved kcat/KM by 70-fold and 60-fold compared with wtOXA-48, respectively. Positive epistasis shaped the evolution of kcat/KM, culminating in a 470-fold and 800-fold improvement in Q3 and Q4, respectively. The burst-phase kcat/KM values excellently correlate with the in vivo IC50 values (R2 = 0.97; Pearson correlation, d.f. 7, P < 0.001; Fig. 2c), Notably, that correlation is much stronger than the corresponding effect observed for the steady-state kcat/KM values (R2 = 0.85; Pearson correlation, d.f. 7, P < 0.001; Table 1 and Supplementary Figs. 3–5). Fast degradation of CAZ is crucial to bacterial fitness, although in vivo resistance could also be affected by other factors, such as in vivo stability, protein production and translocation into the periplasm32,33. The excellent correlation of the kcat/KM and IC50 values in the burst phase, however, suggests that the resistance increases are dominated by enhanced burst-phase activity.

Table 1 IC50 and burst-phase catalytic parameters of OXA-48 variantsa,b
Full size table

The kinetics of the burst phase unveiled fundamentally different catalytic effects between F72L and either the S212A or T213A variant (Fig. 2d and Supplementary Fig. 3). Michaelis–Menten kinetics of wtOXA-48 did not show saturation due to its high KM (≫300 µM). Introduction of F72L led to saturation kinetics and substantially decreased the KM to 165 ± 15 µM, suggesting an improvement in substrate binding. In contrast, the Michaelis–Menten plots remained linear upon introducing S212A or T213A and only showed marginal increases in kcat/KM (1.6-fold). Interestingly, the double mutants F72L/S212A and F72L/T213A maintained the relatively low KM value observed for F72L while also increasing kcat by five- to sevenfold (Table 1 and Supplementary Fig. 3). Thus, epistasis between F72L and S212A or T213A probably originates from an interplay between improved binding and catalysis.

Burst phases have been previously observed in β-lactamases, where fast formation of an acyl–enzyme intermediate during the first turnover is frequently followed by rate-limiting deacylation34,35,36. Surprisingly, the burst phases observed here have amplitudes much larger than a single turnover. For example, Q4 displayed a burst-phase amplitude corresponding to 3.8 turnovers at 400 µM CAZ (Fig. 2b). Such super-stoichiometric bursts cannot be explained by a simple two-step reaction mechanism comprising a fast followed by a slow step34,35,36. Instead, super-stoichiometric bursts are probably caused by the inactivation of the enzyme during substrate conversion over several catalytic cycles34,35,36. We note that Q4 did not show any evidence for a classical burst stemming from fast acylation and slow deacylation at 1,200 µM CAZ, even at the millisecond timescale. While the exact mechanism underlying the biphasic behaviour is unclear, substrate-induced inactivation becomes faster during evolution, indicating that the mutations acquired during evolution increase conformational flexibility (Supplementary Fig. 2). The inactivation rate also increases with substrate concentrations, which argues for an induced-fit mechanism in which the presence of the substrate or enzyme-bound intermediate triggers inactivation37. Notably, substrate-induced inactivation appeared fully reversible, as indicated by activity assays after incubating Q4 with CAZ (Supplementary Fig. 6a). In addition, inactivation is not caused by a change in the oligomeric state, as shown by both dynamic light scattering and size exclusion chromatography (Supplementary Fig. 6b,c). Thus, turnover apparently triggers a reversible conformational change to a less active state, resulting in the observed burst.

The dependence of the burst phase on the CAZ concentration provides important insights into its role in evolutionary adaptation. With decreasing CAZ concentrations, the turnover rates of the burst and steady-state phases become similar, and the burst phase amplitude decreases. The burst thus became less pronounced at substrate concentrations down to 50 µM at 25 °C (Fig. 2b). This is probably because, at lower concentrations, substrate-induced inactivation is rarer, and any enzyme that is deactivated has more time to recover to the active state before encountering another substrate molecule. Notably, selection was performed under even lower substrate concentrations (≤14 µM) compared with the in vitro kinetic analysis. Thus, the enzymes probably remained predominantly in the burst state under the selection conditions. Taken together, while evolution boosted resistance at the evolutionarily relevant CAZ concentrations, a catalytic bottleneck emerged at high and physiologically irrelevant concentrations that limits activity after the initial burst.

Faster binding and reaction drive positive epistasis

Based on our Michaelis–Menten kinetics, we hypothesized that epistasis in OXA-48 resulted from an interplay between binding and the chemical reaction. F72L apparently unlocked the accessibility of the fitness landscape by substantially accelerating substrate binding, which allowed S212A and T213A to take effect and further boost catalysis (Fig. 2). We dissected this relationship by assaying CAZ binding (k1 and k−1; Fig. 2a) and solvent isotope effects in wtOXA-48, F72L and Q4. By focusing on wtOXA-48 and F72L, we aimed to reveal the mechanistic role of F72L and its ability to recruit subsequent mutations. Comparison with the evolved variant Q4 then allowed us to shed light on the combinatorial effect of F72L with the other mutations.

CAZ binding was assayed by monitoring changes in protein tryptophan fluorescence. In agreement with our burst phase results, CAZ binding progressed with biphasic kinetics (Fig. 2e), where the fast phase probably reflects binding to the burst phase ensemble. The slow phase indicates the presence of a second state that is probably related to the steady-state ensemble. Since the slow phase occurs on a similar timescale as enzyme deactivation and the chemical reaction (Fig. 2b,e), it probably reflects a combination of these processes and binding to the less active state. Given that the burst phase is most probably the physiologically relevant phase, we decided to focus our analysis on the fast-binding phase (Table 2). For wtOXA-48, CAZ binding (k1) is 1.5-fold slower than the burst-phase kcat/KM, which suggests that binding—and not bond-breaking or product release—limits catalytic efficiency in the wild type. Notably, k1 constantly increases during evolution (3.2-fold in F72L and 6,000-fold in Q4). In Q4, k1 is tenfold faster than kcat/KM, indicating that evolution shifted the catalytic bottleneck from substrate binding to either the chemical steps or product release (Fig. 2f). Counterintuitively, the KD determined from k−1/k1 did not follow the trends observed for KM and was 2.6-fold higher for F72L than for wtOXA-48. This observation serves as a reminder that KM is an imperfect approximation for substrate affinity and supports our hypothesis that epistasis in OXA-48 has kinetic and not thermodynamic origins. Epistasis is most probably driven by a change in the rate-limiting step resulting from faster substrate binding (k1) and not by improved substrate affinity (KD).

Table 2 Burst-phase binding and catalysis kineticsa
Full size table

To further test our hypothesis that evolution shifted the rate-limiting step, we determined solvent isotope effects by assaying product formation at 400 μM CAZ in 80% D2O (Table 2). Solvent isotope effects >1 indicate rate-limiting deprotonation during the hydrolysis of the acyl–enzyme complex (k3). The burst-phase isotope effects for wtOXA-48 and Q4 are close to 1, whereas a small isotope effect was observed for F72L (1.4). These isotope effects suggest that hydrolysis of the acyl–enzyme affects the overall rate in F72L, while deacylation is not rate-limiting in wtOXA-48 and Q4. In contrast to the burst-phase effect, F72L and Q4 had pronounced isotope effects of 2.8 and 2.5 in their steady state. Since our data on the evolution of super-stoichiometric burst phases suggest that a conformational change triggers inactivation, the observed increase in isotope effects from the burst phase to the steady state probably indicates that hydrolysis becomes rate-limiting in the less active state. Alternatively, the conformational equilibrium itself could have an isotope effect causing the differences in activity.

Our combined data on binding and activity demonstrate how evolution successively optimized the catalytic cycle and gradually changed the rate-limiting step. In the burst phase, which is probably the physiologically relevant phase, the slowest step appears to be substrate binding in wtOXA-48, binding or deacylation in F72L, and acylation in Q4. Our analysis indicates that the apparent change in rate-limiting step causes the observed epistatic effects. In other words, because substrate binding is slow in wtOXA-48, it is likely that neither S212A nor T213A substantially affect activity in the wild-type background and require F72L to accelerate substrate binding for their improvements to take effect fully.

F72L and alanine mutations orthogonally tune dynamics

To understand how evolution modulated the conformational dynamics of OXA-48, we determined the crystal structures of F72L, Q5, and Q5 covalently bound to CAZ (PDB IDs: 8PEA, 8PEB and 8PEC; Supplementary Table 6 and Supplementary Fig. 7). Comparison of the apo structures of F72L and Q5 with wtOXA-48 (PDB ID: 4S2P (ref. 38)) provided insights into the structural role of F72L during evolution (Supplementary Fig. 8). In wtOXA-48, F72 is embedded within an aromatic pocket formed by Y144, F156, and W157, in close proximity to the active site S70. The F72L mutation disturbs this aromatic network, which probably reduces the thermostability (Supplementary Fig. 1d) and allows the Ω loop harbouring F156 and W157 to adapt a different conformation. Despite the introduction of S212A and T213A in the β5–β6 loop, the overall shape of this loop is preserved in Q5. Similarly, the overall active site architecture (for example, S70 and K73) was maintained in all variants (Supplementary Fig. 7). We note that upon CAZ binding in Q5 (Q5-CAZ), the Ω loop becomes highly disordered to the point that it could not be refined anymore (Supplementary Fig. 7d), which suggests that incubation with CAZ leads to population of a state that could potentially be less active.

To estimate how the conformational dynamics changed during evolution, we performed ensemble refinements based on the crystal structures of wtOXA-48, F72L and Q5 (Fig. 3a). Despite their differences in resolution (Supplementary Table 6), both F72L and Q5 revealed an increased flexibility primarily resulting in the Ω loop adopting various alternate conformations. In addition, we performed atomistic molecular dynamics (MD) simulations of apo wtOXA-48 and Q4 that support this change in dynamics and reproduce the enhanced flexibility, particularly of the Ω loop (Supplementary Fig. 9). This increase in Ω-loop dynamics probably accelerates substrate binding and caused the observed change in rate-limiting step23,39.

Fig. 3: Evolution of a catalytically superior ensemble.
figure 3

a, Ensemble refinement of wtOXA-48, F72L and Q5 reveals increased mobility of the Ω loop. b, ΔRMSF values relative to wtOXA-48, F72L and Q4 from MD simulations reproduce the increased flexibility of the Ω-loop region (see Supplementary Fig. 11 for other variants). Error bands indicate the standard error of the mean. c, PC and cluster analysis show that F72L modulates the conformational landscape in ways likely to accelerate binding displayed as the population shaded in blue (arrow added to highlight change in populations; see Supplementary Figs. 12 and 13 for other variants). d, Cluster representatives indicate that evolution (Q4 variant in blue versus wtOXA-48 in grey) displaced the Ω loop and adjacent α-helix as indicated by the arrow. e, Dynamical correlation analysis reveals that the movement of the acylated S70 becomes tightly coupled with the protein scaffold, particularly the oxyanion hole, by means of the alanine mutations. In contrast, F72L predominantly decreases the interaction of S70 with the Ω loop (increased and decreased correlations relative to wtOXA-48 are shown in blue and red lines, respectively). Line width corresponds to the strength of the correlation. Only statistically significant changes compared with wtOXA-48 are shown (t-test, α = 0.05; see Supplementary Fig. 15 for other variants).

Source data

Full size image

To understand how the conformational dynamics affected the substrate-bound state along the evolutionary trajectory, we performed MD simulations of various OXA-48 variants in the acyl–enzyme complex. These models were constructed on the basis of the Q5-CAZ crystal structure with the unresolved Ω loop taken from the apo wtOXA-48 structure (Supplementary Fig. 7). To avoid biases from the loop conformation, we also analysed Q4 using the Ω loop from the apo Q5 structure, resulting in virtually identical results (Supplementary Figs. 9–15). To explore how evolution affected the enzyme flexibility, we determined per-residue root-mean-square fluctuations (RMSFs) for all variants (Fig. 3b and Supplementary Fig. 11). As expected from the MD simulations of the apo state, introduction of F72L predominantly increased the RMSF values of the Ω-loop region. Interestingly, the α-helix preceding the Ω loop increased in rigidity, indicating that the added space provided by F72L allowed that helix to pack tighter against the protein. Furthermore, the solvent-exposed loops covering the active site became more rigid upon introduction of the alanine mutations, which signals higher organization that probably aids transition-state stabilization.

To further dissect how the mutations affect the protein dynamics, we analysed the conformational landscape of the enzyme by principal component (PC) and cluster analysis (Fig. 3c and Supplementary Figs. 12 and 13). As expected, F72L significantly modulated the conformational landscape, while introducing the S212A or T213A mutations barely changed the sampled conformational space (Fig. 3c and Supplementary Fig. 13). Notably, F72L allows the protein to populate a different conformational ensemble, in which the space created by the loss of the phenylalanine sidechain allows the neighbouring α-helix to slide toward the centre of the protein while increasing the conformational freedom of the Ω loop (Fig. 3d and Supplementary Fig. 14). This movement is accompanied by a tightening of an H-bonding network below the Ω loop involving T71, Y144 and Q169 (Supplementary Fig. 14a).

While our MD analysis supports that F72L accelerates substrate binding and thereby shifts the rate-limiting step, the structural role of S212A and T213A is more elusive. Since these mutations did not affect the sampled conformational space, we hypothesized that dynamical cross-correlation analysis might allow us to identify regions that show correlated movements with the substrate covalently bound to S70 (Fig. 3e and Supplementary Fig. 15). Dynamical correlations signal communication between residues, with an increase in correlation strength suggesting a stronger interaction network. When analysing the changes in the correlation of the acylated S70 with the rest of the protein, we observed distinct dynamical changes wrought by either F72L or the alanine mutations on the protein scaffold. F72L primarily decreased the dynamical correlation between S70 and the Ω loop, which agrees with the increased flexibility and improved substrate binding conferred by F72L. In contrast to F72L, S212A and T213 had an entirely different effect on the protein dynamics. Notably, S212A and T213A enhanced the correlation of S70 with the neighbouring β-strands that harbour the oxyanion hole (Fig. 3e and Supplementary Figs. 14b and 15)40. This increased correlation could be rooted in the increased rigidification of the β5–β6 loop, harbouring the S212A and T213A mutations, which in turn allows for a tighter interaction of the loop with the substrate (Supplementary Figs. 11 and 15). Dynamical correlations between the nucleophile and the oxyanion hole may explain increases in kcat, because they signify better preorganization and the ability to stabilize the transition state19,41,42. Overall, the effects of F72L and the alanine mutations on the correlations were largely orthogonal, like their effects on conformational sampling. Such orthogonal dynamical relationships are probably central to the epistatic effects not only in OXA-48, but also for the evolution of other natural and designer enzymes29,43.

Discussion

Unraveling the molecular mechanisms underlying epistasis is vital for understanding enzyme evolution1,2,3,4. Many previous studies on epistatic enzyme evolution focused on exploring epistatic mechanisms from a structural perspective, for instance, by demonstrating direct interactions between mutations or synergistic effects on conformational dynamics9,13,44. In contrast, our observations highlight that the molecular basis for intramolecular epistasis can also be rooted in changes in the catalytic cycle. We describe how distinct mutations shifted the rate-limiting step from substrate binding to the chemical reaction in the β-lactamase OXA-48, thereby causing strong phenotypic epistasis (Figs. 1 and 2). Intriguingly, the adaptive mutations introduced into OXA-48 are structurally orthogonal, but mechanistically epistatic: F72L induces dynamical and structural perturbations that are largely independent of the S212A and T213A, and vice versa (Fig. 3). Since CAZ binding is rate-limiting in wtOXA-48, F72L must be incorporated first to accelerate binding and unlock the effect of S212A and T213A on the chemical step. Here, evolution shifted the slowest step from binding to chemistry. We note that shifting the rate-limiting step between other stages of the catalytic cycle could likewise result in epistasis. While genetic context-dependent effects on kcat and KM have been previously reported45, deciphering their underlying epistatic relationship has remained challenging. Here, gaining detailed insights into the origins of epistasis was only possible by in-depth characterizations of the burst-phase and steady-state kinetics, isotope effects and dynamical analysis. To the best of our knowledge, this is the first report on how changing the reaction bottleneck from binding to the chemical step leads to synergy in evolution. We hypothesize that similar effects are likely to be of vast importance in the evolution of other biocatalysts.

Catalytically perfect enzymes, in which the reaction rate is limited only by substrate diffusion, display catalytic efficiencies (kcat/KM) of >108 M−1 s−1. Although the kcat/KM of wtOXA-48 is orders of magnitude below the diffusion limit, our analysis unexpectedly revealed that CAZ binding is rate-limiting in this variant (Fig. 2). This stands in contrast to many in other serine β-lactamases where acylation and deacylation are often catalytic bottlenecks25,46. The extended size of CAZ has led to the hypothesis that its binding is more challenging than that of other β-lactams16,46. OXA-48 furthermore differs from other serine β-lactamases in that it accommodates the carboxylate group adjacent to the oxyimino-moiety of CAZ within its active site (Supplementary Fig. 16)24,28, which probably additionally slows down CAZ binding resulting in rate-limiting complex formation. Understanding how individual steps affect overall catalytic activity is crucial for understanding and predicting enzyme evolution and designing more efficient drugs that exploit mechanistic bottlenecks.

Our findings show the selection of a super-stoichiometric burst phase that correlates with increasing levels of antibiotic resistance (Fig. 2). The origin of such a burst phase cannot be described by a simple two-step reaction mechanism and is probably driven by the selection of pre-existing conformational sub-states and changes in conformational dynamics (Fig. 3)13,18,47,48. The intrinsic limitations of the OXA-48 system, comprising slow substrate binding and various possible substrate-bound states that can undergo inactivation, prevented the exact determination of the catalytic mechanism. Nonetheless, our kinetic and structural analysis revealed that increased loop flexibility and tighter substrate interactions drove the evolutionary improvements. While fully understanding these structural effects probably requires advanced NMR or mass spectroscopy, we conclusively show that the evolution of burst-phase kinetics introduces a kinetic bottleneck that restricts efficiency at elevated CAZ concentrations (Fig. 2b). Importantly, this bottleneck is irrelevant at the comparably low substrate concentrations used during selection, where the burst behavior vanishes. We hypothesize that offsetting activity enhancements under physiologically irrelevant and high substrate concentrations probably relieved selection pressure and facilitated improvements at the CAZ concentrations during selection.

In conclusion, our study demonstrates how mutations that orthogonally affect different steps along the catalytic cycle can cause positive epistasis by shifting the catalytic bottleneck. Traditionally, epistasis analysis has focused on structural and dynamical interactions. Our research underscores the importance of considering mechanistic effects throughout the catalytic cycle to grasp the origins of epistasis fully. Especially with regard to β-lactamases and other resistance-mediating enzymes, a detailed comprehension of epistasis is essential to anticipate the emergence of novel pathogens and predict resistance phenotypes. Such a detailed understanding will not only help fight antimicrobial resistance, but also aid the design of drugs exploiting catalytic bottlenecks and advance enzyme engineering.

Methods

General material

Luria-Bertani (LB) agar, LB broth, chloramphenicol, ampicillin and CAZ were purchased from Sigma-Aldrich. Primers (P) used for this study are presented in Supplementary Table 7. All cloning enzymes were purchased from Thermo Fisher Scientific, if not stated otherwise. The E. coli E. cloni 10G (MP21-5) was obtained from Lucigen. All strains used and constructed in this study are presented in Supplementary Table 8. Kinetic data were fitted using Prism v. 9.0 (GraphPad Software).

Directed evolution and cloning

The construction of the low copy number vector pUN-blaOXA-48 (pA15 origin; 10–20 copies per cell) was previously published49. Error-prone PCR was performed using 10 ng pUNE-4-blaOXA-48, GoTag (Promega), 25 mM MgCl2 (Promega), 10 µM P7/P8 and either 50 µM oxo-dGTP or 1 µM dPTP (Jena Bioscience). PCR products were DpnI digested for 1 h and 37 °C, and 5 ng of each product was used for a second PCR, which was performed as described above, but without mutagenic nucleotides. PCR amplicons were digested using NcoI, XhoI and DpnI for 1 h at 37 °C, purified for ligation with the vector backbone and transformed into MP21-5. To insure a sufficient mutational depth, we aimed for library sizes of at least 5,000 colonies, which was determined by plating on agar plates supplemented with 25 mg l−1 chloramphenicol LB and one to two amino acid changes per round of evolution, which was determined by Sanger sequencing (Azenta).

The fitness landscape was constructed in the pUN vector background using Goldengate cloning and the corresponding primers in Supplementary Table 7. In short, we performed whole vector amplification followed by digestions with LguI and DpnI for 1 h 37 °C. Ligations were performed using 10–20 ng of DNA for 1 h at room temperature using T4 ligase and transformed into MP21-5, and clones were grown on 25 mg l−1 chloramphenicol LB agar plates and verified using Sanger sequencing (Azenta).

For protein expression, OXA-48 variants were subcloned into a pDEST-17 (pURR) expression vector without the leader sequence and with a 6-His-tag using P2/P37 and P35/P36. Amplicons were digested using NotI and XhoI, ligated as described above and transformed into MP21-5. Vectors were isolated using the plasmid miniprep kit (Qiagen). pURR expression vectors were transformed into E. coli BL21 AI. Clones were selected on agar containing ampicillin 100 mg l−1 and verified using Sanger sequencing (Azenta).

Selective plating

MP21-5 cultures harbouring either pUNE-4-blaOXA-48 or a library of OXA-48 were plated on LB agar plates containing increasing concentrations of CAZ and grown overnight at 37 °C (Supplementary Table 1). Up to eight colonies grown on the highest concentrations were recovered and their genotype characterized by Sanger sequencing (Azenta). Before determining their IC50 values, the corresponding mutant alleles were subcloned into an isogenic pUN vector backbone and transformed into MP21-5, to exclude mutational effects outside of the target gene.

IC50 determination

IC50 values were determined as described previously23,49. In brief, cultures were grown to full density under 700 rpm shaking overnight at 37 °C. Overnight cultures were diluted in phosphate-buffered saline to a density of 106 cells ml−1 and used to inoculate a 384-well plate (Thermo Fisher Scientific) with a CAZ gradient (0 to 32 mg l−1) at a final cell density of 105 cells ml−1. Plates were incubated statically at 37 °C for 20 h. The absorbance was determined as OD600 using an Epoch spectrophotometer (Biotek). Dose–response curves and their IC50 value were determined on the basis of a non-linear fit.

Protein expression and purification

Cultures of E. coli BL21AI harbouring modified pURR expression vector with blaOXA-48 or mutant alleles were grown in terrific broth supplemented with 100 mg l−1 ampicillin at 30 °C and 220 rpm. Protein expression was induced by adding l-arabinose (Sigma-Aldrich) to a final concentration of 0.2% when the cultures reached an OD600 of 0.4. Cultures were expressed for 16 h at 15 °C and centrifuged at 4 °C for 30 min, and the cell pellets were stored at −20 °C for purification. Protein purification was performed using HisPur Ni-NTA spin columns (Thermo Fisher Scientific) as published previously23,50. For crystallization, the His-tag was cleaved overnight at 4 °C using in-house produced tobacco etch virus (TEV) protease and the TEV-cleaved product was purified through an additional round of HisPur Ni-NTA columns (Thermo Fisher Scientific).

Fluorescence-based thermostability

Thermostability was performed as previously published using purified OXA-48 enzymes, containing 6-His-tag and TEV cleaving site23. In 50 mM HEPES (VWR), pH 7.5 including 50 mM potassium sulfate (Honeywell), enzymes were diluted to 0.2 mg ml−1 and mixed with 5xSYPRO orange (Sigma-Aldrich). Using an MJ minicycler (Bio-Rad), a temperature gradient (25 °C to 70 °C) was performed with a heating rate of 1 °C min−1. All experiments were performed in triplicates, and the melting temperatures (TM) were determined as the inflection point of the melting transition found from the first derivative.

Burst-phase and steady-state enzyme kinetics

Room-temperature burst kinetics were obtained under burst-phase conditions using an SX20 stopped flow (Applied Photophysics) by monitoring substrate depletion by absorbance at 260 nm. Enzyme and substrate were mixed 1:1 at 25 °C and in 0.1 M phosphate buffer (Sigma-Aldrich, pH 7.2) supplemented with 50 mM NaHCO3 (Sigma-Aldrich). Burst kinetics were assayed at final enzyme concentrations of 10 μM and final substrate concentrations varying between 50 and 400 μM (equation (1)).

$$\frac{P}{{E}_{0}}={v}_{{{\mathrm{steady}}}}\times t-\left({v}_{{{\mathrm{steady}}}}-{v}_{{{\mathrm{burst}}}}\right)\times \left(1-{e}^{-k\times t}\right)/k$$
(1)

Catalytic parameters (kcat, KM and kcat/KM) were determined under burst-phase and steady-state conditions using CAZ (Δξ = −9,000 M−1 cm−1) at 260 nm by measuring the initial enzymatic reaction rate in an Epoch plate-reader (Biotek). Burst phase rates were determined at 4 °C, and steady-state parameters were determined at 25 °C. Reactions rates were obtained in at least duplicates at a final enzyme concentration of 1 μM (final assay volume of 100 μl). Ultraviolet-transparent 96-well plates (Corning) were used. Assays were performed in 0.1 M phosphate buffer (Sigma-Aldrich, pH 7.2), supplemented with 50 mM NaHCO3 (Sigma-Aldrich).

Kinetic isotope effects

Kinetic isotope effects (KIE) were determined at 25 °C using an SX20 stopped flow (Applied Photophysics) from burst kinetics obtained at 400 μM CAZ. KIEs were calculated from the ratio of the rate in 80% D2O (Sigma-Aldrich; kD) and water (kH; equation (2)).

$${{\mathrm{KIE}}}=\frac{{k}_{{\rm{H}}}}{{k}_{{\rm{D}}}}$$
(2)

Binding kinetics

CAZ binding kinetics were assayed using an SX20 stopped flow (Applied Photophysics) by tryptophane fluorescence, with an excitation wavelength of 280 nm, a 305 nm lower cutoff emission filter and a 0.2 mm excitation pathlength. Enzyme and substrate were mixed 1:1 at a final enzyme concentration of 1 μM, and the final substrate concentration was varied between 100 and 1,200 μM. Binding was assayed at 25 °C and in 0.1 M phosphate buffer (pH 7.2) supplemented with 50 mM NaHCO3 (Sigma Aldrich). Binding was recorded on a log timescale, and rates were obtained by global fitting of the observed biphasic curves to a double-exponential decay (equation (3)).

$${{\mathrm{Fluorescence}}}={A}_{1}\times {\exp }^{-{k}_{1}* t}+{A}_{2}\times {\exp }^{-{k}_{2}* t}+c$$
(3)

where k1 and k2 are the observed binding rates and A1 and A2 are the amplitudes of the two signals with an offset of c. To enable global fitting, k1 and k2 were fitted to a linear equation each (equations (4) and (5)), where k1,on, k2,on, k1,off and k2,off were shared between all datasets.

$${k}_{1}={k}_{1,{\rm{on}}}\times c\left({\rm{CAZ}}\right)+{k}_{1,{\rm{off}}}$$
(4)
$${k}_{2}={k}_{2,{\rm{on}}}\times c\left({\rm{CAZ}}\right)+{k}_{2,{\rm{off}}}$$
(5)

Size exclusion chromatography

Changes in the molecular size were studied using size exclusion chromatography: 50 nM, 10 μM wtOXA-48 and 10 μM of the mutant R189A/R206A mutants. R189A/R206A has been previously described to disrupt the dimer interface of OXA-48 and elutes later51. Separation was performed in 0.1 M phosphate buffer (Sigma-Aldrich, pH 7.0) using a Superdex 200 10/300 GL column with a flow rate of 0.5 ml min−1 at 4 °C. Elution was monitored by recording the absorbance at 280 nm.

Dynamic light scattering

The hydrodynamic radius of Q4 was determined using a zetasizer (Malvern Panalytical). To that end, 1 ml of 10 μM Q4 was assayed in either the presence or absence of 400 μM CAZ in 0.1 M phosphate buffer (Sigma-Aldrich, pH 7.0) supplemented with 50 mM NaHCO3 (Sigma-Aldrich). The experiment with CAZ was repeated after a 15 min incubation at room temperature.

Sequential mixing

To assess the reversibility of the burst phase, 20 μM Q4 was pre-incubated with 400 μM CAZ, and mixed with a second batch of 400 μM CAZ after a delay of 1,000, 2,000 and 3,000 s in an SX20 stopped flow (Applied Photophysics) by monitoring substrate depletion by absorbance at 260 nm.

Crystallography, structure determination and refinement

Crystallization was performed in a 1 μl hanging drop containing 10 mg ml−1 enzyme and mixed 1:1 with reservoir solution containing 0.1 M Tris, pH 9.0 (Sigma-Aldrich), and 28–30% polyethylene glycol mono ethylene ether 500 (Sigma-Aldrich) at 4 °C. Crystals were cryoprotected by using 15% ethylene glycol (Sigma-Aldrich) in addition to the reservoir solution, and subsequently frozen in liquid nitrogen. Diffraction data were collected on ID23-EH2 (F72L) and ID30B (Q5 and Q5-CAZ), European Synchrotron Radiation Facility, France, at 100 K, wavelengths 0.8731 Å (F72L) and 0.9763 Å (Q5 and Q5-CAZ), and the diffraction images were indexed and integrated using x-ray detector software52. For data scaling, AIMLESS was used53 and an overall high completeness and CC1/2 >0.5 and a mean intensity above 1.0 in the outer resolution shell was aimed for (Supplementary Table 6). Molecular replacement was performed using chain A of PDB ID 6Q5F (ref. 23) and the program PHENIX 1.12 (ref. 54). Parts of the model were manually rebuilt using Coot55. Average structure refinement and ensemble refinement was performed using PHENIX 1.12. PyMOL 1.8 was used for illustrations (Schrödinger).

MD simulations: system setup

MD simulations were set up analogously to our previous work23,28. Simulations of the apo enzymes were set up as described above for wtOXA-48 (PDB ID: 4S2P (ref. 38)) and Q4 based on the structure of Q5 (PDB ID: 8PEB). No additional restraints were employed for the apo simulations (see Supplementary Figs. 9 and 10 for all simulated structures). Acyl–enzyme structures of the OXA-48 variants with covalently bound CAZ were built on the basis of the structure of apo wtOXA-48 (PDB ID: 4S2P (ref. 38)), with CAZ added from the holo structure of Q5-CAZ (PDB ID: 8PEC; Supplementary Fig. 7). All OXA-48 variants were modelled on the basis of this structure using the mutagenesis tool in PyMOL, choosing the rotamer with the least steric clashes with surrounding atoms. For comparison, MD simulations for Q4 were also performed based on the holo structure of Q5-CAZ (PDB ID: 8PEC) and the Ω loop of apo wtOXA-48 (residues D143–I164). This variant is indicated as Q4 (template Q5) in Supplementary Figs. 9–13 and 15. The results from those simulations agreed qualitatively with those obtained from the Q4. The system was parametrized using tleap56, and enzymes were solvated in a 10.0 Å octahedral box of TIP3P water57,58 with net charge neutralized using the Amber uniform neutralizing plasma56. The ff14SB force field59 was used to describe the protein. Parameters for the carbamylated lysine (KCX) were previously obtained28 from restrained electrostatic potential fitting as implemented in the RED Server60. Parameters for the CAZ-acetylated serine were likewise obtained with the RED Server.

Several restraints were applied during the simulations to maintain a productive conformation. The restraints include a ≤4.0 Å distance restraint from the nucleophilic water to either the KCX base oxygen or CAZ carbonyl carbon using one-sided harmonic potentials. The distance of the KCX Nζ and catalytic Ser70 Oγ was likewise restrained to ≤4.0 Å. Lastly, flat-bottom potentials were applied to the CδCεNζCη (≤−130° and ≥−80°) and CγCδCεNζ (≤45° and ≥5 = 95°) dihedral angles. All restraint force constants were 10 kcal mol Å−2 during the equilibration MD and 100 kcal mol Å−2 during minimization.

MD simulations: simulations

All systems were minimized using 10,000 steps of steepest descent followed by 10,000 steps of conjugate gradient. During both minimization steps, the position of all protein atoms was restrained with a weight of 10 kcal mol Å−2. The minimization was subsequently repeated without positional restraints. Subsequently, the system was heated from 50 K to 300 K in 20 ps, and then simulated for 50 ns in the NPT ensemble saving a frame every 100 ps. Langevin dynamics were used with a collision frequency of 0.2 and a 2 fs time step. The Berendsen barostat was used with isotropic position scaling. All bonds involving hydrogens were constrained using the SHAKE algorithm. Twenty independent simulations were run per enzyme variant (for a total of 1.0 µs per variant). All calculations were performed with the Amber18 program package (sander.MPI for minimization and pmemd.cuda for MD simulations)56.

MD simulations: analysis

MD simulations were analysed using CPPTRAJ61. All analyses were based on Cα positions. The first 10 ns of each production MD run were excluded to allow sufficient time for system equilibration. Root mean square deviation values were calculated compared with the minimized starting structures. RMSFs were determined by first calculating an average structure for each replicate, aligning the trajectory against the rigid core of the average structure (residues 42–83, 112–131, 162–211 and 220–239), and then calculating the RMSF for each protein residue. Errors indicate the standard error of the 40 independent replicates.

For each variant, cluster and PC analysis were performed on their combined 40 replicates. To that end, a global average structure over all variants was first determined and all trajectories were aligned by Cα to the rigid core of the average structure (residues 42–83, 112–131, 162–211 and 220–239). Dynamical cross-correlation analysis was performed with CCPTRAJ by calculating the correlation matrix on the basis of the aligned structures. After the alignment, PC and cluster analyses were performed for all Cα atoms without re-aligning the structures. To perform the analyses in the same space, both the cluster and PC analysis were performed simultaneously for all variants. Clustering based on Cα root mean square deviation of the loops was performed with CPPTRAJ using the kmeans algorithm to split each trajectory into two clusters. PC analysis was performed using mdtraj62 and sklearn63.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.