In silico screening of 393 mutants facilitates enzyme engineering of amidase activity in CalB

Our previously presented method for high throughput computational screening of mutant activity (Hediger et al., 2012) is benchmarked against experimentally measured amidase activity for 22 mutants of Candida antarctica lipase B (CalB). Using an appropriate cutoff criterion for the computed barriers, the qualitative activity of 15 out of 22 mutants is correctly predicted. The method identifies four of the six most active mutants with ≥3-fold wild type activity and seven out of the eight least active mutants with ≤0.5-fold wild type activity. The method is further used to screen all sterically possible (386) double-, triple- and quadruple-mutants constructed from the most active single mutants. Based on the benchmark test at least 20 new promising mutants are identified.


1
Introduction 2 Numerous methods are currently being proposed and developed for the description of enzyme activities, the theoretical background of which ranges from phenomenological and bioinformatics based approaches [9][10][11][12][13] to quantum mechanics based ab initio descriptions [14][15][16][17][18][19][20][21][22]. However one can expect that methods which are highly demanding in terms of set-up efforts and computational time are less likely to be employed in industrial contexts where qualitative or semi-quantitative conclusions can be of sufficient use in the beginning and planning phase of a wet-lab study. Few approaches, while taking into account a number of approximations and limitations in accuracy, aim at being used in parallel or prior to experimental work [23,24] and are not designed to be used for high throughput fashion.
Hediger et al. have recently published a computational method for high throughput computational screening of mutant activity [1] and in this paper we benchmark the method against experimentally measured amidase activity for mutants of Candida antarctica lipase B (CalB) and apply the method to identify additional promising mutants.

Methods
We introduce the experimental set-up and the methodology for comparing experimental and computational data. We describe a benchmarking and a combinatorial study of CalB mutant activity.
Experimentally, variants of Candida Antarctica lipase B (CalB) were either produced in Pichia pastoris with C-terminal His6-tag for subsequent affinity purification or expressed in Aspergillus oryzae without terminal tag followed by a three-step purification procedure.
It is generally accepted that in serine protease like enzymes, the formation of the tetrahedral intermediate (TI, Fig. 1) is rate determining [25][26][27][28] and throughout this work we assume that a lower barrier for this reaction correlates to increased overall activity of the enzyme. 3 remove parent templates, they were methylated in vitro prior to PCR with CpG methyltransferase (from NEB) and digested in vivo after transformation of competent E.coli DH5 α cells (TaKaRa) according to the instructions from the manufacturer. Plasmid DNA was isolated from transformed E.coli strains, and sequenced to verify the presence of the desired substitutions. Confirmed plasmid variants were used to transform an Aspergillus oryzae strain that is negative in pyrG (orotidine-5'-phosphate decarboxylase), proteases pepC (aserine protease homologous to yscB), alp (an alkaline protease) NpI (a neutral metalloprotease I) to avoid degradation of the lipase variants during and after fermentation.
The transformed Aspergillus strains were fermented as submerged culture in shake flasks and the lipase variants secreted into the fermentation medium. After the fermentation, the lipase variants were purified from the sterile filtered fermentation medium in a 3 step procedure with 1) hydrophobic interaction chromatography on decylamine-agarose, 2) buffer exchange by gel filtration and 3) ion exchange chromatography with cation exchange on SP-sepharose at pH 4.5. The lipase variant solutions were stored frozen.

