Random Saturation Mutagenesis to Generate Highly Diverse Libraries for Directed Evolution

Random Saturation Mutagenesis to Generate Highly Diverse Libraries for Directed Evolution | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Random Saturation Mutagenesis to Generate Highly Diverse Libraries for Directed Evolution David Weitz, Xinge Zhang, Anqi Chen, Liyin Chen, Rohan Thakur, and 4 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-5583255/v1 This work is licensed under a CC BY 4.0 License Status: Under Review Version 1 posted You are reading this latest preprint version Abstract Directed evolution improves the function of enzymes by diversifying the wild-type to generate variants and searching for the desired functions. When prior knowledge about the wild-type is limited, random mutagenesis is employed to introduce mutations at arbitrary positions. However, conventional methods exhibit significant bias in the distribution of amino acid substitutions, resulting in inaccessible chemical properties, limiting the chemical diversity, and hindering discovery of substantially improved variants. To overcome such limitations, we develop Random Saturation Mutagenesis (RSM) that uniformly introduces all amino acid substitutions across the gene at controlled rates by using overlapping site-mutagenesis oligo pools. Sequencing results illustrate the high quality and diversity of RSM libraries as compared to conventional libraries. We demonstrate the utility of RSM by improving thermotolerance of LipaseA. With only one round of mutagenesis, we identify two variants with over 40°C improvement in thermotolerance, demonstrating the potential of RSM to transform and accelerate directed evolution. Biological sciences/Biochemistry/Biocatalysis Biological sciences/Chemical biology/Enzymes Biological sciences/Biological techniques/High-throughput screening Directed evolution Mutagenesis library generation Enzyme engineering Ultra-high throughput screening Synthetic biology Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Introduction Enzymes are biological catalysts that facilitate biochemical reactions with efficiency and specificity. However, their application in industrial processes often demands adaptations to maintain stability and catalytic performance under challenging conditions. These adaptations can be achieved through a process called directed evolution. Directed evolution comprises two steps: introducing mutations to create a mutagenesis library and screening for variants with desired functional improvements. For a typical-sized enzyme, the sequence space is vast, far exceeding the throughput of existing screening methods. Therefore, when a library is constructed, a subspace in the sequence space is determined, and the choice of the sequence subspace defines the range for the screen. Mutagenesis libraries are designed to maximize the chance of discovering improved variants. For well-studied enzymes, mutagenesis libraries can be designed rationally based on prior structural and mechanistic information. For enzymes lacking sufficient prior information, random mutagenesis is used. Random mutagenesis introduces mutations across the entire sequence, allowing exploration of a broader sequence space 1 . Traditionally, random mutagenesis is achieved using error-prone polymerases, either in vitro through polymerase chain reactions (PCRs) or in vivo through continuous evolution methods 2 , 3 . Error-prone polymerases are strongly biased towards single nucleotide changes per codon, resulting in a high frequency of silent and redundant mutations. Moreover, with a single nucleotide mutation per codon, only 5.7 amino acid substitutions are accessible, on average, to each wild-type amino acid 4 . Furthermore, even these substitutions are often physiochemically conserved due to the evolutionary organization of codons to preserve activity 5 . This narrows the chemical diversity to an evolutionarily favored subspace and results in limited deviations from the wild-type function; many amino acid substitutions leading to large functional improvements can be excluded from the library 6 – 10 . Sequence Saturation Mutagenesis (SeSaM) pioneers the efforts to improve chemical diversity in random mutagenesis by allowing for some consecutive nucleotide exchanges 11 – 14 . In SeSaM, DNA templates are randomly fragmented; the fragments are tailed with one or multiple consecutive universal bases, and then used as primers to synthesize full-length DNA. These products are further amplified by PCR where universal bases are replaced by one of the four standard bases, introducing mutations. This workflow is rather complex and time consuming, yet still approximately 70% of the mutations remain single nucleotide substitutions 14 , and only 35% of all amino acid substitutions 15 are enabled. To achieve greater chemical diversity at specific amino acid positions, site saturation mutagenesis (SSM) is often employed. This approach allows for the systematic introduction of all possible amino acid substitutions at specific sites. However, these sites are chosen either by the screening outcome of random mutagenesis libraries 16 , where important residues may be missed because they were not saturated 8 , 15 , or by mechanistic hypotheses based on prior knowledge, reducing the chance to discover unexpected mechanisms 17 . Therefore, there exists the need for a random mutagenesis method that can achieve high chemical diversity across the entire gene without requiring prior knowledge. A method that introduces all distinct amino acid substitutions randomly throughout an enzyme sequence, while controlling the mutation rate, would allow for a more comprehensive exploration of the sequence space. In this paper, we present Random Saturation Mutagenesis (RSM), a method that introduces mutations across the entire protein sequence, evenly sampling all chemically distinct conversions at a controlled mutation rate. This method uses multiple pools of overlapping mutagenic primers targeting subregions of the gene, allowing for the control of mutation rates by varying the number of primer pools used in each reaction. Sequencing results confirm that RSM libraries achieve controlled mutation rates, high variant quality, and high chemical diversity. To demonstrate the utility of our library preparation method, we use RSM to drastically improve the thermotolerance of Bacillus subtilis Lipase A (BsLipA). We construct an RSM mutagenesis library and use droplet screening to identify two vastly improved variants after only one generation of mutagenesis. Each of the variants has seven single amino acid substitutions and shows over 40°C increase in T 50 values, the temperature at which they retain 50% of their activity after heat incubation at that temperature. To understand the roles of the newly discovered amino acid substitutions in improving thermotolerance, we conduct molecular dynamics simulations. The identified key substitutions that result in large T 50 improvements are not obtainable by single nucleotide changes in a trinucleotide codon. Results Construction of RSM Libraries To demonstrate the utility of Random Saturation Mutagenesis (RSM) to create a random library that is highly chemically diverse, we choose the enzyme BsLipA as a candidate to improve its thermotolerance. Thermotolerance of enzymes can be improved by mutations scattered across the entire gene 18 ; therefore, random mutagenesis is particularly well-suited for this problem. When there is no prior information about the structure or function, an ideal random mutagenesis library should be as chemically diverse as possible to facilitate a comprehensive and unbiased exploration of the sequence subspace 6 – 10 . We use site-directed mutagenesis primers to introduce mutations. To ensure annealing of primers during mutagenesis, we choose a primer length of 30 nucleotides, covering a 10-residue region as shown in Fig. 1 a. We divide the 540 bp BsLipA gene into 18 regions, each targeted by a set of fully overlapping 30-nucleotide primers, as shown in 1b. For simplicity, the schematic illustrates only 5 primer pools, each covering 5 residues. For each residue, we design a set of 19 mutagenesis primers to introduce all possible single amino acid substitutions, saturating each residue. Each primer is complementary to the target gene except for the three nucleotides encoding the amino acid substitution. We select a single codon for each unique amino acid substitution to ensure even distribution at the mutated site. Moreover, the explicit choice of 19 codons prevents the introduction of rare codons and stop codons. The primers targeting all 10 residues in each region are purchased as a primer pool of 190 one-directional primers (10 positions x 19 substitutions per position). To cover the full gene, 18 such primer pools are purchased. To introduce mutations throughout the gene, we combine all 18 primer pools into a single mutagenesis reaction using the wild-type plasmid as the template. During a mutagenesis reaction, only one-directional primers are used. The polymerase amplifies the gene from the 3’ end of the first primer to the 5’end of the next primer leaving a nick in between. A ligase fixes this nick to form a double-stranded plasmid incorporating the designed mutations in one strand; the wild-type strand is then digested by DpnI. In this reaction, a maximum of 18 primers can anneal to the template DNA simultaneously. We transform the single-stranded DNA product from the mutagenesis reaction into bacterial hosts to synthesize double-stranded plasmids. At this stage, the number of mutations per region is limited to one, as shown in Fig. 1 c. To increase mutation diversity, we use staggered extension PCR (StEP) 19 . In a StEP reaction, the extension step is set deliberately to be short, synthesizing only partial genes per cycle. During the denaturing step of the next cycle, these incomplete strands unbind from their original templates and anneal to different ones, continuing extension and effectively recombining mutations from various templates. The result is a final RSM library with random mutations throughout the gene, as illustrated for 5 primer pools in Fig. 1 d. We sequence the library using PacBio technology to obtain high-quality full-length reads and analyze approximately 10,000 of these by aligning them to the wild-type BsLipA gene to identify mutations. On average, this RSM library contains 7.85 single amino acid substitutions per gene. The high mutation rate obtained in this library can result in a predominant population of inactive mutants, which may not be suitable for standard directed evolution applications. To achieve controlled mutation rates, we modify the RSM workflow by using only a subset of all primer pools for each mutagenesis reaction, and performing multiple separate reactions to cover the whole gene. This approach limits the number of primers that can anneal to the DNA template in a single reaction. For example, we can prepare two separate mutagenesis reactions, each containing an equimolar mix of 9 primer pools, as shown in Fig. 1 e. Reaction 1 uses all odd-numbered primer pools and reaction 2 uses all even-numbered primer pools. To confirm the mutations introduced in either reaction, we sequence the product from these two reactions and count the amino acid substitutions at each residue position. To show the distribution of substitutions at each position, we compile a heatmap of which each column represents a residue in the BsLipA sequence, and each row shows the mutations introduced. In both reaction 1 and 2, the mutations are confined within the boundaries of the designed primer pools as shown in Fig. 2 a and 2 b. To remove such confinement and to increase diversity of the final RSM library, the products from the two reactions are recombined through StEP. The final product of this library exhibits uniform coverage of single amino acid substitutions across the gene as shown in Fig. 2 c. The average number of single amino acid substitutions of this library is measured to be 4.64 by PacBio sequencing. Since we use 9 primer pools in 2 mutagenesis reactions, we label this library “9×2”, where 9 denotes the number of pools used per reaction, and 2 denotes the number of reactions. Similarly, the previous library is labeled “18×1”. Following the same procedure, we create a “6×3” library by using 6 primer pools in 3 separate reactions, and a “4×3 + 3×2” library by combining 3 or 4 primer pools per reaction. The average mutation rates for these libraries are measured by PacBio sequencing to be 5.09 and 1.97 amino acid substitutions per protein sequence, respectively. Moreover, the RSM approach can be adapted to construct single saturation mutagenesis libraries by using only one primer pool in 18 separate mutagenesis reactions, without the subsequent StEP recombination (Extended Data Fig. 1 ). The overlapping primer design ensures that the annealing of a primer completely prevents other primers in this region from binding. The average mutation rate for this single saturation library is confirmed to be 0.99 measured by PacBio sequencing, and we label this library as “Single”. The mutation distribution is shown in Fig. 3 d. By performing RSM with varying numbers of primer pools per reaction, mutations can be introduced at controlled rates and random positions throughout the gene. While existing methods like error-prone PCR can achieve similar controlled mutation rates across the gene, by comparison to RSM, they produce low-quality and low-diversity libraries. To illustrate this comparison, we construct three error-prone mutagenesis libraries using a commercially available mutazyme kit (GeneMorph II, Agilent), rendering libraries with high (EPhigh, 2.70 substitutions per protein sequence), medium (EPmid, 2.04), and low (EPlow, 1.18) mutation rates. For a variant in a library to be improved, it must translate to the full length and must differ from the wild-type protein sequence. Silent mutations do not contribute to the diversity of the library. All RSM libraries consistently exhibit less than ~ 5% silent mutations, significantly lower than the error-prone libraries, each of which contains ~ 30% silent mutations, as shown in Fig. 3 b. Additionally, variants that translate into truncated proteins cannot function, hence their presence also negatively impacts library quality. Such truncations can occur when a trinucleotide codon is mutated to either TAG, TGA, or TAA, which are stop codons signaling the termination of the translation process. In all RSM libraries, premature stop codons are found to exist in less than 1% of all variants, much lower than in error-prone libraries, where 5%-14% of variants contain premature stop codons, as shown in Fig. 3 a. Similarly, insertions and deletions (indels) of nucleotides can disrupt protein function by causing frameshifts that lead to incorrect and truncated amino acid sequences if they do not occur in multiples of three. When comparing an RSM library to an error-prone library with similar mutation rates, such as single versus EPlow and 4×3 + 3×2 versus EPmid, the proportion of variants that contain indels in RSM libraries is significantly lower than in error-prone libraries, as shown in Fig. 3 c. Moreover, in error-prone libraries, the proportion of variants with premature stop codons and with indels both increase with mutation rate. Thus, if a mutation rate even higher than EPhigh is required, the proportion of truncated variants will inevitably increase further. By contrast, in RSM libraries, the proportions of these incorrect and truncated variants are not only low but also independent of the mutation rate; this demonstrates the significant advantage of the RSM libraries over the error-prone libraries particularly when high mutation rates are required, as shown in Fig. 3 b and 3 c. Wild-type BsLipA contains 19 of the 20 standard amino acids (excluding cysteine), each of which can be mutated into 19 other possible amino acids, resulting in 361 possible distinct amino acid substitutions (DAASs). An ideal random mutagenesis method for our application of interest should uniformly sample DAASs to ensure high chemical diversity of the library. To evaluate the uniformity of sampling of these 361 DAASs by RSM and by error-prone PCR, we choose an RSM library, 9×2, and an error-prone library, EPhigh, with similar mutation rates and sample the same number of amino acid substitutions from each library. We count the occurrences for each DAAS and rank them by abundance to visualize the uniformity of sampling in each library. The 9×2 RSM library samples DAASs much more uniformly than the EPhigh library. In 9×2, all 361 DAASs are enabled with a mean of 21 and a standard deviation of 10. By contrast, in EPhigh, the top ~⅓ most abundant DAASs comprise over 99.5% of the total count, whereas the remaining ~⅔ DAASs each occur less than or equal to once, as shown in Fig. 3 d. While the mean of DAAS occurrence is the same as 9×2, the standard deviation is 38, four times that of 9×2. To evaluate the deviation of either library from an ideal uniform distribution, we calculate the Gini coefficient for the normalized counts of DAASs. The 9×2 RSM library exhibits a Gini coefficient of 0.27, lower than 0.78 for EPhigh, indicating a much more uniform distribution of DAASs in the RSM library (Extended Data Fig. 2 ). To determine the source of the bias in EPhigh, we investigate the codon change required for each DAAS. We arrange the counts of DAASs in a 20x19 table where each column represents an amino acid used in the wild-type BsLipA, and each row indicates an observed substitution from the wild-type amino acid. The color reflects the counts of individual DAASs observed in sequencing. The color-annotated table shows the much lower variation in the distribution of DAASs sampled in the RSM library as compared to the EPhigh library, as shown in Fig. 4 a,b. Moreover, in the EPhigh library, the top ⅓ most abundant substitutions, which account for > 99.5% of all occurrences, are precisely those enabled by single nucleotide changes per codon, as highlighted by the white frames in Fig. 4 a. Certain amino acid substitutions occur with anomalously high frequency, such as the W-to-R substitution and the F-to-L substitution, likely due to the intrinsic bias of the error-prone polymerase. By contrast, the DAASs outside of the white frames, which require multiple nucleotide changes, are only detected at most once; hence, they are effectively inaccessible by error-prone PCR. On average, each amino acid only has access to 6.6 DAASs, and these DAASs cannot enable all categories of sidechain physicochemical properties. To access the accessibility of each wild-type amino acid to different sidechain properties, we classify all 20 amino acids into eight categories: non-polar (NP, M, I, L, V, A), polar uncharged (PU, S, T, Q, N), positively charged (PC, R, K, L), negatively charged (NC, D, E), aromatic (Ar, F, T, Y), and three special groups, P, C, G, each consisting of only one amino acid due to their distinct properties. Through error-prone PCR, each wild-type amino acid has a few inaccessible property categories, as shown in Fig. 4 c. Furthermore, in this analysis, we aggregate the same wild-type amino acids at different positions which are encoded by different codons. Therefore, when each residue position is individually considered, the categorical accessibility of each wild-type residue is even more limited. In contrast, the 9x2 RSM library uniformly samples all possible DAASs, and thereby enables all chemical properties at each residue, as shown in Fig. 4 b, d. By uniformly sampling all DAASs, RSM enables much more chemical diversity in the library than traditional random mutagenesis methods, enhancing the likelihood to discover dramatically improved variants. Functional Screening of RSM library To illustrate the potential of RSM, we evolve BsLipA to substantially improve its thermotolerance. Lipases such as BsLipA are versatile enzymes that catalyze the hydrolysis and synthesis of esters, playing crucial roles in various industrial applications, including biodiesel production, and use in detergents. However, their practical utility is often limited by the high temperatures used in industrial processes (Vivek, 2022). The wild-type BsLipA exhibits a relatively low melting temperature (T m ), and irreversibly loses activity upon heating. Thus, the development of lipases resistant to heat inactivation will significantly expand their industrial applications, enabling their use across a broader range of high-temperature processes. To make significant improvement in the thermotolerance of BsLipA, we maximize the diversity of our library in two ways: we use RSM to ensure high chemical diversity at each residue position, and choose a high mutation rate so each pair of variants, on average, are different by more residues. However, as the mutation rate of a library increases, the theoretical sequence space grows rapidly, and the best variants can become too rare to be discovered by the experimentally practical screening throughput. With the screening throughput of droplet sorting, which is ~ 10 7 , there exists a maximum mutation rate with which the best mutants can still be discovered. To estimate this maximum mutation rate for our problem, we assume that the most thermotolerant mutants in a library can result from the combination of single mutations that individually improve thermotolerance. The frequency of such variants in a high mutation rate library can then be estimated by the proportion of single mutants that exhibit higher thermal stability than the wild-type. This proportion can be inferred by constructing a single saturation library using RSM and performing a droplet cytometry analysis on this library. For BsLipA, a deep mutational scan of all single mutants has been done to measure the change in thermotolerance over the wild-type 9 . Therefore, we count the single mutants in this dataset that show improved thermotolerance over the wild-type, which comprises ~ 0.1 of all single mutants. Using this proportion, we estimate the frequency of mutants in an n-mutation RSM library that combine all individually beneficial mutations to be (0.1) n . For more than one of such mutants to be discovered, we require this frequency to be higher than one over the screening throughput of 10 7 . Therefore, we choose a mutation rate of ~ 7. To identify BsLipA variants that retain activity after thermal exposure, we express the RSM library in E. coli host and use a fluorescence-activated droplet sorting (FADS) assay to screen for thermotolerant candidates. Each E. coli cell is co-encapsulated with lysis reagents and the fluorogenic substrate, resorufin acetate, in a 20 µm droplet. The droplets are incubated at 65°C for one hour, during which time the cells are lysed and the BsLipA proteins are released. Variants that remain active at 65°C catalyze the conversion of resorufin acetate into the fluorescent product, resorufin. After incubation, the droplets are cooled to 25°C and reinjected to a sorter for high-throughput screening, as shown in Fig. 5 a. The sorting experiment, conducted at 25°C, can also capture variants that might regain activity during the cooling phase, further converting resorufin acetate. We screen ~ 10 7 droplets and sort the 0.04% of droplets with the highest fluorescence intensity. The sorted droplets are merged, and the DNA is amplified and used as input for a second round of cell sorting; the repeated sort reduces the effects of variations in cell expression. In the second round, we screen ~ 3x10 6 droplets and sort the 0.04% of droplets. To isolate individual variants, we transform the sorted plasmids into E. coli and plate the cells on agar. Single colonies are then picked and grown in separate wells of a 96-well plate. To coarsely assess the thermotolerance of the sorted variants, we lyse the cells using lysis reagents, incubate the lysates at 80°C for one hour, and then perform a fluorescence-based assay on a plate reader to measure residual activity at 25°C from the cell lysates (Fig. S6.). We identify two variants with highest residual activity, named mutant 5 (MT5) and mutant 20 (MT20). Both MT5 and MT20 contain seven amino acid substitutions, as shown in Fig. 5 b. Each single amino acid substitution, except for the S163P substitution in MT20, individually improves thermotolerance 9 , consistent with the assumption that very good mutants can result from the combination of beneficial single substitutions. These two variants as well as the wild-type (WT) are subsequently purified using NiNTA mini-spin columns for further characterization. Characterization of Thermotolerant BsLipA Variants To evaluate the improvement of heat tolerance of MT5 and MT20, we incubate the purified enzymes at various temperatures for 20 minutes and measure the residual activity after cooling them back to room temperature. To enable comparison, we normalize all residual activity values to the WT activity after incubation at 25°C. The WT exhibits constant residual activity after incubation at temperatures from 25°C to 45°C. However, at higher temperatures, the residual activity declines rapidly and remains low from 55°C to 95°C, as shown in Fig. 5 c. To characterize this sharp drop, we interpolate the data to estimate T 50 , the temperature at which 50% activity is lost as compared to the lowest incubation temperature of 25°C. The value of T 50 is commonly used to compare the thermotolerance of enzymes 16 , 17 , 20 . For the WT, T 50 ≈ 51.5°C, as shown in Fig. 5 c. In contrast, MT5 maintains high residual activity over a wide temperature range, showing only a gradual decrease at the highest temperatures as shown in Fig. 5 c. Similarly, MT20 maintains high residual activity up to 75°C, with a decline of approximately 35% at 95°C as shown in Fig. 5 c. Thus, the T 50 values of both variants are beyond the range of our measurement. In addition, MT20 also demonstrates nearly two-fold higher enzymatic activity than WT at room temperature. To understand the origin of the different thermotolerance of the WT and the variants, we measure the structural stability at elevated temperatures using Nano Differential Scanning Fluorimetry (NanoDSF). During a temperature ramp from 20°C to 95°C, we monitor the red shift of intrinsic fluorescence by measuring the change in R f , the fluorescence intensity ratio between 350 nm and 330nm. An increase in R f reflects the increase of polarity in the chemical microenvironment of tyrosines and tryptophans, thereby indicating unfolding of the protein. For the WT, we observe a single inflection in the temperature dependence of R f indicating an unfolding transition at T m = 54.7°C, as shown in Fig. 5 d. The value of T m of the WT roughly coincides with the value of T 50 , suggesting that the enzyme activity is irreversibly lost following thermal unfolding. This behavior is attributed to irreversible aggregation upon thermal unfolding 17 ; this behavior is irreversible, as shown by the NanoDSF, where R f exhibits no defined transition upon cooling from 95°C to 20°C, indicating the absence of refolding, as shown in Fig. 5 e. As the temperature approaches the T m of the WT, its native tertiary structure progressively unfolds. Unlike the native state (N) where the hydrophobic residues are buried within the protein structure, certain intermediate denatured states (I) expose the hydrophobic residues on the surface and lead to irreversible aggregation (A), where the aggregates remain stable even when the temperature decreases to room temperature. This irreversible aggregation inhibits the refolding of the native structure and prevents the recovery of enzymatic activity, as shown schematically in Fig. 5 f. Therefore, the heat inactivation of the WT involves two key steps: thermal unfolding and irreversible aggregation of unfolded states. Since MT5 and MT20 both show exceptional heat tolerance as compared to the WT, they must counteract at least one of the two steps in the pathway of the WT heat inactivation. We first examine using NanoDSF if the thermal unfolding of the variants is shifted to much higher temperatures. Both variants exhibit one thermal transition, with T m values of 61.2°C for MT5 and 57.4°C for MT20, as shown in Fig. 5 d. The increase in T m demonstrates the higher structural stability of the variants as compared to the WT. However, this moderate improvement in structural stability is not sufficient to explain the > 40°C difference in T 50 between the variants and the WT. Unlike the WT, which irreversibly aggregates and loses activity after incubation around or above its T m , both variants can experience temperatures > 30°C above their T m values, and still recover more than half of their activity upon temperature reduction. Therefore, the variants must disturb the second step of the WT heat inactivation process. At temperatures above T m , where the native structure is lost, the variants likely adopt mechanisms to avoid aggregation-prone conformations and remain monomeric. The structural reversibility of both variants is verified by the defined transitions in their NanoDSF refolding ramps, as shown in Fig. 5 e. By circumventing the stable aggregated state, they retain the ability to refold into functional structures upon temperature reduction, as shown schematically in Fig. 5 g. To confirm the aggregation resistance of the variants, we incubate both variants at 65°C for 10 minutes, centrifuge each sample to remove the aggregated fraction and measure the protein concentration that is maintained in each solution. Each variant retains > 90% protein mass in the supernatant. In contrast, the WT retains only ~ 50% protein mass in the supernatant, as shown in Fig. 5 h. These characterizations confirm two mechanisms for heat tolerance of MT5 and MT20: improved thermal stability of the native structure, and reduced propensity for aggregation in thermally unfolded conformations. The first mechanism accounts for the T m improvements of 6.5°C and 2.7°C, whereas the second mechanism is the major contribution to the > 40°C improvements in T 50 . Molecular Dynamics Simulations To understand the role of these substitutions on the thermotolerance, we perform molecular dynamics (MD) simulations for the WT and for MT20. The catalytic triad of BsLipA is formed by S77, D133, and H156 located in flexible loops enclosing the active pocket, as shown in Fig. 6 a. Both thermally stable variants contain mutations near the catalytic residues 133 and 156, as well as mutations at position 17, located on the loop considered to mimic a typical lipase lid absent in BsLipA 21 , as shown in Fig. 6 a. For the WT, structure 5CRI from Protein Data Bank (PDB) is used as the starting structure and equilibrated at 25°C in the MD simulations. To obtain a starting structure for the MT20, we introduce the seven amino acid substitutions into the equilibrated WT structure using SCWRL4 and re-equilibrate it. We perform 200 ns of MD simulations at 25°C for both the WT and the MT20 and use trajectory information to create representative 3D structures. The representative structure of MT20 closely resembles that of the WT as indicated by a high degree of structural overlap, as shown in Fig. 6 b. This similarity indicates that the mutations do not significantly alter the overall protein fold. However, by aligning the WT and MT structures and visualizing interactions present only in MT20, we identify new hydrogen bonds stabilizing the MT20 structure, as indicated by light blue lines in Fig. 6 b. These include a hydrogen bond between V154 and S130 each next to a catalytic residue, stabilizing the catalytically active conformation, and another hydrogen bond between K170 and the mutated S174 residue, stabilizing the terminal helix. These interactions may contribute to the thermal stabilization of the MT20, which exhibits a melting temperature 2.7°C higher than the WT. To investigate the conformational change of the WT and MT20 at elevated temperatures, we increase the simulation temperature to 100°C. This impractically high temperature is chosen to capture trends of significant conformational changes within the simulation time 22 – 25 . We use the first 20 ns and the last 20 ns from a 200-ns production run at 100°C to create representative 3D structures. For the WT, we observe a significant conformational change at 100°C as compared to the 25°C structure. At room temp, the WT has its hydrophobic residues enclosing the active site facing inward, as shown in Fig. 6 c. When the temperature increases, the flexibility of two loops (131–137 and 152–156) that enclose the active pocket increases. As a result, highly hydrophobic residues M134, I135, V136, M137 and Y139 are exposed to the surface of the protein, increasing the aggregation propensity, as shown in Fig. 6 d. Moreover, the loop between β7 and αE, which contains D133, folds into the extension of αE that now includes residues 134–141, as shown in Fig. 6 e. The extended αE turns D133 away from H156 and results in a catalytically inactive conformation. This newly formed helix persists in the structure until the temperature is increased to 180°C and the time of simulation is prolonged to 1 µs, suggesting it may represent a metastable intermediate state in the thermal denaturation pathway of the enzyme. This conformation may be a representative intermediate state that leads to the irreversible aggregation of the WT upon heating. However, for MT20, an increase in temperature during the simulation only increases the movement of flexible regions. No significant structural changes around the active pocket are identified, as shown in Fig. 6 f, and no additional hydrophobic residues are exposed as compared to the native structure, as shown in Fig. 6 g and 6 h. The lack of αE extension at 134–137 suggests that the mutations near this loop, including M134P, M137A and Y139S, disturb its formation. As proline has the lowest α-helix propensity of all amino acids 26 , the M134P substitution potentially perturbs the α-helix extension of residues 134–137. Further stabilization of this loop is achieved by substituting large nonpolar residues M137 and Y139 with small or polar amino acids A and S. By rejecting the intermediate state in the WT where a large hydrophobic patch is exposed, MT20 avoids aggregation during thermal unfolding. We therefore conclude that the heat tolerance of BsLipA in MT20 is achieved through two molecular mechanisms: the formation of new hydrogen bonds and the absence of aggregation-prone intermediates during thermal unfolding. The latter can be considered a predominant stabilization factor for activity retention after exposure to high temperatures 27 . While we have not run an MD simulation for MT5, we postulate that MT5 adopts aggregation resistance from its four substitutions around the loop between β7 and αE (M134E, I135E, M137D, and Y139D). These substitutions replace nonpolar residues with negatively charged ones, likely generating electronic repulsion to prevent the formation of αE extension and the aggregation of denatured conformers 16 , 17 . In addition to MT5 and MT20, six other variants of BsLipA have been reported to exhibit aggregation resistance upon thermal unfolding 16 , 17 , 27 . Each of these eight variants carries at least one methionine substitution at positions 134 or 137. Variants 5A and 5B 17 , each carries 10 substitutions, but only the M134E or M137P substitution results in the resistance to thermal aggregation 17 . Such experimental evidence, along with our MD simulation demonstrates the critical role of residues 134 and 137 in conferring aggregation resistance upon thermal unfolding. Notably, among the substitutions at residues 134 and 137 that occur in these eight variants, M134D, M134E, M134P, M137P and M137A, none could be achieved by epPCR with single NT changes, as shown in Fig. 4 a. In fact, when starting from a wild-type methionine, the categories of negatively charged residues, aromatic residues, and special cases G, P, C are all barely accessible by epPCR, as shown in Fig. 4 c. The discovery of each one of these eight variants relied on saturation mutagenesis (Extended Data Table 1). In particular, the parent of 5A and 5B, 4D3, was discovered after four rounds of epPCR, four rounds of screening for thermotolerance, and three rounds of recombination. However, 4D3 only achieved a T 50 improvement of 15°C, and did not acquire resistance to thermal aggregation when heated to around its melting temperature. Only through site-saturation mutagenesis at residues 134 and 137 of 4D3 were the aggregation-resistant variants 5A and 5B identified. The mutagenesis trajectory leading to these variants underscores the importance of saturating the chemical diversity at key residues in mutagenesis libraries. By using RSM to prepare the library, identification of key residue positions and the saturation of their chemical diversity can be achieved simultaneously. Discussion Mutagenesis libraries are the starting point of directed evolution campaigns and thereby determine the range of the improvement. With no prior information about the wild-type, a random mutagenesis method that uniformly samples all possible amino acid substitutions across all residue positions is critical to preserve the accessibility to beneficial properties. For an enzyme with N residues, the number of all possible single amino acid substitutions, 19xN, is typically a few thousand. Therefore, the number of unique variants with n substitutions is ~(19xN) n . With only three random substitutions per gene, the theoretical sequence space exceeds the screening throughput of most current screening techniques 28 . For a given screen to effectively cover a meaningful subspace of an n-mutation sequence space, it is therefore important that the mutagenesis method allows mutation rate control to adjust the expected hit rate to match the screening throughput. All three critical features of a random mutagenesis method: mutating across the gene, saturating the diversity of each residue, and controlling mutation rates, are simultaneously provided by RSM. When the screening throughput is limited to a few thousand, a single or low mutation rate RSM library can be prepared to enable mutational analysis of single amino acid substitutions and to guide recombination. If screening with even higher throughput is available, a high mutation rate RSM library can be implemented, directly enabling the discovery of significantly improved variants combining multiple beneficial mutations. The seven-mutation library used in our study has a theoretical diversity of (181 choose 7) x 20 7 = 1.4 x 10 21 . With a throughput of 10 7 , our screen samples only one in 10 15 in this sequence space. The discovery of significantly improved variants with such low sample rate suggests that similar high-performing variants are likely prevalent in the sequence space; therefore, it is possible to capture them with substantial down sampling. However, the sampling scheme, determined by the library preparation method, must retain the accessibility of each residue to all alternative chemical properties. If the library preparation method inherently causes significant loss of residue properties, such as with epPCR, the crucial substitutions may be eliminated from the accessible sequence space once the library is constructed. By enforcing uniform sampling, RSM reduces this bias and preserves a better representation of the total chemical space. The combination of a high mutation rate RSM library with high-throughput droplet screening enables us to more uniformly sample the vast sequence space and identify two variants with T 50 values improved by > 40°C compared to the wild-type enzyme in a single mutagenesis round. By contrast, using other methods, similar variants can only be achieved after 4–5 rounds of mutagenesis, recombination and screening (Extended Data Table 1). When some prior information is available, RSM can be adjusted to reflect this knowledge and better facilitate the search of improved variants. For example, if previous studies have identified certain regions in the wild-type enzyme in conferring certain properties, oligo pools that cover those regions can be used in one mutagenesis reaction while the remaining oligo pools can be used in other mutagenesis reactions. When mixing mutagenesis products for StEP, the product of the first reaction can be used with a much higher molar ratio than the product of the second reaction, resulting in higher mutation rates in the previously-identified regions. If prior information further identifies beneficial properties at certain positions, the set of primers to mutate these positions can be designed to target substituted amino acids with those properties. These adjustments of RSM can focus the mutations to specific regions and specific substitutions, further improving the likelihood of discovering significantly improved variants. These results demonstrate the potential of RSM to accelerate directed evolution. Material and Methods Microfluidics device fabrication. We design the devices (droplet generator and droplet sorter; see Fig. S12) using AutoCAD and producing them as photomasks (CAD/Art Services, Inc.). We fabricate the devices through the well-established techniques of soft lithography, employing SU8-on-silicon-wafer masters and PDMS-on-glass devices. We pour polydimethylsiloxane (PDMS) (Sylgard 184) onto the masters and bake them at 65°C overnight to cure the PDMS. Subsequently, we cautiously peel each PDMS device from the master and seal it to a pristine glass slide (Corning, 2947). For devices incorporating electrodes (droplet sorter), we integrate the electrodes into the design as channels within the microfluidic devices. We fill these channels with a low melting point metal alloy (Indalloy 19, 0.020 in. diameter) while heating the devices on a hot plate. We add terminal blocks to the punched holes in the devices to facilitate electrical connections during experiments. Following the fabrication process, we inject Aquapel (fluoroalkylsilanes) through the punched holes of the devices and use pressurized air to impel it through the microfluidic channel walls, rendering the internal surfaces of the channels hydrophobic. We expel excess Aquapel using compressed air and bake the device at 65°C overnight. Library Construction and Sequencing For the single library, we prepare 18 parallel multisite-directed mutagenesis reactions in PCR strips, each containing one primer pool as the primers. After thermocycling, we purify the DNA product and perform overlapping PCR with wild-type BsLipA on the pET28a + vector backbone as the template. We then incubate the PCR product with DpnI to cleave methylated wild-type plasmid templates and column purify it using the Monarch DNA purification kit (New England Biolabs). We transform the plasmid library into competent cells (Turbo chemically competent cells, New England Biolabs), plate a portion to verify successful transformation, and transfer the remainder into liquid culture for plasmid extraction using the Monarch Plasmid extraction kit (New England Biolabs). For the 4×3 + 3×2 library, we prepare five parallel multisite-directed mutagenesis reactions, each containing three or four primer groups. Following thermocycling, we perform a StEP reaction on the purified DNA products to induce recombination of the resultant mutants. We prepare the reaction mix for StEP, consisting of a DNA shuffling mixture with 2.5 µL of template (1 ng/µL), 50 µL of 2× Taq Master Mix (New England Biolabs), 1.5 µL of 10 µM forward primer (TAATACGACTCACTATAGGG), 1.5 µL of 10 µM reverse primer (GCTAGTTATTGCTCAGCGG), and MilliQ water to a total volume of 100 µL. We use a thermocycling program of 94°C for 30 seconds and 55°C for 5 seconds for 99 cycles. After thermocycling, we use a portion of the DNA product for gel electrophoresis to confirm its size and purify the remaining products using a DNA cleanup kit. The purified PCR products serve as megaprimers in an overlap PCR reaction with 250 ng of PCR products, 300 ng of wild-type BsLipA plasmids, 25 µL of 2× Q5 Master Mix, and MilliQ water to a total volume of 50 µL. We use a thermocycling program of 98°C for 2 minutes, followed by 98°C for 10 seconds, 61°C for 20 seconds, 72°C for 3 minutes, and a final extension at 72°C for 5 minutes. After the StEP and overlapping PCR, we add 1 µL of DpnI and incubate at 37°C for 3 hours. We purify the PCR product using a PCR cleanup kit, elute it into 10 µL, and transform it into Agilent electrocompetent cells. After adding 900 µL of SOC and incubating at 37°C and 250 rpm for 1 hour, we plate 20 µL on LB-trimethoprim plates and add 3 mL of LB-trimethoprim to a liquid culture for overnight growth at 37°C. The isolated plasmids are ready for further analysis. Similarly, for the 6×3 library, we pool six primer groups at a time, creating a maximum of six mutations per round, followed by StEP. We purify the DNA product and perform overlap PCR with the wild-type plasmid as the template, which we later remove by DpnI digestion after the reaction. We then column purify the resultant DNA product, transform it into competent cells, plate a portion to verify transformation, and transfer the rest into liquid culture to grow overnight for plasmid extraction. For the 9×2 library, we add odd-numbered primer pools first, then even-numbered pools, followed by StEP, overlap PCR, DpnI digestion, DNA column purification, transformation, and overnight growth in liquid culture for plasmid extraction. For the 18×1 library, we add all primer pools together at once, creating a maximum of 18 mutations per round, followed by StEP, overlap PCR, DpnI digestion, DNA column purification, transformation, and overnight growth in liquid culture for plasmid extraction. For thermotolerance variant screening, we prepare a library with a higher mutation rate by using two consecutive 9×2 RSM protocols. We also construct three error-prone libraries using the GeneMorph II kit (Agilent) with high, medium, and low mutation rates, resulting in 2.70, 2.04, and 1.18 mutations per gene, respectively, following the manufacturer’s manual. We transform the plasmid library into competent cells, plate a portion to verify successful transformation, transfer the remainder into liquid culture, and extract plasmids. We send all the libraries (single, 4×3 + 3×2, 6×3, 9×2, 18×1, EPhigh, EPmid, EPlow) for sequencing using the PacBio platform for high-fidelity sequencing of long reads. Sequencing Result Analysis We filter sequencing reads based on the quality score of the lowest quality base in the open reading frame with a threshold of 90, eliminating any imperfect sequencing result. After filtering, we acquire approximately 10 5 reads for each library. We align the reads to the wild-type open reading frame, translate the aligned sequences into polypeptide sequence and compare the translation to wild-type translation. All analysis is conducted using a customized python package based on open source libraries such as Biopython (github.com/angiechen918/Mutan.git) Statistical Analysis on Library quality To analyze the percentage of reads with insertions and deletions, we filter all reads with aligned length of wild-type length +/- 3 and count the fraction of reads with either insertions or deletions. To analyze the percentage of reads with in-frame stop codons, we filter all reads with the correct length and count the fraction of reads that translates into a polypeptide with shorter length than the wild-type enzyme. To analyze the percentage of silenced mutations, we filter all reads with the correct length and identify the codon positions that contain nucleotide changes. For all codon positions that contain nucleotide changes, we count the fraction of codon positions that does not lead to a different amino acid than the wild-type enzyme at this position. Statistical Analysis on Distinct Amino Acid Substitution Abundance in Different Libraries We filter all reads with the correct length and compare the translation to the wild-type enzyme to identify single amino acid substitutions. For each library, we analyze ~ 72k single amino acid substitutions and categorize them based on the amino acid that they “mutate from” in the wild-type (19 amino acids, 20 - cysteine) and the amino acid that they “mutate to” in the variant (19 amino acids, 20 - the “mutate from” amino acid). We each of these 19x19 = 361 categories a chemically distinct substitution. As each amino acid may be used by the wild-type enzyme multiple times, we compile a table where the rows list the 19 amino acids in the wild-type BsLipA and the columns list the 19 possible mutations, resulting in a heatmap of the 361 chemically distinct amino acid substitutions. The color gradient reflects the abundance of the mutations, with warmer colors representing higher mutation frequencies. The numbers in each square of the heatmap represent the normalized count of each amino acid substitution. A number of 0 indicates that the substitution is enabled but occurs so infrequently that it rounds to 0 after normalization. Amino acid substitutions that can be enabled by a single nucleotide change in a trinucleotide codon are highlighted by white-bordered squares. There are 124 such squares, indicating that only 34.4% of possible substitutions are enabled by single nucleotide changes. Statistical Analysis on Impact of Staggered Extension Process (StEP) on Mutation Distribution and Correlation To decrease the positional correlation and increase the diversity of the final RSM library, we mix the products from the two reactions in equimolar ratios as templates for StEP. Before StEP, the number of mutations per region is limited to one, and the mutations are fully confined within the regions covered by the primer pools of choice (Sup.Fig. S1 ). In the StEP thermocycling program, we limit the extension step to a few seconds, preventing the DNA polymerase from synthesizing the full gene in each cycle. During the denaturing step of the next cycle, incomplete strands unbind from the original template and anneal to different template strands, effectively combining mutations from various templates. Therefore, after StEP, the positions for mutations are no longer fully correlated. We reflect this change by counting the number of reads where a specific pair of mutation positions appear on the same read, N(A,B), choosing A > B. For example, if a variant in the library contains three mutations at positions 1, 21, and 31, then this read contributes counts to N(1,21), N(1,31), and N(21,31). When mutations at position A coincide frequently with mutations at position B, we expect a high N(A,B) value. By contrast, if the mutations at positions A and B never coincide on the same read, N(A,B) should be 0. The abundance of individual positions A and B also strongly influences N(A,B). To calculate the association between positions without considering the abundance of individual positions, we divide the number of position pairs by the occurrence of either position in the pair and define this value as Mutual Association, M(A,B) = N(A)×N(B)N(A,B)​. This value is similar to the exponential of pointwise mutual information (PMI), a concept broadly used in Natural Language Processing (NLP) to evaluate the association of two words in a corpus (Sup.Fig. S2). In the 9×2 library before StEP, M(A,B) values are high if positions A and B are in pools with indices of the same parity, resulting from the grouping of pools in the mutagenesis reaction. By contrast, if A and B are in pools with indices of opposite parity, M(A,B) values are close to 0, as shown in Fig. 2 e. However, this pattern in the association of positions becomes relaxed after StEP. In the 9×2 library after StEP, as long as the distance between positions A and B is longer than ~ 30 nucleotides, M(A,B) values are uniform, indicating low bias in position-wise association, as shown in Fig. 2 f. For position pairs with short distances, positions in pools with the same indices of opposite parity are still more associated than positions in pools with opposite indices of opposite parity. We expect such bias to be further eliminated with shorter annealing times during StEP or more rounds of StEP. Notably, despite the dramatic shift in mutation distribution with StEP, the average mutation rate of the 9×2 library remains consistent before (4.87 ± 2.57) and after StEP (4.64 ± 2.40), ensuring predictability when creating a mutagenesis library with a target mutation rate. Droplet Sorting We co-encapsulate each bacterial host cell expressing the mutant protein with resorufin acetate substrate and a bacterial lysis reagent (BugBuster Protein Extraction Reagent, Sigma-Aldrich) using a dual-inlet droplet maker (Fig. S7). We collect the emulsion under a layer of mineral oil into an Eppendorf tube sitting in a heat block set to 65°C and incubate it for 1 hour to inactivate the non-thermostable variants. Next, we carefully remove all the mineral oil with a pipette and re-inject the emulsions into a microfluidic sorting device (Fig. S7). With an additional stream of spacing oil, the closely packed droplets are evenly spaced and enter the sorter periodically. We place a 532 nm excitation laser at the entrance of the sorter. If the fluorescent reaction product in a droplet exceeds the set fluorescence threshold, we use an electric field to pull it into an adjacent channel and collect it through polyethylene micro-tubing into an Eppendorf tube placed on ice. We operate the chip at 600–1000 drops per second, probing approximately 3 million cells per hour. By applying electric fields of 50 V at a frequency of 25 kHz, we screen more than 3 million droplets and sort 1,260 droplets. We break the emulsion by adding 200 µL of 20% (v/v) PFO (1H,1H,2H,2H-perfluoro-1-octanol, Sigma-Aldrich) in HFE 7500 (3M) and purify the aqueous phase using a DNA cleanup kit (Monarch® PCR & DNA Cleanup Kit, New England Biolabs), eluting the DNA in 5 µL of ddH₂O. We then transform the plasmids into high-efficiency competent cells (Turbo Competent E. coli , New England Biolabs) for plasmid replication and extraction. Residual Activity Measurement in Cell Lysate We transform the plasmids from the sorted library into BL21(DE3) competent cells (NEB C2527H) and plate 200 µL of the culture onto agarose LB-kanamycin plates, incubating them at 37°C overnight for colony growth. We pick 94 colonies and transfer them to a 96-well plate, including one TBS plus substrate blank and a wild-type control. We shake the plate at 37°C for 12 hours to allow the colonies to grow overnight in 200 µL of LB-kanamycin liquid culture. We then dilute the culture 20-fold, with triple replicates, and add BugBuster reagent to each well. We incubate the plate at 37°C for one hour, then increase the temperature to 70°C and incubate for an additional hour to inactivate non-thermostable variants. Subsequently, we add 7.5 µM resorufin acetate while keeping the plate on ice throughout the procedure. We measure the relative fluorescence units (RFU) at 549 nm excitation and 585 nm emission over time. Two variants are identified with residual activity significantly higher than the wild type after incubation at 70°C for one hour. Protein expression and purification We grow the identified variants in liquid culture and induce them using Invitrogen's MagicMedia E. coli Expression Medium at 30°C overnight. We purify the protein using Ni-NTA spin columns (Thermo Scientific, 88224). To remove imidazole, we dialyze the eluted sample against 1X DPBS at 4°C using for 4 hours a 7k MWCO dialysis cassette (Thermo Scientific 66373). After dialysis, we use SDS-PAGE to confirm the molecular weight of the protein at the correct size. To preserve the purified enzymes, we standardize their concentrations to 0.4 mg/mL in 50% glycerol in 1x DPBS and store them at -20°C before use. We then characterize these purified proteins for their T₅₀, residual activity, and melting temperature. T 50 Measurement To measure the T 50 of enzyme mutants, we dilute the enzyme stocks 100-fold by mixing 10 µL of the enzyme with 1,490 µL of 1× TBS buffer. We then incubate the diluted enzyme solution on three PCR machines simultaneously at various temperatures (gradient 35–60°C, gradient 65–90°C, and 95°C) for 20 minutes, with a total volume of 120 µL in each PCR tube, selecting seven temperature points: 36, 45, 55, 66, 75, 85, 95°C, and leaving one sample at room temperature (25°C). Following incubation, we further dilute the enzyme mutants 10-fold by adding 12 µL of the enzyme solution to 108 µL of TBS buffer. We then add this diluted enzyme to each well of a 96-well plate. We prepare the substrate by diluting 20 µL of 10 mM resorufin acetate (RA) in 10 mL of 1× TBS buffer, freshly prepared for each group. In each well of the 96-well plate, we add 12 µL of the 1,000-fold diluted enzyme and 12 µL of the substrate, ensuring thorough mixing. We initiate the measurement program at 25°C on a qPCR machine (Bio-Rad CFX96 Real-Time PCR System), with each point measured for 15 seconds using the Texas Red channel. The groups tested include mutants 5 and 20, wild type, and a glycerol-only control. Residual Activity at High Temperature Measurement with Purified Protein To assess both the residual activity and the high-temperature enzyme activity of the identified Bacillus subtilis Lipase A (BsLipA) variants, we use the Bio-Rad CFX96 real-time PCR system to measure fluorescent products at various temperatures. The groups tested include mutants 5 and 20, wild type, and a glycerol-only control. First, we dilute purified enzyme stocks (0.4 mg/mL) in 1× TBS buffer (pH 7.4) at a ratio of 1:500 (2 µL in 1,000 µL), adding an additional 10 µL of 50% glycerol to match the glycerol concentration. For the WT enzyme, which has an initial concentration of 0.5 mg/mL, we first dilute it at a ratio of 4:5, then follow the same dilution protocol as for the mutants. Next, we subject the enzyme solutions to heat inactivation at 80°C for 1 hour using Eppendorf tubes on a heat block. Immediately after heating, we place the tubes on ice to stop the heat inactivation process. We prepare the substrate by diluting 20 µL of 25 mM fluorescein diacetate (FDA) in 10 mL of 1× TBS buffer, freshly prepared for each group. In each well of a 96-well plate, we add 12 µL of the 1,000-fold diluted enzyme and 12 µL of the substrate, ensuring thorough mixing. We initiate the measurement program at 41, 45.6, 50.2, 56, 60.8, and 65°C on a qPCR machine (Bio-Rad CFX96 Real-Time PCR System), with each point measured for 15 seconds using the SYBR channel. Melting Temperature ( T m ) Measurement We use NanoDSF to measure the melting points of the purified proteins MT5, MT20, and WT. We load prepared samples into a 384-well plate and then into NanoTemper NT.Plex Capillary Chips. We perform the experiments in duplicates. We heat and cool the samples from 20°C to 95°C to 20°C at a rate of 1.5°C/min and collect the data using the NanoTemper PR.ThermControl software. NanoDSF is an advanced differential scanning fluorimetry technique that monitors changes in intrinsic fluorescence of tryptophan and tyrosine residues as proteins unfold. As proteins are exposed to a temperature gradient, their local chemical environment changes, causing variations in fluorescence intensity and wavelength, particularly when aromatic side chains become solvent-exposed upon unfolding. NanoDSF measures fluorescence at 330 nm and 350 nm and calculates the ratio as a function of temperature to determine the melting temperature (T m ), which is the temperature at which 50% of the protein population is unfolded. NanoDSF offers higher throughput and lower sample consumption compared to traditional methods like differential scanning calorimetry (DSC) or circular dichroism (CD), with comparable accuracy for determining folding free energies and unfolding temperatures. NanoDSF requires the presence of aromatic amino acids in the protein of interest. Measurement of Soluble Protein Fraction after Heat Inactivation To assess protein solubility upon heat incubation, we dilute purified enzyme stocks MT5 and MT20, along with wild-type (WT) protein (0.4 mg/ml in 50% glycerol 1xDPBS), to 0.025 mg/ml in 1x DPBS buffer. We subject the diluted enzyme solutions to heat inactivation at 65°C for 10 minutes in Eppendorf tubes placed on a heat block. After incubation, we centrifuge the heated samples along with non-heated controls at 12,000 × g for 10 minutes to separate soluble proteins from precipitated aggregates. We carefully collect the supernatants to determine the concentration of soluble protein using Micro BCA (Thermo Scientific, 23235) following the manufacturer instruction. The absorption at 562 nm is measured using a NanoDrop spectrophotometer. Molecular Dynamics Simulations of Protein Folding at Low and High Temperatures For molecular dynamics simulations we use crystal structure of BsLipA PDB ID: 5CRI. For molecular dynamics simulations, protein preparation and analysis we use Amber22 29 software with GPU support and AmberTools23 30 suite of programs. We use ff19SB protein force field 31 and OPC 4-point water model 32 . We initially minimize the system for 2000 steps using Newton-Raphson algorithm, and then slowly heat it from 100 K to final simulation temperature for 1 ns with time-step of 1 fs, while keeping the volume constant. The initial restraint on protein atoms is 100 kcal/mol-Å 2 . We then relax the system for total of 6 ns with 1 fs time step, in NPT ensemble, by slowly lifting restraints (to 10.0, 1.0 and 0.1 kcal/mol-Å 2 ) on protein backbone atoms. Finally, we perform 20 ns of equilibration in NPT ensemble with 1 fs time step and no restraints. We perform 200 ns production runs using NPT ensemble with Langevin thermostat and Berendsen barostat using 2 fs time step. To impose constraints on hydrogen atoms we employ SHAKE algorithm. To visualize MD trajectories, we use VMD 1.9.3 software. We introduce mutations into wild-type protein using SCWRL4.0 33 software. We perform molecular dynamics simulations at 25°C and increased temperatures for both M20 and WT to 100°C. We perform all simulations in triplicates. We use 200 ns trajectory information to create representative 3D structures for MT20 and WT using SINAPs extension in UCSF Chimera 1.17.3 34 . We align these structures and visualize interactions present only in the mutated protein. Additionally, we use the first 20 ns after equilibration and the last 20 ns from simulations at 100°C to create representative 3D structures for MT20 and WT and overlap them to visualize structural differences. We also use UCSF Chimera to visualize hydrophobic surfaces in representative 3D structures. Declarations Author Contributions X.D.Z., A.C., K.M. and D.A.W. developed the concept. X.D.Z, A.C., R.H., P.E. and K.M. designed the experiments and analysis approaches. X.D.Z, A.C., LC, Y.X., R.H., P.E. and K.M performed the experiments and conducted the analyses. MM performed and analyzed computational simulations. X.D.Z, A.C., D.A.W. and K.M. wrote the manuscript with input from all authors. Acknowledgments This work was supported in part by the National Science Foundation (CBET-2103538), the Harvard MRSEC (DMR-2011754), the Wyss Institute and Science Fund of the Republic of Serbia, #7750288. Conflict of Interest David A. Weitz and Karla Milcic are co-founders of RevivBio Inc. Data Availability Statement The data that support the findings of this study are available in the supplementary material and extended data of this article and upon request from corresponding authors. Code Availability Statement All analysis of PacBio sequencing data is conducted using a customized python package based on open-source libraries such as Biopython and is available at: github.com/angiechen918/Mutan.git References Packer, M. S. & Liu, D. R. Methods for the directed evolution of proteins. Nature Reviews Genetics 16 , 379-394 (2015). https://doi.org:10.1038/nrg3927 Esvelt, K. M., Carlson, J. C. & Liu, D. R. A system for the continuous directed evolution of biomolecules. 472 , 503 (2011). Rosenthal, R. G., Zhang, X. D., Đurđić, K. I., Collins, J. J. & Weitz, D. A. Controlled Continuous Evolution of Enzymatic Activity Screened at Ultrahigh Throughput Using Drop-Based Microfluidics. 62 (2023). Drummond, D. A., Iverson, B. L., Georgiou, G. & Arnold, F. H. Why high-error-rate random mutagenesis libraries are enriched in functional and improved proteins. 350 , 816 (2005). Wong, T. S., Roccatano, D., Zacharias, M. & Schwaneberg, U. A Statistical Analysis of Random Mutagenesis Methods Used for Directed Protein Evolution. Journal of Molecular Biology 355 , 858-871 (2006). https://doi.org:10.1016/j.jmb.2005.10.082 Frauenkron-Machedjou, V. J. et al. Towards Understanding Directed Evolution: More than Half of All Amino Acid Positions Contribute to Ionic Liquid Resistance of Bacillus subtilis Lipase A. 16 , 945 (2015). Frauenkron-Machedjou, V. J. et al. Exploring the full natural diversity of single amino acid exchange reveals that 40–60% of BSLA positions improve organic solvents resistance. 5 (2018). Markel, U. et al. Are directed evolution approaches efficient in exploring nature’s potential to stabilize a lipase in organic cosolvents? 7 (2017). Nutschel, C., Fulton, A., Zimmermann, O., Schwaneberg, U., Jaeger, K. E. & Gohlke, H. Systematically Scrutinizing the Impact of Substitution Sites on Thermostability and Detergent Tolerance for Bacillus subtilis Lipase A. 60 , 1568-1584 (2020). Cheng, F., Zhu, L. & Schwaneberg, U. Directed evolution 2.0: Improving and deciphering enzyme properties. 51 , 9772 (2015). Wong, T. S., Tee, K. L., Hauer, B. & Schwaneberg, U. Sequence saturation mutagenesis (SeSaM): a novel method for directed evolution. Nucleic Acids Res 32 , e26 (2004). https://doi.org:10.1093/nar/gnh028 Ruff, A. J., Kardashliev, T., Dennig, A. & Schwaneberg, U. The Sequence Saturation Mutagenesis (SeSaM) method. Methods Mol Biol 1179 , 45-68 (2014). https://doi.org:10.1007/978-1-4939-1053-3_4 Wong, T. S. et al. Transversion-enriched sequence saturation mutagenesis (SeSaM-Tv+): a random mutagenesis method with consecutive nucleotide exchanges that complements the bias of error-prone PCR. Biotechnol J 3 , 74-82 (2008). https://doi.org:10.1002/biot.200700193 Zhao, J., Kardashliev, T., Joelle Ruff, A., Bocola, M. & Schwaneberg, U. Lessons from diversity of directed evolution experiments by an analysis of 3,000 mutations. Biotechnol Bioeng 111 , 2380-2389 (2014). https://doi.org:10.1002/bit.25302 Zhao, J. et al. Amino acid substitutions in random mutagenesis libraries: lessons from analyzing 3000 mutations. 101 , 3187 (2017). Reetz, M. T., Carballeira, J. D. & Vogel, A. Iterative saturation mutagenesis on the basis of b factors as a strategy for increasing protein thermostability. 45 , 7751 (2006). Ahmad, S. & Rao, N. M. Thermally denatured state determines refolding in lipase: Mutational analysis. 18 , 1196 (2009). Sellés Vidal, L., Isalan, M., Heap, J. T. & Ledesma-Amaro, R. A primer to directed evolution: current methodologies and future directions. 4 , 291 (2023). Zhao, H. & Zha, W. In vitro 'sexual' evolution through the PCR-based staggered extension process (StEP). Nat Protoc 1 , 1865-1871 (2006). https://doi.org:10.1038/nprot.2006.309 Ahmad, S., Kamal, M. Z., Sankaranarayanan, R. & Rao, N. M. Thermostable Bacillus subtilis Lipases: In Vitro Evolution and Structural Insight. 381 , 340 (2008). Mala, J. G. S. T., Satoru. . Understanding Structural Features of Microbial LipasesaAn Overview. Analytical Chemistry Insights 3 , 9-19 (2008). Singh, B., Bulusu, G. & Mitra, A. Understanding the thermostability and activity of bacillus subtilis lipase mutants: Insights from molecular dynamics simulations. 119 , 409 (2015). Huang, X. & Zhou, H. X. Similarity and difference in the unfolding of thermophilic and mesophilic cold shock proteins studied by molecular dynamics simulations. 91 , 2463 (2006). Meharenna, Y. T. & Poulos, T. L. Using molecular dynamics to probe the structural basis for enhanced stability in thermal stable cytochromes P450. 49 , 6686 (2010). Purmonen, M., Valjakka, J., Takkinen, K., Laitinen, T. & Rouvinen, J. Molecular dynamics studies on the thermostability of family 11 xylanases. 20 , 559 (2007). Scholtz, C. N. P. J. M. A Helix Propensity Scale Based on Experimental Studies of Peptides and Proteins. Biophysical Journal 75 , 422-427 (1998). Kamal, M. Z. et al. In vitro evolved non-aggregating and thermostable lipase: Structural and thermodynamic investigation. 413 , 741 (2011). Autour, A. & Ryckelynck, M. Ultrahigh-throughput improvement and discovery of enzymes using droplet-based microfluidic screening. 8 , (2017). Case, D. A. et al. Amber 2022, University of California, San Francisco. (2022). Case, D. A. et al. AmberTools. J Chem Inf Model 63 , 6183-6191 (2023). Tian, C. et al. ff19SB: Amino-Acid-Specific Protein Backbone Paramters Trained against Quantum Mechanics Energy Surfaces in Solution. J Chem Theory Comput 16 , 528-552. Izadi, S. Anandakrishnan, R. & Onufriev, A. V. Building Water Models: A Different Approach. J Phys Chem Lett 5 , 3863– 3871 (2014). Krivov, G. G. Shapovalov, M. V. & Dunbrack, Jr. R. L. Improved prediction of protein side-chain conformations with SCWRL4. Proteins 77 , 778-795 (2009) Bedart C. Renault N. Chavatte P. Porcherie, A. Lachgar, A. Capron, M. & Farce A. SINAPs: A Software Tool for Analysis and Visualization of Interaction Networks of Molecular Dynamics Simulations. J Chem Inf Model 62 , 1425-1436 (2022) Additional Declarations Yes there is potential Competing Interest. David A. Weitz and Karla Milcic are co-founders of RevivBio Inc. Supplementary Files ExtendedData.docx SUPPLEMENTARY DATA SET 1 SupplementaryMaterials.docx Cite Share Download PDF Status: Under Review Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-5583255","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":405733836,"identity":"34cb4433-46c7-4c8f-bf53-d8b21842aa22","order_by":0,"name":"David Weitz","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA9klEQVRIiWNgGAWjYFACxgYGhgowiw0qksDAwINfS2MDwxnStACtYWwjRYvB8eb2hz/n1cnx9x9ge/Bzz2F58/YExgdv2/BoOXOwsUFy22FjiRsJ7IY9zw4bzjnzgNlwLh4tZjcSGxsMtx1IbLjBwCbBcyCNcYZEAps0Lz4t9x82NiTOqauff/4Am+SfA2n2QC3sv/FquQEMsYMNzAkGB4CG8xywSQTZwoxPi/2ZxMaZDccOG268kdgmLXPAJnkGz8NmyTnncGuRbD/+4OOPmjp5ufOHj0m+OSBhO4M9+eCHN2W4tSABUDJAZYyCUTAKRsEoIBcAAF6qVzy4efgxAAAAAElFTkSuQmCC","orcid":"https://orcid.org/0000-0001-6678-5208","institution":"Harvard University","correspondingAuthor":true,"prefix":"","firstName":"David","middleName":"","lastName":"Weitz","suffix":""},{"id":405733837,"identity":"bed842a6-8748-4213-8eb3-9e5d42c5692b","order_by":1,"name":"Xinge Zhang","email":"","orcid":"","institution":"Harvard Universty","correspondingAuthor":false,"prefix":"","firstName":"Xinge","middleName":"","lastName":"Zhang","suffix":""},{"id":405733838,"identity":"ddfa8231-96c3-4ced-bcd0-05567d664618","order_by":2,"name":"Anqi Chen","email":"","orcid":"","institution":"Harvard Universty","correspondingAuthor":false,"prefix":"","firstName":"Anqi","middleName":"","lastName":"Chen","suffix":""},{"id":405733839,"identity":"d417a6ed-bbcd-46ca-9fbb-e2e59c7e39fe","order_by":3,"name":"Liyin Chen","email":"","orcid":"","institution":"Harvard Universty","correspondingAuthor":false,"prefix":"","firstName":"Liyin","middleName":"","lastName":"Chen","suffix":""},{"id":405733840,"identity":"6afd1479-80c0-4af0-b8ce-445688e1afaf","order_by":4,"name":"Rohan Thakur","email":"","orcid":"","institution":"Harvard Universty","correspondingAuthor":false,"prefix":"","firstName":"Rohan","middleName":"","lastName":"Thakur","suffix":""},{"id":405733841,"identity":"f5b12edb-ff27-4ea4-ba70-d8489b395f8b","order_by":5,"name":"Perry Ellis","email":"","orcid":"","institution":"Harvard Universty","correspondingAuthor":false,"prefix":"","firstName":"Perry","middleName":"","lastName":"Ellis","suffix":""},{"id":405733842,"identity":"c354222f-fd65-4826-a4ec-f08a8196adb6","order_by":6,"name":"Yi Xiao","email":"","orcid":"","institution":"Harvard Universty","correspondingAuthor":false,"prefix":"","firstName":"Yi","middleName":"","lastName":"Xiao","suffix":""},{"id":405733843,"identity":"b65fa6f1-4c58-4379-984a-c3aff964a25a","order_by":7,"name":"Miloš Milčić","email":"","orcid":"","institution":"University of Belgrade-Faculty of Chemistry","correspondingAuthor":false,"prefix":"","firstName":"Miloš","middleName":"","lastName":"Milčić","suffix":""},{"id":405733844,"identity":"f0413553-79f0-46d8-88bc-a356dd66b694","order_by":8,"name":"Karla Milcic","email":"","orcid":"https://orcid.org/0000-0002-3672-4575","institution":"Harvard Universty","correspondingAuthor":false,"prefix":"","firstName":"Karla","middleName":"","lastName":"Milcic","suffix":""}],"badges":[],"createdAt":"2024-12-05 03:00:24","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-5583255/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-5583255/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":74631403,"identity":"18da0799-de24-4081-b247-005fbbd9bd43","added_by":"auto","created_at":"2025-01-24 07:52:32","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":209456,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eWorkflow of the Functional Screening and Construction of Thermotolerant BsLipA libraries \u003c/strong\u003e(a) Site-directed mutagenesis using overlapping primer pools. Each position in the gene is subjected to 19 possible chemically diverse substitutions, covering the entire gene. (b) All primer pools are combined in one site-mutagenesis reaction in equimolar ratio. In total 18 primer pools are used to cover the entire gene. The schematic illustration only shows 5 primer pools for simplicity. (c) The resultant plasmids are recombined through staggered extension PCR (StEP), redistributing the mutations across templates. (d) The final RSM library contains random mutations throughout the gene. (e) Primer pools can be combined in different numbers in separate site-directed mutagenesis reactions to limit mutation number.\u003c/p\u003e","description":"","filename":"image1.png","url":"https://assets-eu.researchsquare.com/files/rs-5583255/v1/2ee0afe50423a54f14081598.png"},{"id":74631404,"identity":"bbe19934-d888-4cab-ac70-29567994cbc2","added_by":"auto","created_at":"2025-01-24 07:52:32","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":153414,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAnalysis of Mutagenesis Library Diversity and Positional Distribution.\u003c/strong\u003e (a) Positional distribution of mutations for even-numbered primer pools in the 9×2 library, illustrating precision and controllability. (b) Positional distribution of mutations for odd-numbered primer pools in the 9×2 library. (c) Positional distribution of mutations in the 9×2 library after StEP. (d) Mutation distribution of RSM libraries and error-prone libraries.\u003c/p\u003e","description":"","filename":"image2.png","url":"https://assets-eu.researchsquare.com/files/rs-5583255/v1/15cc2e4e2f0d1f7c47693a70.png"},{"id":74631008,"identity":"95736472-a33e-454a-ab36-df37d410b76c","added_by":"auto","created_at":"2025-01-24 07:44:32","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":245581,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eLibrary Characterization of RSM and Error-Prone Libraries.\u003c/strong\u003e (a) Percentage of silent mutations. (b) Percentage of full-length genes that contain premature stop codons. (c) Percentage of genes with indels (insertions and deletions ≤ 3 bp). The RSM libraries are shown in shades of blue color and Error-prone libraries are shown in different shades of orange color. (d) Distribution of the frequency of each mutation for the 361 distinct amino acid substitutions, arranged in ascending order. The red dashed line represents the average expected number of mutations.\u003c/p\u003e","description":"","filename":"image3.png","url":"https://assets-eu.researchsquare.com/files/rs-5583255/v1/814b957d21a6596afba0e1ec.png"},{"id":74631010,"identity":"dad1643d-d192-4cbe-aac4-21a4c3d34a45","added_by":"auto","created_at":"2025-01-24 07:44:32","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":992215,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDistinct Amino Acid Substitution (DAAS) and Categorical Accessibility in EPhigh and 9×2 RSM Libraries \u003c/strong\u003e(a) Heatmap depicting distinct amino acid substitutions (DAAS) in the EPhigh library. Rows correspond to the 19 wild-type amino acids in BsLipA, and columns represent the 20 possible substitutions. The color scale indicates the frequency of each substitution, with warmer colors representing higher frequencies. The numbers in each cell represent normalized counts of the substitutions. White-bordered squares highlight substitutions that can occur through a single nucleotide change in the wild-type codon. The data reflects approximately 72,000 sampled single amino acid substitutions and is normalized by the wild-type amino acid usage in BsLipA. (b) Heatmap of distinct amino acid substitutions in the 9×2 library, similar to (a), showing the normalized substitution frequencies. (c) Heatmap for categorical accessibility in the EPhigh library. The wild-type amino acids are grouped into eight categories based on their side chain properties: non-polar (NP: M, I, L, V, A), polar uncharged (PU: S, T, Q, N), positively charged (PC: R, K), negatively charged (NC: D, E), aromatic (Ar: F, W, Y), proline (P), cysteine (C), and glycine (G). Each wild-type amino acid is evaluated for accessibility to these categories, showing that some properties are less accessible to certain residues due to their codon constraints. (d) Heatmap for categorical accessibility in the 9×2 library, following the same classification scheme as (c), illustrating the broader categorical accessibilities or limitations at individual residue positions across the 9×2 RSM library.\u003c/p\u003e","description":"","filename":"image4.png","url":"https://assets-eu.researchsquare.com/files/rs-5583255/v1/929f0c92be0922b0c7af9790.png"},{"id":74631016,"identity":"8233a4b4-a7a6-4a8a-a53b-65217d907089","added_by":"auto","created_at":"2025-01-24 07:44:32","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":632586,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCharacterization of Thermotolerant BsLipA Variants.\u003c/strong\u003e (a) Schematic representation of the droplet-based high-throughput screening process. \u003cem\u003eE. coli\u003c/em\u003e cells expressing BsLipA library variants are encapsulated in droplets with resorufin acetate substrate and BugBuster Protein Extraction Reagent. After incubation at 65°C for 1 hour, droplets are sorted based on fluorescent signals indicating enzyme activity. (b) Mutations identified in MT5 and MT20 by Sanger sequencing. (c) T\u003csub\u003e50\u003c/sub\u003e measurement of the WT, MT5, and MT20. Activity is measured after 20 minutes of incubation at various temperatures (25-95°C). The WT exhibits a T\u003csub\u003e50\u003c/sub\u003e of 50.86°C, while MT5 and MT20 have T\u003csub\u003e50\u003c/sub\u003e values greater than 95°C, indicating significantly enhanced thermal stability. (d) NanoDSF analysis shows the unfolding ramps of the WT, MT5, and MT20. (e) NanoDSF analysis shows the refolding ramps of the WT, MT5, and MT20. (f) Thermal unfolding mechanisms of the WT, where N: native state(s); U: monomeric unfolded state(s); A: aggregate(s); I: intermediate denatured state(s). (g) Thermal unfolding mechanisms of MT5 and MT20. (h) Protein solubility upon heat incubation for WT, MT5, and MT20.\u003c/p\u003e","description":"","filename":"image5.png","url":"https://assets-eu.researchsquare.com/files/rs-5583255/v1/2e77c85b547a971a90f4bf57.png"},{"id":74631011,"identity":"e1ace5f1-5381-4142-8161-e9df8b5e095f","added_by":"auto","created_at":"2025-01-24 07:44:32","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":637027,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eStructure of wild-type BsLipA and identified variants.\u003c/strong\u003e (a) Three dimensional structure of wild-type BsLipA (PDB ID: 5CRI). Catalytic triad residues are labeled with pink color, residues mutated in both variants 5 and 20 are shown in dark blue, residues mutated in variant 5 are shown in purple and residues mutated in MT20 are shown in green (b) Aligned representative 3D structures of MT20 (light blue) and wild-type (beige) showing new hydrogen bonds formed in MT20 (dark blue). Catalytic residues are labeled in pink color. (c) Representative three dimensional structure of wild-type BsLipA generated from molecular dynamics 200-ns production run at 25°C showing hydrophobicity surface of the protein. Catalytic triad residues are labeled with pink color. Hydrophobic surfaces are labeled as gradient of red color and hydrophilic surfaces are labeled as a gradient of blue color, where the in hydrophobicity/hydrophilicity correlates with the increase in color intensity. (d) Representative three dimensional structure of wild-type BsLipA generated from molecular dynamics 200-ns production run at 100°C showing hydrophobicity surface of the protein. Catalytic triad residues are labeled with pink color. Hydrophobic surfaces are labeled as gradient of red color and hydrophilic surfaces are labeled as a gradient of blue color, where the in hydrophobicity/hydrophilicity correlates with the increase in color intensity. (e) Aligned representative 3D structures of wild-type BsLipA generated from the first (beige) and last 20 ns (light blue) of 200-ns molecular dynamics simulation production runs performed at 100°C. Catalytic residues are shown in pink, residues mutated in both variants 5 and 20 are shown in dark blue and residue mutated in MT20 is shown in green. (f) Aligned representative 3D structures of MT20 generated from the first (beige) and last 20 ns (light blue) of 200-ns molecular dynamics simulation production runs performed at 100°C. Catalytic residues are shown in pink, residues mutated in both variants 5 and 20 are shown in dark blue and residue mutated in MT20 is shown in green. (g) Representative three dimensional structure of MT20 generated from molecular dynamics 200-ns production run at 25°C showing hydrophobicity surface of the protein. Catalytic triad residues are labeled with pink color. Hydrophobic surfaces are labeled as gradient of red color and hydrophilic surfaces are labeled as a gradient of blue color, where the in hydrophobicity/hydrophilicity correlates with the increase in color intensity. (h) Representative three dimensional structure of MT20 generated from molecular dynamics 200-ns production run at 100°C showing hydrophobicity surface of the protein. Catalytic triad residues are labeled with pink color. Hydrophobic surfaces are labeled as gradient of red color and hydrophilic surfaces are labeled as a gradient of blue color, where the in hydrophobicity/hydrophilicity correlates with the increase in color intensity.\u003c/p\u003e","description":"","filename":"image6.png","url":"https://assets-eu.researchsquare.com/files/rs-5583255/v1/4214ac204bc0bdf7b59b431c.png"},{"id":74633082,"identity":"d1821ee2-fdf1-4309-bb2e-be5a107408eb","added_by":"auto","created_at":"2025-01-24 08:08:44","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":3901677,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-5583255/v1/5f046936-1eb1-4bdf-9e0a-2b87cd0e3715.pdf"},{"id":74631004,"identity":"f77b16cc-55b2-495a-8396-89042613d777","added_by":"auto","created_at":"2025-01-24 07:44:32","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":503173,"visible":true,"origin":"","legend":"SUPPLEMENTARY DATA SET 1","description":"","filename":"ExtendedData.docx","url":"https://assets-eu.researchsquare.com/files/rs-5583255/v1/7f2a8833faad9f1afafbb364.docx"},{"id":74631405,"identity":"832f44a1-de4e-4c8f-9b80-9d3062858f67","added_by":"auto","created_at":"2025-01-24 07:52:32","extension":"docx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":4141403,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryMaterials.docx","url":"https://assets-eu.researchsquare.com/files/rs-5583255/v1/110f78cbec81b0a8ee8e6d2b.docx"}],"financialInterests":"\u003cb\u003eYes\u003c/b\u003e there is potential Competing Interest.\nDavid A. Weitz and Karla Milcic are co-founders of RevivBio Inc.","formattedTitle":"Random Saturation Mutagenesis to Generate Highly Diverse Libraries for Directed Evolution","fulltext":[{"header":"Introduction","content":"\u003cp\u003eEnzymes are biological catalysts that facilitate biochemical reactions with efficiency and specificity. However, their application in industrial processes often demands adaptations to maintain stability and catalytic performance under challenging conditions. These adaptations can be achieved through a process called directed evolution. Directed evolution comprises two steps: introducing mutations to create a mutagenesis library and screening for variants with desired functional improvements. For a typical-sized enzyme, the sequence space is vast, far exceeding the throughput of existing screening methods. Therefore, when a library is constructed, a subspace in the sequence space is determined, and the choice of the sequence subspace defines the range for the screen. Mutagenesis libraries are designed to maximize the chance of discovering improved variants. For well-studied enzymes, mutagenesis libraries can be designed rationally based on prior structural and mechanistic information. For enzymes lacking sufficient prior information, random mutagenesis is used. Random mutagenesis introduces mutations across the entire sequence, allowing exploration of a broader sequence space\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eTraditionally, random mutagenesis is achieved using error-prone polymerases, either \u003cem\u003ein vitro\u003c/em\u003e through polymerase chain reactions (PCRs) or \u003cem\u003ein vivo\u003c/em\u003e through continuous evolution methods\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e,\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e. Error-prone polymerases are strongly biased towards single nucleotide changes per codon, resulting in a high frequency of silent and redundant mutations. Moreover, with a single nucleotide mutation per codon, only 5.7 amino acid substitutions are accessible, on average, to each wild-type amino acid\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e\u003c/sup\u003e. Furthermore, even these substitutions are often physiochemically conserved due to the evolutionary organization of codons to preserve activity\u003csup\u003e\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e. This narrows the chemical diversity to an evolutionarily favored subspace and results in limited deviations from the wild-type function; many amino acid substitutions leading to large functional improvements can be excluded from the library\u003csup\u003e\u003cspan additionalcitationids=\"CR7 CR8 CR9\" citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e\u003c/sup\u003e. Sequence Saturation Mutagenesis (SeSaM) pioneers the efforts to improve chemical diversity in random mutagenesis by allowing for some consecutive nucleotide exchanges\u003csup\u003e\u003cspan additionalcitationids=\"CR12 CR13\" citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e. In SeSaM, DNA templates are randomly fragmented; the fragments are tailed with one or multiple consecutive universal bases, and then used as primers to synthesize full-length DNA. These products are further amplified by PCR where universal bases are replaced by one of the four standard bases, introducing mutations. This workflow is rather complex and time consuming, yet still approximately 70% of the mutations remain single nucleotide substitutions\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e, and only 35% of all amino acid substitutions\u003csup\u003e\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e are enabled. To achieve greater chemical diversity at specific amino acid positions, site saturation mutagenesis (SSM) is often employed. This approach allows for the systematic introduction of all possible amino acid substitutions at specific sites. However, these sites are chosen either by the screening outcome of random mutagenesis libraries\u003csup\u003e\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u003c/sup\u003e, where important residues may be missed because they were not saturated\u003csup\u003e\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e,\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e, or by mechanistic hypotheses based on prior knowledge, reducing the chance to discover unexpected mechanisms\u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e. Therefore, there exists the need for a random mutagenesis method that can achieve high chemical diversity across the entire gene without requiring prior knowledge. A method that introduces all distinct amino acid substitutions randomly throughout an enzyme sequence, while controlling the mutation rate, would allow for a more comprehensive exploration of the sequence space.\u003c/p\u003e \u003cp\u003eIn this paper, we present Random Saturation Mutagenesis (RSM), a method that introduces mutations across the entire protein sequence, evenly sampling all chemically distinct conversions at a controlled mutation rate. This method uses multiple pools of overlapping mutagenic primers targeting subregions of the gene, allowing for the control of mutation rates by varying the number of primer pools used in each reaction. Sequencing results confirm that RSM libraries achieve controlled mutation rates, high variant quality, and high chemical diversity. To demonstrate the utility of our library preparation method, we use RSM to drastically improve the thermotolerance of \u003cem\u003eBacillus subtilis\u003c/em\u003e Lipase A (BsLipA). We construct an RSM mutagenesis library and use droplet screening to identify two vastly improved variants after only one generation of mutagenesis. Each of the variants has seven single amino acid substitutions and shows over 40\u0026deg;C increase in T\u003csub\u003e50\u003c/sub\u003e values, the temperature at which they retain 50% of their activity after heat incubation at that temperature. To understand the roles of the newly discovered amino acid substitutions in improving thermotolerance, we conduct molecular dynamics simulations. The identified key substitutions that result in large T\u003csub\u003e50\u003c/sub\u003e improvements are not obtainable by single nucleotide changes in a trinucleotide codon.\u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eConstruction of RSM Libraries\u003c/h2\u003e \u003cp\u003eTo demonstrate the utility of Random Saturation Mutagenesis (RSM) to create a random library that is highly chemically diverse, we choose the enzyme BsLipA as a candidate to improve its thermotolerance. Thermotolerance of enzymes can be improved by mutations scattered across the entire gene\u003csup\u003e\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u003c/sup\u003e; therefore, random mutagenesis is particularly well-suited for this problem. When there is no prior information about the structure or function, an ideal random mutagenesis library should be as chemically diverse as possible to facilitate a comprehensive and unbiased exploration of the sequence subspace\u003csup\u003e\u003cspan additionalcitationids=\"CR7 CR8 CR9\" citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e\u003c/sup\u003e. We use site-directed mutagenesis primers to introduce mutations. To ensure annealing of primers during mutagenesis, we choose a primer length of 30 nucleotides, covering a 10-residue region as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea. We divide the 540 bp BsLipA gene into 18 regions, each targeted by a set of fully overlapping 30-nucleotide primers, as shown in 1b. For simplicity, the schematic illustrates only 5 primer pools, each covering 5 residues. For each residue, we design a set of 19 mutagenesis primers to introduce all possible single amino acid substitutions, saturating each residue. Each primer is complementary to the target gene except for the three nucleotides encoding the amino acid substitution. We select a single codon for each unique amino acid substitution to ensure even distribution at the mutated site. Moreover, the explicit choice of 19 codons prevents the introduction of rare codons and stop codons. The primers targeting all 10 residues in each region are purchased as a primer pool of 190 one-directional primers (10 positions x 19 substitutions per position). To cover the full gene, 18 such primer pools are purchased.\u003c/p\u003e \u003cp\u003eTo introduce mutations throughout the gene, we combine all 18 primer pools into a single mutagenesis reaction using the wild-type plasmid as the template. During a mutagenesis reaction, only one-directional primers are used. The polymerase amplifies the gene from the 3\u0026rsquo; end of the first primer to the 5\u0026rsquo;end of the next primer leaving a nick in between. A ligase fixes this nick to form a double-stranded plasmid incorporating the designed mutations in one strand; the wild-type strand is then digested by DpnI. In this reaction, a maximum of 18 primers can anneal to the template DNA simultaneously. We transform the single-stranded DNA product from the mutagenesis reaction into bacterial hosts to synthesize double-stranded plasmids. At this stage, the number of mutations per region is limited to one, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ec. To increase mutation diversity, we use staggered extension PCR (StEP)\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e. In a StEP reaction, the extension step is set deliberately to be short, synthesizing only partial genes per cycle. During the denaturing step of the next cycle, these incomplete strands unbind from their original templates and anneal to different ones, continuing extension and effectively recombining mutations from various templates. The result is a final RSM library with random mutations throughout the gene, as illustrated for 5 primer pools in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ed. We sequence the library using PacBio technology to obtain high-quality full-length reads and analyze approximately 10,000 of these by aligning them to the wild-type BsLipA gene to identify mutations. On average, this RSM library contains 7.85 single amino acid substitutions per gene.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eThe high mutation rate obtained in this library can result in a predominant population of inactive mutants, which may not be suitable for standard directed evolution applications. To achieve controlled mutation rates, we modify the RSM workflow by using only a subset of all primer pools for each mutagenesis reaction, and performing multiple separate reactions to cover the whole gene. This approach limits the number of primers that can anneal to the DNA template in a single reaction. For example, we can prepare two separate mutagenesis reactions, each containing an equimolar mix of 9 primer pools, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ee. Reaction 1 uses all odd-numbered primer pools and reaction 2 uses all even-numbered primer pools. To confirm the mutations introduced in either reaction, we sequence the product from these two reactions and count the amino acid substitutions at each residue position. To show the distribution of substitutions at each position, we compile a heatmap of which each column represents a residue in the BsLipA sequence, and each row shows the mutations introduced. In both reaction 1 and 2, the mutations are confined within the boundaries of the designed primer pools as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea and \u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb. To remove such confinement and to increase diversity of the final RSM library, the products from the two reactions are recombined through StEP. The final product of this library exhibits uniform coverage of single amino acid substitutions across the gene as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec. The average number of single amino acid substitutions of this library is measured to be 4.64 by PacBio sequencing. Since we use 9 primer pools in 2 mutagenesis reactions, we label this library \u0026ldquo;9\u0026times;2\u0026rdquo;, where 9 denotes the number of pools used per reaction, and 2 denotes the number of reactions. Similarly, the previous library is labeled \u0026ldquo;18\u0026times;1\u0026rdquo;.\u003c/p\u003e \u003cp\u003eFollowing the same procedure, we create a \u0026ldquo;6\u0026times;3\u0026rdquo; library by using 6 primer pools in 3 separate reactions, and a \u0026ldquo;4\u0026times;3\u0026thinsp;+\u0026thinsp;3\u0026times;2\u0026rdquo; library by combining 3 or 4 primer pools per reaction. The average mutation rates for these libraries are measured by PacBio sequencing to be 5.09 and 1.97 amino acid substitutions per protein sequence, respectively. Moreover, the RSM approach can be adapted to construct single saturation mutagenesis libraries by using only one primer pool in 18 separate mutagenesis reactions, without the subsequent StEP recombination (Extended Data Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). The overlapping primer design ensures that the annealing of a primer completely prevents other primers in this region from binding. The average mutation rate for this single saturation library is confirmed to be 0.99 measured by PacBio sequencing, and we label this library as \u0026ldquo;Single\u0026rdquo;. The mutation distribution is shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ed.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eBy performing RSM with varying numbers of primer pools per reaction, mutations can be introduced at controlled rates and random positions throughout the gene. While existing methods like error-prone PCR can achieve similar controlled mutation rates across the gene, by comparison to RSM, they produce low-quality and low-diversity libraries. To illustrate this comparison, we construct three error-prone mutagenesis libraries using a commercially available mutazyme kit (GeneMorph II, Agilent), rendering libraries with high (EPhigh, 2.70 substitutions per protein sequence), medium (EPmid, 2.04), and low (EPlow, 1.18) mutation rates. For a variant in a library to be improved, it must translate to the full length and must differ from the wild-type protein sequence. Silent mutations do not contribute to the diversity of the library. All RSM libraries consistently exhibit less than ~\u0026thinsp;5% silent mutations, significantly lower than the error-prone libraries, each of which contains\u0026thinsp;~\u0026thinsp;30% silent mutations, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb. Additionally, variants that translate into truncated proteins cannot function, hence their presence also negatively impacts library quality. Such truncations can occur when a trinucleotide codon is mutated to either TAG, TGA, or TAA, which are stop codons signaling the termination of the translation process. In all RSM libraries, premature stop codons are found to exist in less than 1% of all variants, much lower than in error-prone libraries, where 5%-14% of variants contain premature stop codons, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea. Similarly, insertions and deletions (indels) of nucleotides can disrupt protein function by causing frameshifts that lead to incorrect and truncated amino acid sequences if they do not occur in multiples of three. When comparing an RSM library to an error-prone library with similar mutation rates, such as single versus EPlow and 4\u0026times;3\u0026thinsp;+\u0026thinsp;3\u0026times;2 versus EPmid, the proportion of variants that contain indels in RSM libraries is significantly lower than in error-prone libraries, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ec. Moreover, in error-prone libraries, the proportion of variants with premature stop codons and with indels both increase with mutation rate. Thus, if a mutation rate even higher than EPhigh is required, the proportion of truncated variants will inevitably increase further. By contrast, in RSM libraries, the proportions of these incorrect and truncated variants are not only low but also independent of the mutation rate; this demonstrates the significant advantage of the RSM libraries over the error-prone libraries particularly when high mutation rates are required, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb and \u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ec.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWild-type BsLipA contains 19 of the 20 standard amino acids (excluding cysteine), each of which can be mutated into 19 other possible amino acids, resulting in 361 possible distinct amino acid substitutions (DAASs). An ideal random mutagenesis method for our application of interest should uniformly sample DAASs to ensure high chemical diversity of the library. To evaluate the uniformity of sampling of these 361 DAASs by RSM and by error-prone PCR, we choose an RSM library, 9\u0026times;2, and an error-prone library, EPhigh, with similar mutation rates and sample the same number of amino acid substitutions from each library. We count the occurrences for each DAAS and rank them by abundance to visualize the uniformity of sampling in each library. The 9\u0026times;2 RSM library samples DAASs much more uniformly than the EPhigh library. In 9\u0026times;2, all 361 DAASs are enabled with a mean of 21 and a standard deviation of 10. By contrast, in EPhigh, the top ~⅓ most abundant DAASs comprise over 99.5% of the total count, whereas the remaining ~⅔ DAASs each occur less than or equal to once, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ed. While the mean of DAAS occurrence is the same as 9\u0026times;2, the standard deviation is 38, four times that of 9\u0026times;2. To evaluate the deviation of either library from an ideal uniform distribution, we calculate the Gini coefficient for the normalized counts of DAASs. The 9\u0026times;2 RSM library exhibits a Gini coefficient of 0.27, lower than 0.78 for EPhigh, indicating a much more uniform distribution of DAASs in the RSM library (Extended Data Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eTo determine the source of the bias in EPhigh, we investigate the codon change required for each DAAS. We arrange the counts of DAASs in a 20x19 table where each column represents an amino acid used in the wild-type BsLipA, and each row indicates an observed substitution from the wild-type amino acid. The color reflects the counts of individual DAASs observed in sequencing. The color-annotated table shows the much lower variation in the distribution of DAASs sampled in the RSM library as compared to the EPhigh library, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea,b. Moreover, in the EPhigh library, the top ⅓ most abundant substitutions, which account for \u0026gt;\u0026thinsp;99.5% of all occurrences, are precisely those enabled by single nucleotide changes per codon, as highlighted by the white frames in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea. Certain amino acid substitutions occur with anomalously high frequency, such as the W-to-R substitution and the F-to-L substitution, likely due to the intrinsic bias of the error-prone polymerase. By contrast, the DAASs outside of the white frames, which require multiple nucleotide changes, are only detected at most once; hence, they are effectively inaccessible by error-prone PCR. On average, each amino acid only has access to 6.6 DAASs, and these DAASs cannot enable all categories of sidechain physicochemical properties. To access the accessibility of each wild-type amino acid to different sidechain properties, we classify all 20 amino acids into eight categories: non-polar (NP, M, I, L, V, A), polar uncharged (PU, S, T, Q, N), positively charged (PC, R, K, L), negatively charged (NC, D, E), aromatic (Ar, F, T, Y), and three special groups, P, C, G, each consisting of only one amino acid due to their distinct properties. Through error-prone PCR, each wild-type amino acid has a few inaccessible property categories, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ec. Furthermore, in this analysis, we aggregate the same wild-type amino acids at different positions which are encoded by different codons. Therefore, when each residue position is individually considered, the categorical accessibility of each wild-type residue is even more limited. In contrast, the 9x2 RSM library uniformly samples all possible DAASs, and thereby enables all chemical properties at each residue, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb, d. By uniformly sampling all DAASs, RSM enables much more chemical diversity in the library than traditional random mutagenesis methods, enhancing the likelihood to discover dramatically improved variants.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eFunctional Screening of RSM library\u003c/h3\u003e\n\u003cp\u003eTo illustrate the potential of RSM, we evolve BsLipA to substantially improve its thermotolerance. Lipases such as BsLipA are versatile enzymes that catalyze the hydrolysis and synthesis of esters, playing crucial roles in various industrial applications, including biodiesel production, and use in detergents. However, their practical utility is often limited by the high temperatures used in industrial processes (Vivek, 2022). The wild-type BsLipA exhibits a relatively low melting temperature (T\u003csub\u003em\u003c/sub\u003e), and irreversibly loses activity upon heating. Thus, the development of lipases resistant to heat inactivation will significantly expand their industrial applications, enabling their use across a broader range of high-temperature processes.\u003c/p\u003e \u003cp\u003eTo make significant improvement in the thermotolerance of BsLipA, we maximize the diversity of our library in two ways: we use RSM to ensure high chemical diversity at each residue position, and choose a high mutation rate so each pair of variants, on average, are different by more residues. However, as the mutation rate of a library increases, the theoretical sequence space grows rapidly, and the best variants can become too rare to be discovered by the experimentally practical screening throughput. With the screening throughput of droplet sorting, which is ~\u0026thinsp;10\u003csup\u003e7\u003c/sup\u003e, there exists a maximum mutation rate with which the best mutants can still be discovered. To estimate this maximum mutation rate for our problem, we assume that the most thermotolerant mutants in a library can result from the combination of single mutations that individually improve thermotolerance. The frequency of such variants in a high mutation rate library can then be estimated by the proportion of single mutants that exhibit higher thermal stability than the wild-type. This proportion can be inferred by constructing a single saturation library using RSM and performing a droplet cytometry analysis on this library. For BsLipA, a deep mutational scan of all single mutants has been done to measure the change in thermotolerance over the wild-type\u003csup\u003e\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u003c/sup\u003e. Therefore, we count the single mutants in this dataset that show improved thermotolerance over the wild-type, which comprises\u0026thinsp;~\u0026thinsp;0.1 of all single mutants. Using this proportion, we estimate the frequency of mutants in an n-mutation RSM library that combine all individually beneficial mutations to be (0.1)\u003csup\u003en\u003c/sup\u003e. For more than one of such mutants to be discovered, we require this frequency to be higher than one over the screening throughput of 10\u003csup\u003e7\u003c/sup\u003e. Therefore, we choose a mutation rate of ~\u0026thinsp;7.\u003c/p\u003e \u003cp\u003eTo identify BsLipA variants that retain activity after thermal exposure, we express the RSM library in \u003cem\u003eE. coli\u003c/em\u003e host and use a fluorescence-activated droplet sorting (FADS) assay to screen for thermotolerant candidates. Each \u003cem\u003eE. coli\u003c/em\u003e cell is co-encapsulated with lysis reagents and the fluorogenic substrate, resorufin acetate, in a 20 \u0026micro;m droplet. The droplets are incubated at 65\u0026deg;C for one hour, during which time the cells are lysed and the BsLipA proteins are released. Variants that remain active at 65\u0026deg;C catalyze the conversion of resorufin acetate into the fluorescent product, resorufin. After incubation, the droplets are cooled to 25\u0026deg;C and reinjected to a sorter for high-throughput screening, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ea. The sorting experiment, conducted at 25\u0026deg;C, can also capture variants that might regain activity during the cooling phase, further converting resorufin acetate. We screen\u0026thinsp;~\u0026thinsp;10\u003csup\u003e7\u003c/sup\u003e droplets and sort the 0.04% of droplets with the highest fluorescence intensity. The sorted droplets are merged, and the DNA is amplified and used as input for a second round of cell sorting; the repeated sort reduces the effects of variations in cell expression. In the second round, we screen\u0026thinsp;~\u0026thinsp;3x10\u003csup\u003e6\u003c/sup\u003e droplets and sort the 0.04% of droplets. To isolate individual variants, we transform the sorted plasmids into \u003cem\u003eE. coli\u003c/em\u003e and plate the cells on agar. Single colonies are then picked and grown in separate wells of a 96-well plate. To coarsely assess the thermotolerance of the sorted variants, we lyse the cells using lysis reagents, incubate the lysates at 80\u0026deg;C for one hour, and then perform a fluorescence-based assay on a plate reader to measure residual activity at 25\u0026deg;C from the cell lysates (Fig. S6.). We identify two variants with highest residual activity, named mutant 5 (MT5) and mutant 20 (MT20). Both MT5 and MT20 contain seven amino acid substitutions, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eb. Each single amino acid substitution, except for the S163P substitution in MT20, individually improves thermotolerance\u003csup\u003e\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u003c/sup\u003e, consistent with the assumption that very good mutants can result from the combination of beneficial single substitutions. These two variants as well as the wild-type (WT) are subsequently purified using NiNTA mini-spin columns for further characterization.\u003c/p\u003e\n\u003ch3\u003eCharacterization of Thermotolerant BsLipA Variants\u003c/h3\u003e\n\u003cp\u003eTo evaluate the improvement of heat tolerance of MT5 and MT20, we incubate the purified enzymes at various temperatures for 20 minutes and measure the residual activity after cooling them back to room temperature. To enable comparison, we normalize all residual activity values to the WT activity after incubation at 25\u0026deg;C. The WT exhibits constant residual activity after incubation at temperatures from 25\u0026deg;C to 45\u0026deg;C. However, at higher temperatures, the residual activity declines rapidly and remains low from 55\u0026deg;C to 95\u0026deg;C, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ec. To characterize this sharp drop, we interpolate the data to estimate T\u003csub\u003e50\u003c/sub\u003e, the temperature at which 50% activity is lost as compared to the lowest incubation temperature of 25\u0026deg;C. The value of T\u003csub\u003e50\u003c/sub\u003e is commonly used to compare the thermotolerance of enzymes\u003csup\u003e\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e,\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e,\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u003c/sup\u003e. For the WT, T\u003csub\u003e50\u003c/sub\u003e\u0026thinsp;\u0026asymp;\u0026thinsp;51.5\u0026deg;C, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ec. In contrast, MT5 maintains high residual activity over a wide temperature range, showing only a gradual decrease at the highest temperatures as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ec. Similarly, MT20 maintains high residual activity up to 75\u0026deg;C, with a decline of approximately 35% at 95\u0026deg;C as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ec. Thus, the T\u003csub\u003e50\u003c/sub\u003e values of both variants are beyond the range of our measurement. In addition, MT20 also demonstrates nearly two-fold higher enzymatic activity than WT at room temperature.\u003c/p\u003e \u003cp\u003eTo understand the origin of the different thermotolerance of the WT and the variants, we measure the structural stability at elevated temperatures using Nano Differential Scanning Fluorimetry (NanoDSF). During a temperature ramp from 20\u0026deg;C to 95\u0026deg;C, we monitor the red shift of intrinsic fluorescence by measuring the change in R\u003csub\u003ef\u003c/sub\u003e, the fluorescence intensity ratio between 350 nm and 330nm. An increase in R\u003csub\u003ef\u003c/sub\u003e reflects the increase of polarity in the chemical microenvironment of tyrosines and tryptophans, thereby indicating unfolding of the protein. For the WT, we observe a single inflection in the temperature dependence of R\u003csub\u003ef\u003c/sub\u003e indicating an unfolding transition at T\u003csub\u003em\u003c/sub\u003e = 54.7\u0026deg;C, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ed. The value of T\u003csub\u003em\u003c/sub\u003e of the WT roughly coincides with the value of T\u003csub\u003e50\u003c/sub\u003e, suggesting that the enzyme activity is irreversibly lost following thermal unfolding. This behavior is attributed to irreversible aggregation upon thermal unfolding\u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e; this behavior is irreversible, as shown by the NanoDSF, where R\u003csub\u003ef\u003c/sub\u003e exhibits no defined transition upon cooling from 95\u0026deg;C to 20\u0026deg;C, indicating the absence of refolding, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ee. As the temperature approaches the T\u003csub\u003em\u003c/sub\u003e of the WT, its native tertiary structure progressively unfolds. Unlike the native state (N) where the hydrophobic residues are buried within the protein structure, certain intermediate denatured states (I) expose the hydrophobic residues on the surface and lead to irreversible aggregation (A), where the aggregates remain stable even when the temperature decreases to room temperature. This irreversible aggregation inhibits the refolding of the native structure and prevents the recovery of enzymatic activity, as shown schematically in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ef. Therefore, the heat inactivation of the WT involves two key steps: thermal unfolding and irreversible aggregation of unfolded states.\u003c/p\u003e \u003cp\u003eSince MT5 and MT20 both show exceptional heat tolerance as compared to the WT, they must counteract at least one of the two steps in the pathway of the WT heat inactivation. We first examine using NanoDSF if the thermal unfolding of the variants is shifted to much higher temperatures. Both variants exhibit one thermal transition, with T\u003csub\u003em\u003c/sub\u003e values of 61.2\u0026deg;C for MT5 and 57.4\u0026deg;C for MT20, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ed. The increase in T\u003csub\u003em\u003c/sub\u003e demonstrates the higher structural stability of the variants as compared to the WT. However, this moderate improvement in structural stability is not sufficient to explain the \u0026gt;\u0026thinsp;40\u0026deg;C difference in T\u003csub\u003e50\u003c/sub\u003e between the variants and the WT. Unlike the WT, which irreversibly aggregates and loses activity after incubation around or above its T\u003csub\u003em\u003c/sub\u003e, both variants can experience temperatures\u0026thinsp;\u0026gt;\u0026thinsp;30\u0026deg;C above their T\u003csub\u003em\u003c/sub\u003e values, and still recover more than half of their activity upon temperature reduction. Therefore, the variants must disturb the second step of the WT heat inactivation process. At temperatures above T\u003csub\u003em\u003c/sub\u003e, where the native structure is lost, the variants likely adopt mechanisms to avoid aggregation-prone conformations and remain monomeric. The structural reversibility of both variants is verified by the defined transitions in their NanoDSF refolding ramps, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ee. By circumventing the stable aggregated state, they retain the ability to refold into functional structures upon temperature reduction, as shown schematically in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eg. To confirm the aggregation resistance of the variants, we incubate both variants at 65\u0026deg;C for 10 minutes, centrifuge each sample to remove the aggregated fraction and measure the protein concentration that is maintained in each solution. Each variant retains\u0026thinsp;\u0026gt;\u0026thinsp;90% protein mass in the supernatant. In contrast, the WT retains only\u0026thinsp;~\u0026thinsp;50% protein mass in the supernatant, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eh. These characterizations confirm two mechanisms for heat tolerance of MT5 and MT20: improved thermal stability of the native structure, and reduced propensity for aggregation in thermally unfolded conformations. The first mechanism accounts for the T\u003csub\u003em\u003c/sub\u003e improvements of 6.5\u0026deg;C and 2.7\u0026deg;C, whereas the second mechanism is the major contribution to the \u0026gt;\u0026thinsp;40\u0026deg;C improvements in T\u003csub\u003e50\u003c/sub\u003e.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e\n\u003ch3\u003eMolecular Dynamics Simulations\u003c/h3\u003e\n\u003cp\u003eTo understand the role of these substitutions on the thermotolerance, we perform molecular dynamics (MD) simulations for the WT and for MT20. The catalytic triad of BsLipA is formed by S77, D133, and H156 located in flexible loops enclosing the active pocket, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea. Both thermally stable variants contain mutations near the catalytic residues 133 and 156, as well as mutations at position 17, located on the loop considered to mimic a typical lipase lid absent in BsLipA\u003csup\u003e\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e\u003c/sup\u003e, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea. For the WT, structure 5CRI from Protein Data Bank (PDB) is used as the starting structure and equilibrated at 25\u0026deg;C in the MD simulations. To obtain a starting structure for the MT20, we introduce the seven amino acid substitutions into the equilibrated WT structure using SCWRL4 and re-equilibrate it. We perform 200 ns of MD simulations at 25\u0026deg;C for both the WT and the MT20 and use trajectory information to create representative 3D structures. The representative structure of MT20 closely resembles that of the WT as indicated by a high degree of structural overlap, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb. This similarity indicates that the mutations do not significantly alter the overall protein fold. However, by aligning the WT and MT structures and visualizing interactions present only in MT20, we identify new hydrogen bonds stabilizing the MT20 structure, as indicated by light blue lines in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb. These include a hydrogen bond between V154 and S130 each next to a catalytic residue, stabilizing the catalytically active conformation, and another hydrogen bond between K170 and the mutated S174 residue, stabilizing the terminal helix. These interactions may contribute to the thermal stabilization of the MT20, which exhibits a melting temperature 2.7\u0026deg;C higher than the WT.\u003c/p\u003e \u003cp\u003eTo investigate the conformational change of the WT and MT20 at elevated temperatures, we increase the simulation temperature to 100\u0026deg;C. This impractically high temperature is chosen to capture trends of significant conformational changes within the simulation time \u003csup\u003e\u003cspan additionalcitationids=\"CR23 CR24\" citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u003c/sup\u003e. We use the first 20 ns and the last 20 ns from a 200-ns production run at 100\u0026deg;C to create representative 3D structures. For the WT, we observe a significant conformational change at 100\u0026deg;C as compared to the 25\u0026deg;C structure. At room temp, the WT has its hydrophobic residues enclosing the active site facing inward, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec. When the temperature increases, the flexibility of two loops (131\u0026ndash;137 and 152\u0026ndash;156) that enclose the active pocket increases. As a result, highly hydrophobic residues M134, I135, V136, M137 and Y139 are exposed to the surface of the protein, increasing the aggregation propensity, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed. Moreover, the loop between β7 and αE, which contains D133, folds into the extension of αE that now includes residues 134\u0026ndash;141, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ee. The extended αE turns D133 away from H156 and results in a catalytically inactive conformation. This newly formed helix persists in the structure until the temperature is increased to 180\u0026deg;C and the time of simulation is prolonged to 1 \u0026micro;s, suggesting it may represent a metastable intermediate state in the thermal denaturation pathway of the enzyme. This conformation may be a representative intermediate state that leads to the irreversible aggregation of the WT upon heating.\u003c/p\u003e \u003cp\u003eHowever, for MT20, an increase in temperature during the simulation only increases the movement of flexible regions. No significant structural changes around the active pocket are identified, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ef, and no additional hydrophobic residues are exposed as compared to the native structure, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eg and \u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eh. The lack of αE extension at 134\u0026ndash;137 suggests that the mutations near this loop, including M134P, M137A and Y139S, disturb its formation. As proline has the lowest α-helix propensity of all amino acids\u003csup\u003e\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e\u003c/sup\u003e, the M134P substitution potentially perturbs the α-helix extension of residues 134\u0026ndash;137. Further stabilization of this loop is achieved by substituting large nonpolar residues M137 and Y139 with small or polar amino acids A and S. By rejecting the intermediate state in the WT where a large hydrophobic patch is exposed, MT20 avoids aggregation during thermal unfolding.\u003c/p\u003e \u003cp\u003eWe therefore conclude that the heat tolerance of BsLipA in MT20 is achieved through two molecular mechanisms: the formation of new hydrogen bonds and the absence of aggregation-prone intermediates during thermal unfolding. The latter can be considered a predominant stabilization factor for activity retention after exposure to high temperatures\u003csup\u003e\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e\u003c/sup\u003e. While we have not run an MD simulation for MT5, we postulate that MT5 adopts aggregation resistance from its four substitutions around the loop between β7 and αE (M134E, I135E, M137D, and Y139D). These substitutions replace nonpolar residues with negatively charged ones, likely generating electronic repulsion to prevent the formation of αE extension and the aggregation of denatured conformers\u003csup\u003e\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e,\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eIn addition to MT5 and MT20, six other variants of BsLipA have been reported to exhibit aggregation resistance upon thermal unfolding\u003csup\u003e\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e,\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e,\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e\u003c/sup\u003e. Each of these eight variants carries at least one methionine substitution at positions 134 or 137. Variants 5A and 5B\u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e, each carries 10 substitutions, but only the M134E or M137P substitution results in the resistance to thermal aggregation\u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e. Such experimental evidence, along with our MD simulation demonstrates the critical role of residues 134 and 137 in conferring aggregation resistance upon thermal unfolding. Notably, among the substitutions at residues 134 and 137 that occur in these eight variants, M134D, M134E, M134P, M137P and M137A, none could be achieved by epPCR with single NT changes, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea. In fact, when starting from a wild-type methionine, the categories of negatively charged residues, aromatic residues, and special cases G, P, C are all barely accessible by epPCR, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ec. The discovery of each one of these eight variants relied on saturation mutagenesis (Extended Data Table\u0026nbsp;1). In particular, the parent of 5A and 5B, 4D3, was discovered after four rounds of epPCR, four rounds of screening for thermotolerance, and three rounds of recombination. However, 4D3 only achieved a T\u003csub\u003e50\u003c/sub\u003e improvement of 15\u0026deg;C, and did not acquire resistance to thermal aggregation when heated to around its melting temperature. Only through site-saturation mutagenesis at residues 134 and 137 of 4D3 were the aggregation-resistant variants 5A and 5B identified. The mutagenesis trajectory leading to these variants underscores the importance of saturating the chemical diversity at key residues in mutagenesis libraries. By using RSM to prepare the library, identification of key residue positions and the saturation of their chemical diversity can be achieved simultaneously.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eMutagenesis libraries are the starting point of directed evolution campaigns and thereby determine the range of the improvement. With no prior information about the wild-type, a random mutagenesis method that uniformly samples all possible amino acid substitutions across all residue positions is critical to preserve the accessibility to beneficial properties. For an enzyme with N residues, the number of all possible single amino acid substitutions, 19xN, is typically a few thousand. Therefore, the number of unique variants with n substitutions is ~(19xN)\u003csup\u003en\u003c/sup\u003e. With only three random substitutions per gene, the theoretical sequence space exceeds the screening throughput of most current screening techniques\u003csup\u003e\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e\u003c/sup\u003e. For a given screen to effectively cover a meaningful subspace of an n-mutation sequence space, it is therefore important that the mutagenesis method allows mutation rate control to adjust the expected hit rate to match the screening throughput. All three critical features of a random mutagenesis method: mutating across the gene, saturating the diversity of each residue, and controlling mutation rates, are simultaneously provided by RSM.\u003c/p\u003e \u003cp\u003eWhen the screening throughput is limited to a few thousand, a single or low mutation rate RSM library can be prepared to enable mutational analysis of single amino acid substitutions and to guide recombination. If screening with even higher throughput is available, a high mutation rate RSM library can be implemented, directly enabling the discovery of significantly improved variants combining multiple beneficial mutations. The seven-mutation library used in our study has a theoretical diversity of (181 choose 7) x 20\u003csup\u003e7\u003c/sup\u003e = 1.4 x 10\u003csup\u003e21\u003c/sup\u003e. With a throughput of 10\u003csup\u003e7\u003c/sup\u003e, our screen samples only one in 10\u003csup\u003e15\u003c/sup\u003e in this sequence space. The discovery of significantly improved variants with such low sample rate suggests that similar high-performing variants are likely prevalent in the sequence space; therefore, it is possible to capture them with substantial down sampling. However, the sampling scheme, determined by the library preparation method, must retain the accessibility of each residue to all alternative chemical properties. If the library preparation method inherently causes significant loss of residue properties, such as with epPCR, the crucial substitutions may be eliminated from the accessible sequence space once the library is constructed. By enforcing uniform sampling, RSM reduces this bias and preserves a better representation of the total chemical space. The combination of a high mutation rate RSM library with high-throughput droplet screening enables us to more uniformly sample the vast sequence space and identify two variants with T\u003csub\u003e50\u003c/sub\u003e values improved by \u0026gt;\u0026thinsp;40\u0026deg;C compared to the wild-type enzyme in a single mutagenesis round. By contrast, using other methods, similar variants can only be achieved after 4\u0026ndash;5 rounds of mutagenesis, recombination and screening (Extended Data Table\u0026nbsp;1). When some prior information is available, RSM can be adjusted to reflect this knowledge and better facilitate the search of improved variants. For example, if previous studies have identified certain regions in the wild-type enzyme in conferring certain properties, oligo pools that cover those regions can be used in one mutagenesis reaction while the remaining oligo pools can be used in other mutagenesis reactions. When mixing mutagenesis products for StEP, the product of the first reaction can be used with a much higher molar ratio than the product of the second reaction, resulting in higher mutation rates in the previously-identified regions. If prior information further identifies beneficial properties at certain positions, the set of primers to mutate these positions can be designed to target substituted amino acids with those properties. These adjustments of RSM can focus the mutations to specific regions and specific substitutions, further improving the likelihood of discovering significantly improved variants. These results demonstrate the potential of RSM to accelerate directed evolution.\u003c/p\u003e"},{"header":"Material and Methods","content":"\u003cp\u003e \u003cb\u003eMicrofluidics device fabrication.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eWe design the devices (droplet generator and droplet sorter; see Fig. S12) using AutoCAD and producing them as photomasks (CAD/Art Services, Inc.). We fabricate the devices through the well-established techniques of soft lithography, employing SU8-on-silicon-wafer masters and PDMS-on-glass devices. We pour polydimethylsiloxane (PDMS) (Sylgard 184) onto the masters and bake them at 65\u0026deg;C overnight to cure the PDMS. Subsequently, we cautiously peel each PDMS device from the master and seal it to a pristine glass slide (Corning, 2947).\u003c/p\u003e \u003cp\u003eFor devices incorporating electrodes (droplet sorter), we integrate the electrodes into the design as channels within the microfluidic devices. We fill these channels with a low melting point metal alloy (Indalloy 19, 0.020 in. diameter) while heating the devices on a hot plate. We add terminal blocks to the punched holes in the devices to facilitate electrical connections during experiments. Following the fabrication process, we inject Aquapel (fluoroalkylsilanes) through the punched holes of the devices and use pressurized air to impel it through the microfluidic channel walls, rendering the internal surfaces of the channels hydrophobic. We expel excess Aquapel using compressed air and bake the device at 65\u0026deg;C overnight.\u003c/p\u003e\n\u003ch3\u003eLibrary Construction and Sequencing\u003c/h3\u003e\n\u003cp\u003eFor the single library, we prepare 18 parallel multisite-directed mutagenesis reactions in PCR strips, each containing one primer pool as the primers. After thermocycling, we purify the DNA product and perform overlapping PCR with wild-type BsLipA on the pET28a\u0026thinsp;+\u0026thinsp;vector backbone as the template. We then incubate the PCR product with DpnI to cleave methylated wild-type plasmid templates and column purify it using the Monarch DNA purification kit (New England Biolabs). We transform the plasmid library into competent cells (Turbo chemically competent cells, New England Biolabs), plate a portion to verify successful transformation, and transfer the remainder into liquid culture for plasmid extraction using the Monarch Plasmid extraction kit (New England Biolabs).\u003c/p\u003e \u003cp\u003eFor the 4\u0026times;3\u0026thinsp;+\u0026thinsp;3\u0026times;2 library, we prepare five parallel multisite-directed mutagenesis reactions, each containing three or four primer groups. Following thermocycling, we perform a StEP reaction on the purified DNA products to induce recombination of the resultant mutants. We prepare the reaction mix for StEP, consisting of a DNA shuffling mixture with 2.5 \u0026micro;L of template (1 ng/\u0026micro;L), 50 \u0026micro;L of 2\u0026times; Taq Master Mix (New England Biolabs), 1.5 \u0026micro;L of 10 \u0026micro;M forward primer (TAATACGACTCACTATAGGG), 1.5 \u0026micro;L of 10 \u0026micro;M reverse primer (GCTAGTTATTGCTCAGCGG), and MilliQ water to a total volume of 100 \u0026micro;L. We use a thermocycling program of 94\u0026deg;C for 30 seconds and 55\u0026deg;C for 5 seconds for 99 cycles. After thermocycling, we use a portion of the DNA product for gel electrophoresis to confirm its size and purify the remaining products using a DNA cleanup kit. The purified PCR products serve as megaprimers in an overlap PCR reaction with 250 ng of PCR products, 300 ng of wild-type BsLipA plasmids, 25 \u0026micro;L of 2\u0026times; Q5 Master Mix, and MilliQ water to a total volume of 50 \u0026micro;L. We use a thermocycling program of 98\u0026deg;C for 2 minutes, followed by 98\u0026deg;C for 10 seconds, 61\u0026deg;C for 20 seconds, 72\u0026deg;C for 3 minutes, and a final extension at 72\u0026deg;C for 5 minutes.\u003c/p\u003e \u003cp\u003eAfter the StEP and overlapping PCR, we add 1 \u0026micro;L of DpnI and incubate at 37\u0026deg;C for 3 hours. We purify the PCR product using a PCR cleanup kit, elute it into 10 \u0026micro;L, and transform it into Agilent electrocompetent cells. After adding 900 \u0026micro;L of SOC and incubating at 37\u0026deg;C and 250 rpm for 1 hour, we plate 20 \u0026micro;L on LB-trimethoprim plates and add 3 mL of LB-trimethoprim to a liquid culture for overnight growth at 37\u0026deg;C. The isolated plasmids are ready for further analysis.\u003c/p\u003e \u003cp\u003eSimilarly, for the 6\u0026times;3 library, we pool six primer groups at a time, creating a maximum of six mutations per round, followed by StEP. We purify the DNA product and perform overlap PCR with the wild-type plasmid as the template, which we later remove by DpnI digestion after the reaction. We then column purify the resultant DNA product, transform it into competent cells, plate a portion to verify transformation, and transfer the rest into liquid culture to grow overnight for plasmid extraction. For the 9\u0026times;2 library, we add odd-numbered primer pools first, then even-numbered pools, followed by StEP, overlap PCR, DpnI digestion, DNA column purification, transformation, and overnight growth in liquid culture for plasmid extraction. For the 18\u0026times;1 library, we add all primer pools together at once, creating a maximum of 18 mutations per round, followed by StEP, overlap PCR, DpnI digestion, DNA column purification, transformation, and overnight growth in liquid culture for plasmid extraction. For thermotolerance variant screening, we prepare a library with a higher mutation rate by using two consecutive 9\u0026times;2 RSM protocols.\u003c/p\u003e \u003cp\u003eWe also construct three error-prone libraries using the GeneMorph II kit (Agilent) with high, medium, and low mutation rates, resulting in 2.70, 2.04, and 1.18 mutations per gene, respectively, following the manufacturer\u0026rsquo;s manual. We transform the plasmid library into competent cells, plate a portion to verify successful transformation, transfer the remainder into liquid culture, and extract plasmids.\u003c/p\u003e \u003cp\u003eWe send all the libraries (single, 4\u0026times;3\u0026thinsp;+\u0026thinsp;3\u0026times;2, 6\u0026times;3, 9\u0026times;2, 18\u0026times;1, EPhigh, EPmid, EPlow) for sequencing using the PacBio platform for high-fidelity sequencing of long reads.\u003c/p\u003e\n\u003ch3\u003eSequencing Result Analysis\u003c/h3\u003e\n\u003cp\u003eWe filter sequencing reads based on the quality score of the lowest quality base in the open reading frame with a threshold of 90, eliminating any imperfect sequencing result. After filtering, we acquire approximately 10\u003csup\u003e5\u003c/sup\u003e reads for each library. We align the reads to the wild-type open reading frame, translate the aligned sequences into polypeptide sequence and compare the translation to wild-type translation. All analysis is conducted using a customized python package based on open source libraries such as Biopython (github.com/angiechen918/Mutan.git)\u003c/p\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eStatistical Analysis on Library quality\u003c/h2\u003e \u003cp\u003eTo analyze the percentage of reads with insertions and deletions, we filter all reads with aligned length of wild-type length +/- 3 and count the fraction of reads with either insertions or deletions. To analyze the percentage of reads with in-frame stop codons, we filter all reads with the correct length and count the fraction of reads that translates into a polypeptide with shorter length than the wild-type enzyme. To analyze the percentage of silenced mutations, we filter all reads with the correct length and identify the codon positions that contain nucleotide changes. For all codon positions that contain nucleotide changes, we count the fraction of codon positions that does not lead to a different amino acid than the wild-type enzyme at this position.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eStatistical Analysis on Distinct Amino Acid Substitution Abundance in Different Libraries\u003c/h2\u003e \u003cp\u003eWe filter all reads with the correct length and compare the translation to the wild-type enzyme to identify single amino acid substitutions. For each library, we analyze\u0026thinsp;~\u0026thinsp;72k single amino acid substitutions and categorize them based on the amino acid that they \u0026ldquo;mutate from\u0026rdquo; in the wild-type (19 amino acids, 20 - cysteine) and the amino acid that they \u0026ldquo;mutate to\u0026rdquo; in the variant (19 amino acids, 20 - the \u0026ldquo;mutate from\u0026rdquo; amino acid). We each of these 19x19\u0026thinsp;=\u0026thinsp;361 categories a chemically distinct substitution. As each amino acid may be used by the wild-type enzyme multiple times, we compile a table where the rows list the 19 amino acids in the wild-type BsLipA and the columns list the 19 possible mutations, resulting in a heatmap of the 361 chemically distinct amino acid substitutions. The color gradient reflects the abundance of the mutations, with warmer colors representing higher mutation frequencies. The numbers in each square of the heatmap represent the normalized count of each amino acid substitution. A number of 0 indicates that the substitution is enabled but occurs so infrequently that it rounds to 0 after normalization. Amino acid substitutions that can be enabled by a single nucleotide change in a trinucleotide codon are highlighted by white-bordered squares. There are 124 such squares, indicating that only 34.4% of possible substitutions are enabled by single nucleotide changes.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003eStatistical Analysis on Impact of Staggered Extension Process (StEP) on Mutation Distribution and Correlation\u003c/h2\u003e \u003cp\u003eTo decrease the positional correlation and increase the diversity of the final RSM library, we mix the products from the two reactions in equimolar ratios as templates for StEP. Before StEP, the number of mutations per region is limited to one, and the mutations are fully confined within the regions covered by the primer pools of choice (Sup.Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). In the StEP thermocycling program, we limit the extension step to a few seconds, preventing the DNA polymerase from synthesizing the full gene in each cycle. During the denaturing step of the next cycle, incomplete strands unbind from the original template and anneal to different template strands, effectively combining mutations from various templates. Therefore, after StEP, the positions for mutations are no longer fully correlated.\u003c/p\u003e \u003cp\u003eWe reflect this change by counting the number of reads where a specific pair of mutation positions appear on the same read, N(A,B), choosing A\u0026thinsp;\u0026gt;\u0026thinsp;B. For example, if a variant in the library contains three mutations at positions 1, 21, and 31, then this read contributes counts to N(1,21), N(1,31), and N(21,31). When mutations at position A coincide frequently with mutations at position B, we expect a high N(A,B) value. By contrast, if the mutations at positions A and B never coincide on the same read, N(A,B) should be 0. The abundance of individual positions A and B also strongly influences N(A,B). To calculate the association between positions without considering the abundance of individual positions, we divide the number of position pairs by the occurrence of either position in the pair and define this value as Mutual Association, M(A,B)\u0026thinsp;=\u0026thinsp;N(A)\u0026times;N(B)N(A,B)​. This value is similar to the exponential of pointwise mutual information (PMI), a concept broadly used in Natural Language Processing (NLP) to evaluate the association of two words in a corpus (Sup.Fig. S2).\u003c/p\u003e \u003cp\u003eIn the 9\u0026times;2 library before StEP, M(A,B) values are high if positions A and B are in pools with indices of the same parity, resulting from the grouping of pools in the mutagenesis reaction. By contrast, if A and B are in pools with indices of opposite parity, M(A,B) values are close to 0, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ee. However, this pattern in the association of positions becomes relaxed after StEP. In the 9\u0026times;2 library after StEP, as long as the distance between positions A and B is longer than ~\u0026thinsp;30 nucleotides, M(A,B) values are uniform, indicating low bias in position-wise association, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ef. For position pairs with short distances, positions in pools with the same indices of opposite parity are still more associated than positions in pools with opposite indices of opposite parity. We expect such bias to be further eliminated with shorter annealing times during StEP or more rounds of StEP. Notably, despite the dramatic shift in mutation distribution with StEP, the average mutation rate of the 9\u0026times;2 library remains consistent before (4.87\u0026thinsp;\u0026plusmn;\u0026thinsp;2.57) and after StEP (4.64\u0026thinsp;\u0026plusmn;\u0026thinsp;2.40), ensuring predictability when creating a mutagenesis library with a target mutation rate.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eDroplet Sorting\u003c/h2\u003e \u003cp\u003eWe co-encapsulate each bacterial host cell expressing the mutant protein with resorufin acetate substrate and a bacterial lysis reagent (BugBuster Protein Extraction Reagent, Sigma-Aldrich) using a dual-inlet droplet maker (Fig. S7). We collect the emulsion under a layer of mineral oil into an Eppendorf tube sitting in a heat block set to 65\u0026deg;C and incubate it for 1 hour to inactivate the non-thermostable variants.\u003c/p\u003e \u003cp\u003eNext, we carefully remove all the mineral oil with a pipette and re-inject the emulsions into a microfluidic sorting device (Fig. S7). With an additional stream of spacing oil, the closely packed droplets are evenly spaced and enter the sorter periodically. We place a 532 nm excitation laser at the entrance of the sorter. If the fluorescent reaction product in a droplet exceeds the set fluorescence threshold, we use an electric field to pull it into an adjacent channel and collect it through polyethylene micro-tubing into an Eppendorf tube placed on ice. We operate the chip at 600\u0026ndash;1000 drops per second, probing approximately 3\u0026nbsp;million cells per hour. By applying electric fields of 50 V at a frequency of 25 kHz, we screen more than 3\u0026nbsp;million droplets and sort 1,260 droplets. We break the emulsion by adding 200 \u0026micro;L of 20% (v/v) PFO (1H,1H,2H,2H-perfluoro-1-octanol, Sigma-Aldrich) in HFE 7500 (3M) and purify the aqueous phase using a DNA cleanup kit (Monarch\u0026reg; PCR \u0026amp; DNA Cleanup Kit, New England Biolabs), eluting the DNA in 5 \u0026micro;L of ddH₂O.\u003c/p\u003e \u003cp\u003eWe then transform the plasmids into high-efficiency competent cells (Turbo Competent \u003cem\u003eE. coli\u003c/em\u003e, New England Biolabs) for plasmid replication and extraction.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003eResidual Activity Measurement in Cell Lysate\u003c/h2\u003e \u003cp\u003eWe transform the plasmids from the sorted library into BL21(DE3) competent cells (NEB C2527H) and plate 200 \u0026micro;L of the culture onto agarose LB-kanamycin plates, incubating them at 37\u0026deg;C overnight for colony growth. We pick 94 colonies and transfer them to a 96-well plate, including one TBS plus substrate blank and a wild-type control. We shake the plate at 37\u0026deg;C for 12 hours to allow the colonies to grow overnight in 200 \u0026micro;L of LB-kanamycin liquid culture. We then dilute the culture 20-fold, with triple replicates, and add BugBuster reagent to each well. We incubate the plate at 37\u0026deg;C for one hour, then increase the temperature to 70\u0026deg;C and incubate for an additional hour to inactivate non-thermostable variants. Subsequently, we add 7.5 \u0026micro;M resorufin acetate while keeping the plate on ice throughout the procedure. We measure the relative fluorescence units (RFU) at 549 nm excitation and 585 nm emission over time. Two variants are identified with residual activity significantly higher than the wild type after incubation at 70\u0026deg;C for one hour.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec16\" class=\"Section2\"\u003e \u003ch2\u003eProtein expression and purification\u003c/h2\u003e \u003cp\u003eWe grow the identified variants in liquid culture and induce them using Invitrogen's MagicMedia \u003cem\u003eE. coli\u003c/em\u003e Expression Medium at 30\u0026deg;C overnight. We purify the protein using Ni-NTA spin columns (Thermo Scientific, 88224). To remove imidazole, we dialyze the eluted sample against 1X DPBS at 4\u0026deg;C using for 4 hours a 7k MWCO dialysis cassette (Thermo Scientific 66373). After dialysis, we use SDS-PAGE to confirm the molecular weight of the protein at the correct size. To preserve the purified enzymes, we standardize their concentrations to 0.4 mg/mL in 50% glycerol in 1x DPBS and store them at -20\u0026deg;C before use. We then characterize these purified proteins for their T₅₀, residual activity, and melting temperature.\u003c/p\u003e \u003cp\u003eT\u003csub\u003e50\u003c/sub\u003e \u003cb\u003eMeasurement\u003c/b\u003e\u003c/p\u003e \u003cp\u003eTo measure the T\u003csub\u003e50\u003c/sub\u003e of enzyme mutants, we dilute the enzyme stocks 100-fold by mixing 10 \u0026micro;L of the enzyme with 1,490 \u0026micro;L of 1\u0026times; TBS buffer. We then incubate the diluted enzyme solution on three PCR machines simultaneously at various temperatures (gradient 35\u0026ndash;60\u0026deg;C, gradient 65\u0026ndash;90\u0026deg;C, and 95\u0026deg;C) for 20 minutes, with a total volume of 120 \u0026micro;L in each PCR tube, selecting seven temperature points: 36, 45, 55, 66, 75, 85, 95\u0026deg;C, and leaving one sample at room temperature (25\u0026deg;C). Following incubation, we further dilute the enzyme mutants 10-fold by adding 12 \u0026micro;L of the enzyme solution to 108 \u0026micro;L of TBS buffer. We then add this diluted enzyme to each well of a 96-well plate.\u003c/p\u003e \u003cp\u003eWe prepare the substrate by diluting 20 \u0026micro;L of 10 mM resorufin acetate (RA) in 10 mL of 1\u0026times; TBS buffer, freshly prepared for each group. In each well of the 96-well plate, we add 12 \u0026micro;L of the 1,000-fold diluted enzyme and 12 \u0026micro;L of the substrate, ensuring thorough mixing. We initiate the measurement program at 25\u0026deg;C on a qPCR machine (Bio-Rad CFX96 Real-Time PCR System), with each point measured for 15 seconds using the Texas Red channel. The groups tested include mutants 5 and 20, wild type, and a glycerol-only control.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec17\" class=\"Section2\"\u003e \u003ch2\u003eResidual Activity at High Temperature Measurement with Purified Protein\u003c/h2\u003e \u003cp\u003eTo assess both the residual activity and the high-temperature enzyme activity of the identified \u003cem\u003eBacillus subtilis\u003c/em\u003e Lipase A (BsLipA) variants, we use the Bio-Rad CFX96 real-time PCR system to measure fluorescent products at various temperatures. The groups tested include mutants 5 and 20, wild type, and a glycerol-only control. First, we dilute purified enzyme stocks (0.4 mg/mL) in 1\u0026times; TBS buffer (pH 7.4) at a ratio of 1:500 (2 \u0026micro;L in 1,000 \u0026micro;L), adding an additional 10 \u0026micro;L of 50% glycerol to match the glycerol concentration. For the WT enzyme, which has an initial concentration of 0.5 mg/mL, we first dilute it at a ratio of 4:5, then follow the same dilution protocol as for the mutants.\u003c/p\u003e \u003cp\u003eNext, we subject the enzyme solutions to heat inactivation at 80\u0026deg;C for 1 hour using Eppendorf tubes on a heat block. Immediately after heating, we place the tubes on ice to stop the heat inactivation process. We prepare the substrate by diluting 20 \u0026micro;L of 25 mM fluorescein diacetate (FDA) in 10 mL of 1\u0026times; TBS buffer, freshly prepared for each group. In each well of a 96-well plate, we add 12 \u0026micro;L of the 1,000-fold diluted enzyme and 12 \u0026micro;L of the substrate, ensuring thorough mixing. We initiate the measurement program at 41, 45.6, 50.2, 56, 60.8, and 65\u0026deg;C on a qPCR machine (Bio-Rad CFX96 Real-Time PCR System), with each point measured for 15 seconds using the SYBR channel.\u003c/p\u003e \u003cp\u003e \u003cb\u003eMelting Temperature (\u003c/b\u003eT\u003csub\u003em\u003c/sub\u003e\u003cb\u003e) Measurement\u003c/b\u003e\u003c/p\u003e \u003cp\u003eWe use NanoDSF to measure the melting points of the purified proteins MT5, MT20, and WT. We load prepared samples into a 384-well plate and then into NanoTemper NT.Plex Capillary Chips. We perform the experiments in duplicates. We heat and cool the samples from 20\u0026deg;C to 95\u0026deg;C to 20\u0026deg;C at a rate of 1.5\u0026deg;C/min and collect the data using the NanoTemper PR.ThermControl software.\u003c/p\u003e \u003cp\u003eNanoDSF is an advanced differential scanning fluorimetry technique that monitors changes in intrinsic fluorescence of tryptophan and tyrosine residues as proteins unfold. As proteins are exposed to a temperature gradient, their local chemical environment changes, causing variations in fluorescence intensity and wavelength, particularly when aromatic side chains become solvent-exposed upon unfolding. NanoDSF measures fluorescence at 330 nm and 350 nm and calculates the ratio as a function of temperature to determine the melting temperature (T\u003csub\u003em\u003c/sub\u003e), which is the temperature at which 50% of the protein population is unfolded. NanoDSF offers higher throughput and lower sample consumption compared to traditional methods like differential scanning calorimetry (DSC) or circular dichroism (CD), with comparable accuracy for determining folding free energies and unfolding temperatures. NanoDSF requires the presence of aromatic amino acids in the protein of interest.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec18\" class=\"Section2\"\u003e \u003ch2\u003eMeasurement of Soluble Protein Fraction after Heat Inactivation\u003c/h2\u003e \u003cp\u003eTo assess protein solubility upon heat incubation, we dilute purified enzyme stocks MT5 and MT20, along with wild-type (WT) protein (0.4 mg/ml in 50% glycerol 1xDPBS), to 0.025 mg/ml in 1x DPBS buffer. We subject the diluted enzyme solutions to heat inactivation at 65\u0026deg;C for 10 minutes in Eppendorf tubes placed on a heat block. After incubation, we centrifuge the heated samples along with non-heated controls at 12,000 \u0026times; \u003cem\u003eg\u003c/em\u003e for 10 minutes to separate soluble proteins from precipitated aggregates. We carefully collect the supernatants to determine the concentration of soluble protein using Micro BCA (Thermo Scientific, 23235) following the manufacturer instruction. The absorption at 562 nm is measured using a NanoDrop spectrophotometer.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec19\" class=\"Section2\"\u003e \u003ch2\u003eMolecular Dynamics Simulations of Protein Folding at Low and High Temperatures\u003c/h2\u003e \u003cp\u003eFor molecular dynamics simulations we use crystal structure of BsLipA PDB ID: 5CRI. For molecular dynamics simulations, protein preparation and analysis we use Amber22\u003csup\u003e29\u003c/sup\u003e software with GPU support and AmberTools23\u003csup\u003e30\u003c/sup\u003e suite of programs. We use ff19SB protein force field\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e and OPC 4-point water model\u003csup\u003e\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u003c/sup\u003e. We initially minimize the system for 2000 steps using Newton-Raphson algorithm, and then slowly heat it from 100 K to final simulation temperature for 1 ns with time-step of 1 fs, while keeping the volume constant. The initial restraint on protein atoms is 100 kcal/mol-\u0026Aring;\u003csup\u003e2\u003c/sup\u003e. We then relax the system for total of 6 ns with 1 fs time step, in NPT ensemble, by slowly lifting restraints (to 10.0, 1.0 and 0.1 kcal/mol-\u0026Aring;\u003csup\u003e2\u003c/sup\u003e) on protein backbone atoms. Finally, we perform 20 ns of equilibration in NPT ensemble with 1 fs time step and no restraints. We perform 200 ns production runs using NPT ensemble with Langevin thermostat and Berendsen barostat using 2 fs time step. To impose constraints on hydrogen atoms we employ SHAKE algorithm. To visualize MD trajectories, we use VMD 1.9.3 software. We introduce mutations into wild-type protein using SCWRL4.0\u003csup\u003e33\u003c/sup\u003e software. We perform molecular dynamics simulations at 25\u0026deg;C and increased temperatures for both M20 and WT to 100\u0026deg;C. We perform all simulations in triplicates. We use 200 ns trajectory information to create representative 3D structures for MT20 and WT using SINAPs extension in UCSF Chimera 1.17.3\u003csup\u003e34\u003c/sup\u003e. We align these structures and visualize interactions present only in the mutated protein. Additionally, we use the first 20 ns after equilibration and the last 20 ns from simulations at 100\u0026deg;C to create representative 3D structures for MT20 and WT and overlap them to visualize structural differences. We also use UCSF Chimera to visualize hydrophobic surfaces in representative 3D structures.\u003c/p\u003e \u003c/div\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAuthor Contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eX.D.Z., A.C., K.M. and D.A.W. developed the concept. X.D.Z, A.C., R.H., P.E. and K.M. designed the experiments and analysis approaches. X.D.Z, A.C., LC, Y.X., R.H., P.E. and K.M performed the experiments and conducted the analyses. MM performed and analyzed computational simulations. X.D.Z, A.C., D.A.W. and K.M. wrote the manuscript with input from all authors.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgments\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis work was supported in part by the National Science Foundation (CBET-2103538), the Harvard MRSEC (DMR-2011754), the Wyss Institute and Science Fund of the Republic of Serbia, #7750288.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConflict of Interest\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eDavid A. Weitz and Karla Milcic are co-founders of RevivBio Inc.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData Availability Statement\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe data that support the findings of this study are available in the supplementary material and extended data of this article and upon request from corresponding authors.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCode Availability Statement\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll analysis of PacBio sequencing data is conducted using a customized python package based on open-source libraries such as Biopython and is available at: github.com/angiechen918/Mutan.git\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003ePacker, M. S. \u0026amp; Liu, D. R. Methods for the directed evolution of proteins. \u003cem\u003eNature Reviews Genetics\u003c/em\u003e \u003cstrong\u003e16\u003c/strong\u003e, 379-394 (2015). https://doi.org:10.1038/nrg3927\u003c/li\u003e\n\u003cli\u003eEsvelt, K. M., Carlson, J. C. \u0026amp; Liu, D. R. A system for the continuous directed evolution of biomolecules. \u003cstrong\u003e472\u003c/strong\u003e, 503 (2011). \u003c/li\u003e\n\u003cli\u003eRosenthal, R. G., Zhang, X. D., Đurđić, K. I., Collins, J. J. \u0026amp; Weitz, D. A. Controlled Continuous Evolution of Enzymatic Activity Screened at Ultrahigh Throughput Using Drop-Based Microfluidics. \u003cstrong\u003e62\u003c/strong\u003e (2023). \u003c/li\u003e\n\u003cli\u003eDrummond, D. A., Iverson, B. L., Georgiou, G. \u0026amp; Arnold, F. H. Why high-error-rate random mutagenesis libraries are enriched in functional and improved proteins. \u003cstrong\u003e350\u003c/strong\u003e, 816 (2005). \u003c/li\u003e\n\u003cli\u003eWong, T. S., Roccatano, D., Zacharias, M. \u0026amp; Schwaneberg, U. A Statistical Analysis of Random Mutagenesis Methods Used for Directed Protein Evolution. \u003cem\u003eJournal of Molecular Biology\u003c/em\u003e \u003cstrong\u003e355\u003c/strong\u003e, 858-871 (2006). https://doi.org:10.1016/j.jmb.2005.10.082\u003c/li\u003e\n\u003cli\u003eFrauenkron-Machedjou, V. J.\u003cem\u003e et al.\u003c/em\u003e Towards Understanding Directed Evolution: More than Half of All Amino Acid Positions Contribute to Ionic Liquid Resistance of Bacillus subtilis Lipase A. \u003cstrong\u003e16\u003c/strong\u003e, 945 (2015). \u003c/li\u003e\n\u003cli\u003eFrauenkron-Machedjou, V. J.\u003cem\u003e et al.\u003c/em\u003e Exploring the full natural diversity of single amino acid exchange reveals that 40\u0026ndash;60% of BSLA positions improve organic solvents resistance. \u003cstrong\u003e5\u003c/strong\u003e (2018). \u003c/li\u003e\n\u003cli\u003eMarkel, U.\u003cem\u003e et al.\u003c/em\u003e Are directed evolution approaches efficient in exploring nature\u0026rsquo;s potential to stabilize a lipase in organic cosolvents? \u003cstrong\u003e7\u003c/strong\u003e (2017). \u003c/li\u003e\n\u003cli\u003eNutschel, C., Fulton, A., Zimmermann, O., Schwaneberg, U., Jaeger, K. E. \u0026amp; Gohlke, H. Systematically Scrutinizing the Impact of Substitution Sites on Thermostability and Detergent Tolerance for Bacillus subtilis Lipase A. \u003cstrong\u003e60\u003c/strong\u003e, 1568-1584 (2020). \u003c/li\u003e\n\u003cli\u003eCheng, F., Zhu, L. \u0026amp; Schwaneberg, U. Directed evolution 2.0: Improving and deciphering enzyme properties. \u003cstrong\u003e51\u003c/strong\u003e, 9772 (2015). \u003c/li\u003e\n\u003cli\u003eWong, T. S., Tee, K. L., Hauer, B. \u0026amp; Schwaneberg, U. Sequence saturation mutagenesis (SeSaM): a novel method for directed evolution. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e32\u003c/strong\u003e, e26 (2004). https://doi.org:10.1093/nar/gnh028\u003c/li\u003e\n\u003cli\u003eRuff, A. J., Kardashliev, T., Dennig, A. \u0026amp; Schwaneberg, U. The Sequence Saturation Mutagenesis (SeSaM) method. \u003cem\u003eMethods Mol Biol\u003c/em\u003e \u003cstrong\u003e1179\u003c/strong\u003e, 45-68 (2014). https://doi.org:10.1007/978-1-4939-1053-3_4\u003c/li\u003e\n\u003cli\u003eWong, T. S.\u003cem\u003e et al.\u003c/em\u003e Transversion-enriched sequence saturation mutagenesis (SeSaM-Tv+): a random mutagenesis method with consecutive nucleotide exchanges that complements the bias of error-prone PCR. \u003cem\u003eBiotechnol J\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 74-82 (2008). https://doi.org:10.1002/biot.200700193\u003c/li\u003e\n\u003cli\u003eZhao, J., Kardashliev, T., Joelle Ruff, A., Bocola, M. \u0026amp; Schwaneberg, U. Lessons from diversity of directed evolution experiments by an analysis of 3,000 mutations. \u003cem\u003eBiotechnol Bioeng\u003c/em\u003e \u003cstrong\u003e111\u003c/strong\u003e, 2380-2389 (2014). https://doi.org:10.1002/bit.25302\u003c/li\u003e\n\u003cli\u003eZhao, J.\u003cem\u003e et al.\u003c/em\u003e Amino acid substitutions in random mutagenesis libraries: lessons from analyzing 3000 mutations. \u003cstrong\u003e101\u003c/strong\u003e, 3187 (2017). \u003c/li\u003e\n\u003cli\u003eReetz, M. T., Carballeira, J. D. \u0026amp; Vogel, A. Iterative saturation mutagenesis on the basis of b factors as a strategy for increasing protein thermostability. \u003cstrong\u003e45\u003c/strong\u003e, 7751 (2006). \u003c/li\u003e\n\u003cli\u003eAhmad, S. \u0026amp; Rao, N. M. Thermally denatured state determines refolding in lipase: Mutational analysis. \u003cstrong\u003e18\u003c/strong\u003e, 1196 (2009). \u003c/li\u003e\n\u003cli\u003eSell\u0026eacute;s Vidal, L., Isalan, M., Heap, J. T. \u0026amp; Ledesma-Amaro, R. A primer to directed evolution: current methodologies and future directions. \u003cstrong\u003e4\u003c/strong\u003e, 291 (2023). \u003c/li\u003e\n\u003cli\u003eZhao, H. \u0026amp; Zha, W. In vitro \u0026apos;sexual\u0026apos; evolution through the PCR-based staggered extension process (StEP). \u003cem\u003eNat Protoc\u003c/em\u003e \u003cstrong\u003e1\u003c/strong\u003e, 1865-1871 (2006). https://doi.org:10.1038/nprot.2006.309\u003c/li\u003e\n\u003cli\u003eAhmad, S., Kamal, M. Z., Sankaranarayanan, R. \u0026amp; Rao, N. M. Thermostable Bacillus subtilis Lipases: In Vitro Evolution and Structural Insight. \u003cstrong\u003e381\u003c/strong\u003e, 340 (2008). \u003c/li\u003e\n\u003cli\u003eMala, J. G. S. T., Satoru.\u0026thinsp;. Understanding Structural Features of Microbial LipasesaAn Overview. \u003cem\u003eAnalytical Chemistry Insights\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 9-19 (2008). \u003c/li\u003e\n\u003cli\u003eSingh, B., Bulusu, G. \u0026amp; Mitra, A. Understanding the thermostability and activity of bacillus subtilis lipase mutants: Insights from molecular dynamics simulations. \u003cstrong\u003e119\u003c/strong\u003e, 409 (2015). \u003c/li\u003e\n\u003cli\u003eHuang, X. \u0026amp; Zhou, H. X. Similarity and difference in the unfolding of thermophilic and mesophilic cold shock proteins studied by molecular dynamics simulations. \u003cstrong\u003e91\u003c/strong\u003e, 2463 (2006). \u003c/li\u003e\n\u003cli\u003eMeharenna, Y. T. \u0026amp; Poulos, T. L. Using molecular dynamics to probe the structural basis for enhanced stability in thermal stable cytochromes P450. \u003cstrong\u003e49\u003c/strong\u003e, 6686 (2010). \u003c/li\u003e\n\u003cli\u003ePurmonen, M., Valjakka, J., Takkinen, K., Laitinen, T. \u0026amp; Rouvinen, J. Molecular dynamics studies on the thermostability of family 11 xylanases. \u003cstrong\u003e20\u003c/strong\u003e, 559 (2007). \u003c/li\u003e\n\u003cli\u003eScholtz, C. N. P. J. M. A Helix Propensity Scale Based on Experimental Studies of Peptides and Proteins. \u003cem\u003eBiophysical Journal\u003c/em\u003e \u003cstrong\u003e75\u003c/strong\u003e, 422-427 (1998). \u003c/li\u003e\n\u003cli\u003eKamal, M. Z.\u003cem\u003e et al.\u003c/em\u003e In vitro evolved non-aggregating and thermostable lipase: Structural and thermodynamic investigation. \u003cstrong\u003e413\u003c/strong\u003e, 741 (2011). \u003c/li\u003e\n\u003cli\u003eAutour, A. \u0026amp; Ryckelynck, M. Ultrahigh-throughput improvement and discovery of enzymes using droplet-based microfluidic screening. \u003cstrong\u003e8\u003c/strong\u003e, (2017). \u003c/li\u003e\n\u003cli\u003eCase, D. A. \u003cem\u003eet al. \u003c/em\u003eAmber 2022, University of California, San Francisco. (2022).\u003c/li\u003e\n\u003cli\u003eCase, D. A. \u003cem\u003eet al.\u003c/em\u003e AmberTools. \u003cem\u003eJ Chem Inf Model\u003c/em\u003e \u003cstrong\u003e63\u003c/strong\u003e, 6183-6191 (2023).\u003c/li\u003e\n\u003cli\u003eTian, C. \u003cem\u003eet al.\u003c/em\u003e ff19SB: Amino-Acid-Specific Protein Backbone Paramters Trained against Quantum Mechanics Energy Surfaces in Solution. \u003cem\u003eJ Chem Theory Comput \u003c/em\u003e\u003cstrong\u003e16\u003c/strong\u003e, 528-552.\u003c/li\u003e\n\u003cli\u003eIzadi, S. Anandakrishnan, R. \u0026amp; Onufriev, A. V. Building Water Models: A Different Approach. \u003cem\u003eJ Phys Chem Lett\u003c/em\u003e \u003cstrong\u003e5\u003c/strong\u003e, 3863\u0026ndash; 3871 (2014).\u003c/li\u003e\n\u003cli\u003eKrivov, G. G. Shapovalov, M. V. \u0026amp; Dunbrack, Jr. R. L. Improved prediction of protein side-chain conformations with SCWRL4. \u003cem\u003eProteins\u003c/em\u003e \u003cstrong\u003e77\u003c/strong\u003e, 778-795 (2009) \u003c/li\u003e\n\u003cli\u003eBedart C. Renault N. Chavatte P. Porcherie, A. Lachgar, A. Capron, M. \u0026amp; Farce A.\u003cstrong\u003e\u003cem\u003e \u003c/em\u003e\u003c/strong\u003eSINAPs: A Software Tool for Analysis and Visualization of Interaction Networks of Molecular Dynamics Simulations. \u003cem\u003eJ Chem Inf Model\u003c/em\u003e \u003cstrong\u003e62\u003c/strong\u003e, 1425-1436 (2022)\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"nature-portfolio","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"","title":"Nature Portfolio","twitterHandle":"","acdcEnabled":false,"dfaEnabled":false,"editorialSystem":"ejp","reportingPortfolio":"","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"Directed evolution, Mutagenesis library generation, Enzyme engineering, Ultra-high throughput screening, Synthetic biology","lastPublishedDoi":"10.21203/rs.3.rs-5583255/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-5583255/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eDirected evolution improves the function of enzymes by diversifying the wild-type to generate variants and searching for the desired functions. When prior knowledge about the wild-type is limited, random mutagenesis is employed to introduce mutations at arbitrary positions. However, conventional methods exhibit significant bias in the distribution of amino acid substitutions, resulting in inaccessible chemical properties, limiting the chemical diversity, and hindering discovery of substantially improved variants. To overcome such limitations, we develop Random Saturation Mutagenesis (RSM) that uniformly introduces all amino acid substitutions across the gene at controlled rates by using overlapping site-mutagenesis oligo pools. Sequencing results illustrate the high quality and diversity of RSM libraries as compared to conventional libraries. We demonstrate the utility of RSM by improving thermotolerance of LipaseA. With only one round of mutagenesis, we identify two variants with over 40\u0026deg;C improvement in thermotolerance, demonstrating the potential of RSM to transform and accelerate directed evolution.\u003c/p\u003e","manuscriptTitle":"Random Saturation Mutagenesis to Generate Highly Diverse Libraries for Directed Evolution","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-01-24 07:44:27","doi":"10.21203/rs.3.rs-5583255/v1","editorialEvents":[],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"nature-chemical-biology","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"nchembio","sideBox":"Learn more about [Nature Chemical Biology](http://www.nature.com/nchembio/)","snPcode":"","submissionUrl":"","title":"Nature Chemical Biology","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"ejp","reportingPortfolio":"Nature Research","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"b67fb3d1-ab6c-4967-8db6-7e32039e8138","owner":[],"postedDate":"January 24th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"under-review","subjectAreas":[{"id":43263065,"name":"Biological sciences/Biochemistry/Biocatalysis"},{"id":43263066,"name":"Biological sciences/Chemical biology/Enzymes"},{"id":43263067,"name":"Biological sciences/Biological techniques/High-throughput screening"}],"tags":[],"updatedAt":"2025-09-20T20:50:10+00:00","versionOfRecord":[],"versionCreatedAt":"2025-01-24 07:44:27","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-5583255","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-5583255","identity":"rs-5583255","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}