MaterialsCommercial reagents are listed in Supplementary Data 2. Oligonucleotides (except for oligo pools) were ordered from Microsynth. Clonal genes and oligo pools were ordered from Twist Bioscience.Synthesis of 5-nitrobenzisoxazole (1)The synthesis of 5-nitrobenzisoxazole (1) was performed according to published protocols51. 1,2-Benzisoxazole (5 ml, 5.87 g) was added to concentrated H2SO4 (20 ml) at 0 °C until the solution turned yellow. A mixture of concentrated HNO3 (3.4 ml) and concentrated H2SO4 (1.3 ml) was added slowly at 0 °C, and the solution was stirred for 30 min. The reaction product was poured onto an ice/water mixture (1:1, 100 ml), and the crystals that formed were collected by filtration, washed with ice-cooled water and dried. The crude product was purified by normal-phase flash chromatography (RediSep column) with cyclohexane and ethyl acetate as mobile phases. The solvent was removed to yield 5-nitrobenzisoxazole (1) (3.05 g) as colorless needles (1H NMR (500 MHz, CDCl3) δ = 8.90 (d, 1 H), 8.73 (d, 1 H), 8.51 (dd, 1 H), and 7.76 (d, 1 H) and 13C NMR (500 MHz, CDCl3) δ = 164.4, 147.1, 144.8, 125.6, 121.9, 119.2, and 110.5).Plasmid constructsThe HG3 and HG3.17 constructs were ordered from Twist Bioscience as cloned genes in a pET28b vector between the NcoI/XhoI restriction sites containing a C-terminal His6-tag with a stop codon. The pelB-HG3 and pelB-HG3.17 variants used for screening were obtained by cloning the genes into the pET22b vector using NcoI/XhoI restriction enzymes, introducing an N-terminal pelB signal peptide and a C-terminal His6-tag. Sequences of all genes are reported in the Supplementary Information.Generation of complex single-site libraries using oligo poolsThe cloning procedure is illustrated in Extended Data Fig. 1.The oligo pool from Twist Bioscience was resuspended in MilliQ-water to yield a DNA concentration of 5 ng µl−1. The pool was amplified by PCR (25 µl) using KAPA HiFi HotStart DNA polymerase according to the manufacturer’s protocol, with the common forward and reverse primers with a 2 µl oligo pool. The individual oligo subpools were amplified analogously with subpool-specific primers (designed with LibGENiE19) and 1 µl purified oligo pool amplification product (~50 ng). The flanking upstream and downstream regions of the fragments were amplified using the parent variant as a template and corresponding primers (designed with LibGENiE19) in combination with T7_fw (TAATACGACTCACTATAGGG) or T7_rv (TGCTAGTTATTGCTCAGCGG).Full genes were reassembled by an overlap extension PCR (25 µl) with KAPA HiFi HotStart DNA polymerase according to the manufacturer’s protocol. In total, 1 µl amplified oligo subpool and 1 µl upstream and downstream PCR were used as template fragments with the T7_fw and T7_rv primers. The protocol included eight amplification cycles without flanking T7 primers, followed by 25 cycles with T7 primers.The reassembled genes were introduced into the pET22b vector containing the pelB signal tag using the MEGAWHOP52 cloning strategy. The PCR mixture (50 µl) contained 10 µl 5× Q5 reaction buffer, 2 µl dNTPs (10 mM each), 6 µl purified overlap extension PCR product (50–150 ng µl−1), 1 µl Q5 high-fidelity DNA polymerase and 1 µl template plasmid (100 ng, pET22b-pelB containing the parent gene of the evolution round). The template DNA was digested by incubation with DpnI (37 °C for 2 h).Generation of combinatorial librariesThe cloning procedure is illustrated in Extended Data Fig. 2.Primers encoding beneficial mutations and the parent amino acid were ordered from Microsynth. If beneficial mutations had appeared at the same site or in close proximity, they were encoded combinatorically on primers spanning the corresponding region using degenerate codons if feasible. Primers for such regions were mixed in equimolar quantities, ensuring an even distribution of the mutations.Gene fragments spanning primer regions were produced with complementary overlaps. The PCR (50 µl) contained 10 µl 5× Q5 reaction buffer, 1 µl dNTPs (10 mM each), 1 µl template plasmid (50–100 ng µl−1 of pET22b-pelB containing the parent gene for the evolution round), 2.5 µl combinatorial primer mix (10 µM), 2.5 µl complementary overlap primer and 1 µl Q5 high-fidelity DNA polymerase.The genes were reassembled by an overlap extension PCR (50 µl) with Q5 high-fidelity DNA polymerase according to the manufacturer’s protocol. The reaction contained 1 µl of each purified fragment as a template and 2.5 µl of the T7_fw and T7_rv primers, respectively. The temperature protocol included eight amplification cycles without flanking T7 primers, followed by 25 cycles including the flanking T7 primers.The PCR product was cloned in a pET22b-pelB vector using restriction enzymes XbaI and XhoI followed by a ligation step with T4 ligase. All procedures were carried out according to the manufacturer’s protocols.Generation of shuffled librariesThe genes of HG3.17 and HG3.R5 were amplified from the pET22b-pelB plasmid with T7_fw and T7_rv primers using the Q5 high-fidelity DNA polymerase according to the manufacturer’s protocol. The purified PCR products were mixed in equal quantities (500 ng each) and digested with 5 mU of DNAseI in 50 µl buffer (20 mM MnCl2 and 100 mM Tris (pH 7.4)) for 2, 3 and 5 min at 15 °C. DNAseI was inactivated (80 °C for 10 min), and the reactions were purified. In total, 4 ng of digested gene fragments were reassembled with a step-down PCR (50 µl) without primers using the Q5 high-fidelity DNA polymerase according to the manufacturer’s protocol. Reassembled genes were introduced into a pET22b vector containing the pelB signal tag with MEGAWHOP52 cloning as described earlier.Screening of the single-site variant, hit combination and shuffled librariesChemically competent E. coli NEB10 cells were transformed with 5 µl of cloning products. Colonies were grown on LB agar plates (100 mg l−1 ampicillin). All colonies were scraped off the plate, and plasmids were isolated with the NucleoSpin plasmid kit. Chemically competent E. coli BL21 (DE3) cells were transformed with 5 µl of the isolated plasmid and grown on LB agar plates (100 mg l−1 ampicillin). Precultures were prepared by inoculating 1 ml LB medium (100 mg ml−1 ampicillin) in 96-well deep-well plates (DWP) with single E. coli BL21 colonies and grown overnight at 30 °C with 300 rpm in a Duetz-System (Adolf Kühner AG). In total, 2 µl of the preculture was used to inoculate 1 ml of ZYM-5052 autoinduction medium (100 mg ml−1 ampicillin) in a 96-well DWP. The main cultures were grown for 15 h at 30 °C and 24 h at 20 °C at 300 rpm in a Duetz-System.The activity assay was performed in 96-well microtiter plates. Main cultures were diluted into reaction buffer (50 mM sodium phosphate (pH 7) and 100 mM NaCl) to obtain 100 µl of final volume. The ratio of culture to reaction buffer was adjusted according to the evolution round (R1: 20:80; R2: 20:80; R3: 10:90; R4: 1:99; R5: 0.5:99.5). The substrate (25 mM 5-nitrobenzisoxazole (1) in acetonitrile) was diluted into the reaction buffer (1:50), and 100 µl of this solution was immediately transferred to the diluted main cultures, initiating the reaction in a total volume of 200 µl containing 250 µM 5-nitrobenzisoxazole (1), 1% acetonitrile and round dependent volume of grown cultures (R1, 10 µl; R2, 10 µl; R3, 5 µl; R4, 0.5 µl and R5 0.25 µl). The reaction was monitored at 380 nm and 35 °C on a Tecan Spark (Tecan Group).Production and purification of the HG3 variantsGenes were cloned into the pET28b plasmid between NcoI/XhoI restriction sites (without the pelB-leader sequence but with the C-terminal His6-tag and stop codon). A preculture of 6 ml LB medium (50 µg ml−1 kanamycin) was grown overnight at 37 °C. In total, 500 ml of LB medium (50 µg ml−1 kanamycin) was inoculated with 5 ml preculture and incubated at 37 °C. Production of the HG3 protein was induced when the optical density at 600 nm (OD600) reached a value of 0.6 by adding 250 µM isopropyl-β-d-thiogalactopyranoside. The culture was then incubated overnight at 18 °C.Cells were collected by centrifugation for 60 min at 3,700g, 4°C, resuspended in 15 ml buffer (50 mM Tris–HCl (pH 7.4) and 500 mM NaCl), lysed by ultrasonic treatment (amplitude, 50%; pulse, 1 s/1 s and time, 1 min; Sonoplus) and clarified by centrifugation (1 h at 21,000g) as well as filtration (0.45 µm). The lysate was purified via nickel affinity chromatography on an ÄKTA Pure FPLC system (GE Healthcare) using a 5 ml HisTrap FF crude column (Cytiva). Loading buffer consisted of 500 mM NaCl and 20 mM imidazole in 50 mM Tris–HCl (pH 7.4), while the elution buffer contained 500 mM NaCl and 300 mM imidazole in 50 mM Tris–HCl (pH 7.4). Samples were desalted with three 5 ml HiTrap desalting columns (Cytiva) into a buffer containing 5 mM sodium phosphate (pH 6.0) and 20 mM NaCl and concentrated by ultrafiltration (Amicon Ultra-4 10K; Merck Millipore).Progress curvesProgress curves to obtain kinetic parameters were recorded according to the protocol given in ref. 35. The 10× concentrated Kemp eliminase variant stocks were prepared in reaction buffer (55.5 mM sodium phosphate (pH 6.8) and 111.1 mM NaCl), and 10× concentrated stocks of 5-nitrobenzisoxazole (1) were prepared in methanol. In total, 120 μl of reaction buffer and 15 μl of the enzyme stock solution were added to a quartz cuvette, and the absorption at 440 nm in a Cary 60 UV–Vis spectrophotometer (Agilent Technologies) was set to zero. The assay was initiated by the addition of 15 μl 10× concentrated substrate stock. The final enzyme concentration in the reaction buffer varied according to the variant (HG3, 19.0 μM; HG3.17, 61.2 mM; HG3.R1, 9.02 μM; HG3.R2, 2.23 μM; HG3.R3, 1.33 μM; HG3.R4, 225 nM; HG3.R5, 126 nM; HG3.R5w17, 3.88 μM and HG3.17wR5, 3.70 μM).Progress curves were measured by observing product formation at 440 nm (25 °C), applying an extinction coefficient of 1,050 M−1 cm−1. The initial substrate concentration was evaluated by diluting the reaction in an alkaline solution and measuring the product absorbance at 380 nm, applying an extinction coefficient of 15,800 M−1 cm−1 (15 μl reaction mixture was diluted in 235 μl of a 0.68 M sodium hydroxide solution). The data were analyzed using KinTek Explorer (6.1)53,54 to extract the steady-state parameters (Supplementary Figs. 1a and 2a).The rate constant for the background reaction (k1) was fixed to the value reported in the literature (k1 = 1.04 × 10−4)35. The reaction rate constants for HG3 and HG3.17 were corrected for a conformational selection step as reported in ref. 35 (k2HG3 = 2.32 × 10−4 s−1, k−2HG3 = 8.12 × 10−5 s−1, k2HG3.17 = 16.7 × 10−4 s−1 and k−2HG3.17 = 6.68 × 10−5 s−1). These reaction rates (k2 and k−2) were not available for other variants, which were assumed to be fully active as isolated to avoid overestimation of their activity. Furthermore, the on-rates for substrate and product (that is, k3 and k−5, respectively) were fixed to a diffusion-limited value of 1,000 μM−1 s−1. This left three parameters for fitting (that is, k−3, k4 and k5). The individual microscopic rate constants cannot be determined by fitting, but the values for kcat and Km,substrate can be estimated reliably using the equations below:$${k}_{{\mathrm{cat}}}=\frac{\,{k}_{4}{k}_{5}}{{k}_{4}+{k}_{-5}+{k}_{5}}=\frac{\,{k}_{4}{k}_{5}}{{k}_{4}+{k}_{5}}$$
(1)
$${K}_{{\mathrm{m}},{\mathrm{substrate}}}=\frac{\,{k}_{4}{k}_{5}+{k}_{-3}({k}_{-4}+{k}_{5})}{\left(\frac{{k}_{2}{k}_{3}}{{k}_{-2}+{k}_{2}}\right)({k}_{4}+{k}_{-4}+{k}_{5})}=\frac{\,{k}_{4}{k}_{5}+{k}_{-3}{k}_{5}}{\left(\frac{{k}_{2}{k}_{3}}{{k}_{-2}+{k}_{2}}\right)({k}_{4}+{k}_{5})}$$
(2)
$${K}_{{\mathrm{m}},{\mathrm{substrate}}}=\frac{\,{k}_{4}{k}_{5}+{k}_{-3}(\,{k}_{-4}+{k}_{5})}{\,{k}_{3}({k}_{4}+{k}_{-4}+{k}_{5})}=\frac{\,{k}_{4}{k}_{5}+{k}_{-3}{k}_{5}}{\,{k}_{3}({k}_{4}+{k}_{5})}$$
(3)
Initial rate analysisThe kcat/Km values shown in Table 1 were determined using initial rates analysis for product formation. Reactions were initiated by adding enzyme (HG3, 750 nM; HG3.R1, 459 nM; HG3.R2, 112 nM; HG3.R3, 65.7 nM; HG3.R4, 25.1 nM; HG3.R5, 6.6 nM; HG3.17, 6.8 nM; HG3.17wR5, 209 nM; HG3.R5w17, 208 nM) to the benzisoxazole substrate (31.5 µM to 2.0 mM final concentration) in 10% methanol, 100 mM NaCl and 50 mM sodium phosphate buffer (pH 7). Product formation was monitored at 380 nm with ε380 (pH 7) = 15,784.9 M−1 cm−1 (\({\varepsilon }_{380}={\varepsilon }_{\max }/\left(1+{10}^{\left({{\rm{p}}K}_{a}-{\rm{pH}}\right)}\right)\); εmax = 15,800 M−1 cm−1) in a Cary 60 ultraviolet/visible spectrometer (Agilent Technologies) at 27 °C using glass cuvettes (1 cm path length). Data were fitted to the linear portion of the Michaelis–Menten model (v0 = (kcat/Km)[E0][S]), and kcat/Km was deduced from the slope. The substrate stock solution concentration was determined by preparing a 2000-fold dilution in 10% methanol containing 10 mM sodium hydroxide and measuring the sample absorption at 380 nm with ε380 (pH > 9) = 15,800 M−1 cm−1.Melting temperature determinationThe melting temperature of the purified variants was determined with a Prometheus Panta (NanoTemper Technologies). The proteins were prepared at approximately 1 mg ml−1 in 5 mM sodium phosphate buffer (pH 6) containing 20 mM NaCl. The thermal unfolding experiment was performed from 20 to 95 °C with a ramp of 1 °C min−1. The derived melting profile was analyzed using the accompanying software (PR.Panta Analysis v1.6.3).Crystallization of HG3.R5 with the TSA 6-nitrobenzotriazole (3)HG3.R5 was produced and purified as described above, but the desalting step was performed on a size-exclusion column (HiLoad 16/600 Superdex 75 prep grade) with a buffer containing 5 mM sodium phosphate (pH 6) and 20 mM NaCl. The fractions containing protein were concentrated by ultrafiltration (Amicon Ultra-15 10K; Merck Millipore). The protein solution was prepared at a concentration of 30 mg ml−1 containing 5 mM 6-nitrobenzotriazole (3) (from a 100 mM stock solution in DMSO), resulting in a solution containing 5% DMSO, 4.75 mM sodium phosphate at pH 6 and 19 mM NaCl. Crystallization was performed using the sitting drop vapor diffusion method at 20 °C in Intelli-Plates 96-3 LVR (Hampton Research) at the Protein Crystallization Center (University of Zurich, Switzerland). A 150 nl drop of enzyme and substrate solution (30 mg ml−1 enzyme with 5 mM 6-nitrobenzotriazole (3)) was mixed with 150 nl precipitants (100 mM Tris–acetate (pH 7.4) and 1.3 M (NH4)2SO4) and equilibrated against a 50 µl reservoir solution (100 mM Tris–acetate (pH 7.4) and 1.3 M (NH4)2SO4). Crystal formation was observed after 25 days of incubation, and crystals were picked after incubation for a total of 41 days. In total, 2 µl of 30% glycerol in a reservoir solution was added to the drop as a cryoprotectant. The crystals were snap-frozen and stored in liquid nitrogen.Diffraction data were collected at the Swiss Light Source (Paul Scherrer Institute, Switzerland) on the beamline X06SA-PX at a temperature of 100 K and wavelength of 1.0000 Å. The data were (1) indexed and integrated with XDS (version 10 January 2022), (2) scaled and merged with Aimless (v0.7.9) and (3) solved using the experimental Kemp structure PDB code 5RGA (chain A) by molecular replacement using MOLREP (v11.0)55. Iterative refinement and manual model building steps were performed with REFMAC (v5.8.0419)56 and Coot (v0.9.8.92)57, respectively, embedded in the CCP4i2 (v8.0.013)58 suite. The final model with a resolution of 1.5 Å and a MolProbity59 score of 1.13 contained two protein chains (A and B) including the TSA (3).Library design based on evolutionary informationAn initial multiple sequence alignment (MSA) was created with the online HHblits tool60 (https://toolkit.tuebingen.mpg.de/tools/hhblits), using the UniRef30_2022_02 database and default parameters. The MSA was further processed with HHfilter61 (https://toolkit.tuebingen.mpg.de/tools/hhfilter) with the following settings: max ident: 90; min seq ident: 30; rest default. Variants were selected based on the consensus and frequency strategy outlined in the HotSpot Wizard overview30, leading to 76 additional variants included in each single-site library.Library design based on computational mutational scanningΔΔG predictions were based on a protocol published by the official Rosetta forums (https://www.rosettacommons.org/node/11126, 10 February 2021), and the relevant Python script can be found in the Supplementary Information. Each variant was predicted three times, and the lowest energy obtained was compared to the wild-type energies to calculate differences in free energy (Python version 3.8, PyRosetta version available at PyRosetta4.Release.python38.linux r275 2021.07+release.c48be26 c48be2695c4ba637c6fa19ee5d289fd9a8aa99ef).Library design based on tunnel and ligand analysisEnzyme tunnel analysis was performed with CAVER 3.0 (ref. 49) using the following parameters: probe_radius 0.9, shell_radius 4.0, shell_depth 5.0, frame_weighting_coefficient 1.0 and frame_clustering_threshold 1.0. The TSA (3) was introduced by aligning the crystal structure of HG3.17 (PDB 4BS0) to an AlphaFold50 model of each round’s parent variant. Sites for randomization were selected based on distances to the ligand and tunnels.Fitness landscapeThe protein sequences were encoded using the esm.pretrained.esm2_t48_15B_UR50D() model39. After extracting each variant’s embedding, these were reduced in dimensionality through principal component analysis (python package: sklearn, default settings, n_components = 2). In the absence of experimental data, the space outside the measured variants was set to 0. Within the area of measured variants, we performed radial basis function interpolation (python package: scipy, smooth = 0.1) and applied a further filter (python package: scipy.ndimage.uniform_filter, size = 2, mode = ‘constant’) to reduce ruggedness.MD simulation of HG3 variantsThe crystal structures of HG3.R5 (PDB 8RD5) were used for the homology modeling of HG3.R1, HG3.R2, HG3.R3 and HG3.R4. Simulations were also carried out with HG3 (PDB 5RGA) and HG3.17 (PDB 5RGE).The cocrystallized TSA (3) was manually replaced in all variants by the substrate 5-nitrobenzisoxazole (1). The PARSE force field was used at a pH 7, and protonation states for the titratable amino acids of all variants were determined using PROPKA3 (ref. 62), except for catalytic base D127 which was set up deprotonated.Substrate parameters were generated using Ambertools (23.3) packages antechamber and parmchk2 (ref. 63) with bond charge correction (AM1-BCC). The resulting mol2 and frcmod files were converted into xml files using ParmEd (4.1.0). The simulations were performed using OpenMM (8.1.0)64 with the ff14SB force field65. The rectangular box with a padding of 1.0 nm was solvated with water (TIP3P water model), and the total charge was neutralized. Energy minimization was conducted until 10 kJ mol−1 tolerance energy was reached. A temperature of 308.15 K and a pH of 7 were set, and the Langevin integrator was used with a friction coefficient of 1 ps−1 and a step size of 2.The long-range Coulomb interactions were computed by the particle mesh Ewald method, with a cutoff of 1 nm for the direct space interactions. The solvent was equilibrated by 5 ns NVT equilibration (using a Langevin integrator for temperature coupling) followed by 5 ns pressure coupled NPT equilibration using a Monte Carlo barostat at 1 atm pressure and 308.15 K reference temperature. The protein backbone and the substrate were restrained with a force of 100 kJ mol−1 Å−2. Finally, all restraints were removed, and a free equilibration of 5 ns was performed. The resulting system was used to produce five replicates of 100 ns each, with periodic boundary conditions. Every 1,000 steps, the trajectories were recorded.Analysis was conducted using MDTraj66. RMSF of the backbone Cα was calculated by superposing all frames onto the first frame as a reference.Contact frequency between protein and 5-nitrobenzisoxazole (1) atoms was calculated with the ContactTrajectory function of the MDTraj (1.9.9)66 derived Contact Map Explorer package. All substrate atoms were used as queries, while the full protein served as a haystack with a cutoff distance of 3 Å. Replicates in which the substrate left the binding site during simulation were not considered for this analysis.QM calculationsRestricted geometry optimizations and transition state searches were carried out with Gaussian 16 (ref. 67) on cluster models. Cluster models were generated from chain A of the crystallographic structure of HG3.R5 (PDB 8RD5) in complex with the TSA (3) and replaced with fully optimized 5-nitrobenzisoxazole (1). Cluster models include the substrate and the following residues: W44, P45, N47 (backbone only), S48 (backbone only), L49 (backbone only), Q50, G82, and G83 and D127 (side chain up to Cβ). All C atoms except substrate positions 3, 3a and 7a were kept fixed at the crystallographic coordinates. Calculations were carried out using the M06-2X68 hybrid functional and 6-31G(d) basis set with ultrafine integration grids. Solvent effects in water were considered implicitly through the IEF-PCM polarizable continuum model69. The energies of the calculated structures are summarized in Supplementary Table 2.QM/MM calculationsThe initial geometry for the HG3.R5–5-nitrobenzisoxazole (1) complex was generated from chain A of the crystallographic structure (PDB 8RD5) of HG3.R5 in complex with the TSA (3) and replacing its coordinates with the optimized coordinates of 5-nitrobenzisoxazole (1) (M06-2X/6-31G(d) level of theory). Crystallographic waters were preserved. The enzyme–substrate complex was then prepared for classical MD relaxation with the Amber 22 (ref. [63) suite using the ff19SB70 force field for the protein and gaff2 (ref. 71) for the substrate. The complex was immersed in a water box with a 10 Å buffer of OPC72 water. A two-stage geometry optimization approach was implemented. The first stage minimizes only the positions of solvent molecules, and the second stage is an unrestrained minimization of all the atoms in the simulation cell. The systems were then heated for 100 ps by incrementing the temperature from 0 to 300 K under a constant pressure of 1 atm and periodic boundary conditions with a timestep of 1 fs. Harmonic restraints of 50 kcal mol−1 Å−2 were applied to the solute, with the Andersen73 temperature coupling scheme to control and equalize the temperature. The last geometry of the equilibration trajectory was extracted, and solvent molecules were trimmed, leaving an 8 Å thick water shell around the solute.Full geometry optimizations and transition state searches were carried out with Gaussian 16 (ref. 67; ONIOM74 QM/MM method, with electrostatic embedding) using the M06-2X hybrid functional and 6-31+G(d) basis set with ultrafine integration grids for the QM part. The MM part was described with the ff14SB65, gaff2 (ref. 71) and TIP3P75 force fields for the protein, ligand and water, respectively. The QM layer includes the whole 5-nitrobenzisoxazole (1), the side chains (up to Cβ) of D127 and Q50 and the two conserved water molecules shown in Extended Data Fig. 4. Only solvent molecules were kept frozen, except for the two water molecules in the QM layer, to allow adaptation of the enzyme along the reaction path. All stationary points were characterized by a frequency analysis performed at the same level of theory, from which thermal corrections were obtained at 298.15 K. ONIOM energies, entropies, enthalpies, Gibbs free energies and lowest frequencies of the calculated structures are summarized in Supplementary Table 3. The effect of crystallographic water we on the activation barrier was computed by manually removing we and reoptimizing the reactant and transition state at the same level of theory.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Enriching productive mutational paths accelerates enzyme evolution