Generation of CalB Variants with His-tags
Variants of CalB carrying the CalB signal peptide and C-terminal His-tags were generated at the DNA level using SOE-PCR and inserted into a dual E.coli /Pichia pastoris expression vector using In-fusion cloning (ClonTech). The SOE-PCR was performed with Phusion DNA polymerase (NEB) and template DNA of the CalB gene. The cloned plasmids were transformed in competent E.coli DH5 α cells (TaKaRa).
Plasmid DNA was isolated from transformed E.coli strains, and sequenced to verify the presence of the desired substitutions. Confirmed plasmid variants were used to transform a Pichia pastoris strain that is Mut(s), Suc(+), His(-). The transformed Pichia strains were fermented as submerged culture in deep well plates and secretion of the lipase variants into the fermentation medium was induced by addition of methanol. After the fermentation, the lipase variants were purified from the cleared supernants using a 4 (THF or DMSO). Reactions containing 5 mM amide substrate, 0.3-3 µM enzyme, and 12 µg/mL BSA were incubated for 18-20 h at 37 • C in a shaker incubator. In a second step, 50 µL of a 20 mM 4nitro-7-chloro-benzo-2-oxa-1,3-diazole (NBD-Cl) solution in 1-hexanol was added and the reaction of NBD-Cl with benzylamine formed during amide hydrolysis proceeded under identical reaction conditions for another hour.
Fluorescence of the final reaction product was determined with excitation at 485 nm and measured emission at 538 nm. Calibration of the amide hydrolysis reaction was performed on each assay plate with benzylamine covering a concentration range between 0.05 and 5 mM. All enzymatic activities were corrected for non-enzymatic background reaction determined under identical conditions without enzyme present.

Computational Details
The computational method used to estimate the reaction barriers of the CalB mutants has been described in detail earlier [1] and is only summarized here.
As described previously [1], in order to make the method computationally feasible, relatively approximate treatments of the wave function, structural model, dynamics and reaction path are used. Given this and the automated setup of calculations, some inaccurate results will be unavoidable. However, the intent of the method is similar to experimental high throughput screens of enzyme activity where, for example, negative results may result from issues unrelated to intrinsic activity of the enzyme such as imperfections in the activity assay, low expression yield, protein aggregation, etc. Just like its experimental counterpart our technique is intended to identify potentially interesting mutants for further study. Reviewing Manuscript 5 being the number of interpolation frames and i the interpolation frame index). In geometry optimization calculations, the gradient convergence criteria is set to 0.5 kcal/(molÅ) and a linear scaling implementation of the PM6 method (MOZYME [31]) together with a NDDO cutoff of 15Å is applied. The energy profile of the reaction barrier at the PM6 level of theory [32] is subsequently mapped out by carrying out conventional SCF calculations of each optimized interpolation frame. All calculations are carried out using the MOPAC suite of programs [33,34]. The molecular models are based on the crystal structure of the CalB enzyme with PDB identifier 1LBS [35]. In order to prevent significant rearrangement of hydrogen bonding network of surface residues during the optimization, a number of additional structural constraints are applied in the geometry optimizations, i.e. the residues S50, P133, Q156, L277 and P280 are kept fixed. These (surface) residues are observed to rearrange and form new hydrogen bonds in optimizations when no constraints are applied. Omitting the constraints leads to unconclusive barrier shapes containing many irregular minima along the reaction coordinate which do not permit to readily define a reaction barrier.
For the analysis, the reaction barrier is defined by the difference between the highest energy point on the reaction profile and the energy corresponding to the enzyme substrate complex. From our calculations (PM6//MOZYME in vacuum), we estimate the wild type (WT) barrier to be 7.5 kcal/mol. Experimentally, specific activity of hydrolysis is determined. Given first order kinetics, saturation of the enzyme with substrate (usual for industrial application) and fast binding and product release, the catalytic rate constant k cat is directly proportional to the specific activity under the assumption that the amount of active enzyme remains constant. This therefore allows the catalytic rate constant k cat and, hence, the barrier height to be compared to the improvement factors reported in the results section. The approximations used here in relating the barrier height on the potential energy surface to k cat have been discussed previously [1].
It is noted that using one CPU per interpolation frame on the reaction barrier, the complete barrier of one mutant can be computed with 10 CPUs usually within less than 12 hours of wall clock time (for a molecular model of the size used in this study). Given a set of molecular models of the enzyme, and 100 available CPUs, it is possible to screen around 1000 mutants within one week.

Combination Mutants
The molecular model of the enzyme and the positions of the point mutations in the enzyme are illustrated in Fig. 2. The point mutations are listed in Table 2. Two sets of mutants are introduced in this section. The point mutations are selected based on different design principles. These are either introduction of structural rearrangements in the active site to change the binding site properties of the active site (residues P38, G39, G41, T42, T103) [2], introduction of space to accomodate the substrate (W104, L278, A282, I285, V286), introduction of dipolar interactions between the enzyme and the substrate (A132, A141, I189) [36] or reduction of polarity in the active site (D223). Of course different heuristic considerations will apply for other enzymes when selecting the single mutations for combinatorial study. The mutants of the benchmarking study are collected in a small set S (22 mutants, Table 1). For the combinatorial study, out of the above we select six residues (G39, T103, W104, A141, I189, L278) which, it is assumed, contribute strongest to increased activity and define the mutations at each position as listed in Table 3.
Given the position i and the number of mutations at each position g i , in general the upper limit for the number of mutants M in a combinatorial study can be calculated by writing a sum term for each type (i.e. "order ") of combination mutant, i.e. single, double, . . . , such that where each sum term consists of N o individual terms (N and o being the number of positions which can be mutated and the order of the mutant, respectively). By this scheme, considering the mutations listed in Table 3, hypothetically 424 (= 13 + 64 + 154 + 193) single to four-fold mutants can be constructed. This number is reduced by applying the restriction that out of the 424 hypothetically possible mutants, 0 single, 2 double, 12 triple and 24 four-fold combination mutants including the pair A141N/Q-I189Y are discarded because in the molecular modeling, these side chains could not be allocated spatially in the same mutant. We further note that 15 out of these remaining 386 mutants (Table 3) are present also in the benchmarking set S and thus the combinatorial study consists of 371 unique mutants. A detailed documentation of the number of screened residues in the combinatorial study is provided in Table 4.  so the most conservative choice is to consider it a non-promising candidate for a more active variant.
Generating plots of the profiles is completely automated and visual inspection can easily be done for hundreds of mutants. Furthermore, out of the mutants with regular reaction barrier shapes, we discard those mutants with barriers >19.0 kcal/mol (i.e. the largest calculated barrier from set S). Following these selection criteria, 61 mutants are discarded because of inconclusive barrier shapes and 47 mutants because the barrier is higher than 19 kcal/mol (a distribution of reaction barriers is shown in Fig. S1).
After these filtering steps, 278 mutants remain in the combinatorial study which we collect in the large set L (out of which 15 are in set S). An overview on the distribution of reaction barriers for the mutants from set L is provided Fig. S2 of the supporting information.
We note that in set S, all barriers appear regular in shape and no mutant contains the A141N/Q and I189Y pair.

Results and Discussion
Set S: Calibration of the accuracy The correspondance of the computed barriers from set S with the experimental assay is shown in Fig. 3.
The exact data is reported in Table 1. A scatterplot of calculated reaction barriers is presented in Fig.

4.
We note that in set S, the highest experimentally observed activity is around 11 times the wild type activity (G39A-T103G-W104F-L278A, Table 1), while roughly ten mutants show no increased activity.
In total, six mutants show 3-fold or higher wild type activity. In the calculations, only one mutant is observed to have a lower barrier than the wild type (7.3 kcal/mol, G39A-T103G-L278A) and the highest observed barrier is 18.9 kcal/mol (I189G).
Given the approximations introduced to make the method sufficiently efficient, it is noted that the intent of the method is not a quantitative ranking of the reaction barriers, but to identify promising mutants for, and to eliminate non-promising mutants from, experimental consideration. Therefore only qualitative changes in overall activity are considered, which are represented by the activity factors (+1/-1).
We categorize the experimentally observed activities and the predicted reaction barriers as follows. Reviewing Manuscript 8 as improving (degrading). Correspondingly, the computed difference in reaction barrier height between a mutant and the wild type is expressed in qualitative terms. For the comparison with the experimental activity assay, we define a barrier cutoff c S = 12.5 kcal/mol to distinguish between potentially improving and degrading mutants in set S. The value of 12.5 kcal/mol is chosen such as to maximize the agreement with experiment, which is 68%, i.e. using a smaller or larger value for the cutoff will decrease this value.
A mutant with a predicted barrier ≥c S (12.5 kcal/mol) is considered to likely have decreased activity compared to the wild type while mutants with reaction barriers <c S are considered likely having increased activity.
We note that defining the cutoff is done purely for a post hoc comparison of experimental and computed data. When using the computed barriers to identify promising experimental mutants, one simply chooses the N mutants with the lowest barriers, where N is the number of mutants affordable to do experimentally (e.g. 20 in the discussion of set L).
Based on this approach, qualitative activity of 15 out of 22 mutants is correctly predicted. It is noted that the correlation is best for mutants with largest activity difference compared to wild type (both positive or negative). For example the method identifies four of the six most active mutants with ≥3-fold wild type activity. Similarly, the method identifies seven out of the eight least active mutants with ≤0.5-fold wild type activity. For mutants with only small differences in activity compared to wild type, the predictions are less accurate.
Set L: Large scale screening study Set L is screened to identify new mutants for which increased activity is predicted. The 20 mutants with the lowest barriers are suggested as candidates for further experimental study in Table 5 In set L, three new mutants are identified with barriers lower than the predicted wild type barrier. Out of the 20 mutants suggested in Table 5, three are double mutants, seven are three-fold and ten are four-fold mutants. No single mutants where found for which increased activity compared to wild type is predicted.
All mutants except one contain the G39A mutation, five contain the T103G mutation, six contain a mutation of W104, 13 contain a mutation of A141, 16 contain a mutation of I189 and eight contain the L278A mutation. From this observation it is likely that mutations of G39, A141 and I189 will likely contribute to an increased activity of the mutant and should thus be included in future experimental activity assays. out of the sixty mutants with lowest barriers (Fig. S3), 33 contain a mutation of W104 out of which 17 are suggested to be W104F, while 14 are suggested to be W104Y (two contain W104Q).
The mutation of I189 is analysed in a similar way. In set L, five different mutations of this residue are screened ( Table 3). The single mutant with the lowest barrier is I189Y and the two mutants with the lowest predicted barrier contain this mutation as well (Table 5). Similarly to above, higher order mutants containing I189A, I189G, I189H or I189N are predicted to have considerably lower barriers than the corresponding single mutants, Fig. 5B. Particularly, out of the mutants listed in Table 5, three contain the I189A, one contains I189G mutation, four contain the I189H mutation and three contain the I189N mutation.
As a special case we highlight that the single mutant I189G has one of the highest calculated barriers (18.9 kcal/mol, Table 1), however the four-fold mutant G39A-A141Q-I189G-L278A has one of the lowest barriers (6.3 kcal/mol, Table 5). Interestingly, the mutant G39A-A141Q-L278A has an intermediate barrier (10.9 kcal/mol). It would appear that I189G as a single mutant is counterproductive (high computed barrier) but lowers the barrier of G39A-A141Q-L278A. This observation is further supported by the observation that the I189G mutation is in spatial proximity to A141Q. While it is difficult to quantify the interaction, it is likely that in the mutant, the rather large side chain of A141Q is better accommodated in the active site and can better interact with the substrate.
Observations as these should be kept in mind when selecting the single mutants to be considered when preparing higher order mutants.

Conclusions
Our previously presented method for high throughput computational screening of mutant activity [1] is benchmarked against experimentally measured amidase activity for 22 mutants of  Experimentally, amidase activity is successfully introduced in 12 mutants, the highest activity is determined to be 11.2 -fold over the wild type activity.
Using an appropriate cutoff criterion for the computed barriers, the qualitative activity of 15 out of 22 mutants is correctly predicted. It is noted that the correlation is best for mutants with largest activity difference compared to wild type (both positive and negative). For example the method identifies four of the six most active mutants with ≥3-fold wild type activity. Similarly, the method identifies seven out of the eight least active mutants with ≤0.5-fold wild type activity.
Thus validated, the computational method is used to screen all sterically possible (386) double-, tripleand quadrupole-mutants constructed from the most active single mutants. Based on the benchmark test at least 20 new promising mutants are identified.
These mutants have so far not been tested experimentally and are thus offered as scientifically testable predictions. Interestingly, we observe that single mutants that are predicted to have low activity appear to have high activity in combination with other mutants. This is illustrated in specific analysis of effects of mutations of two different positions (104 and 189).

NHis224
Ser105 H      Tables   Table 1. Experimental overall activities and calculated reaction barriers of Set S. Activity factors +1/-1 indicate increased/decreased overall activity. Ao and Pp indicating expression in organisms (Org.) Aspergillus oryzae or Pichia pastoris, respectively. The cutoff to distinguish higher and lower activity mutants is defined as 12.5 kcal/mol, see text.