Machine learning-guided evolution of pyrrolysyl-tRNA synthetase for improved incorporation efficiency of diverse noncanonical amino acids

Design of combinatorial variants of PylRS using the FFT-PLSR model

Mutations in the TBD of MmPylRS have been demonstrated to improve the efficiency of ncAA incorporation. It was also shown that N-terminal mutations did not significantly influence the substrate specificity of the PylRS and could be transferred to different variants¹⁹. We here tested four sets of mutations obtained previously in the TBD of MmPylRS to improve catalytic efficiency of IFRS, which include R61K/H63Y/S193R, R19H/H29R/T122S, D2N/K3N/T56P/H62Y, and V31I/T56P/H62Y/A100E^6,7,8,20 (Supplementary Table 1). Based on sequence alignment, we introduced these N-terminal mutations into IFRS and tested their activity for incorporating 3-bromo-Phe (3BrF), a cheaper substrate than 3IF (Fig. 1a, Supplementary Fig. 1). Expression of IFRS was driven by the constitutive, mid-strength E. coli glutaminyl-tRNA synthetase (glnS) promoter, and expression of PylT was controlled by the E. coli lpp gene promoter. Amber suppression of the sfGFPS2TAG gene by 3BrF was investigated by measuring the fluorescence intensity, and the ncAA-containing protein yield was presented by the ratio of fluorescence intensity to optical density OD₆₀₀ (Flu/OD) of cells expressing sfGFP and PylRS. The normalized protein yield was calculated by subtracting the Flu/OD ratio of cells cultured in the presence of ncAA with that of in the absence of ncAA (Supplementary Fig. 2). We found that R19H/H29R/T122S did not achieve the expected increase in stop codon suppression (SCS) efficiency, and the IPYE (V31I/T56P/H62Y/A100E) from chPylRS did not increase the SCS efficiency of IFRS either. Although R61K/H63Y/S193R in IFRS did achieve a modest increase in sfGFPS2TAG expression yield compared to the IFRS alone, the effect was still lower than that observed for the mutations in wild-type MmPylRS. Interestingly, the D2N/K3N/T56P/H62Y from chPylRS increased the SCS efficiency of IFRS by around 7-fold, confirming that tRNA binding domain mutations can indeed enhance the activity of catalytic domain mutants (Supplementary Fig. 2b). We also tested the activity of 12 single-point mutations and found that only D2N, R61K, and H62Y enhanced the activity of IFRS, with D2N being the most active, exhibiting 3-fold improved SCS efficiency compared to the wild type (Fig. 1b, Supplementary Fig. 2a). In the combinatorial mutant D2N/K3N/T56P/H62Y, D2N and H62Y enhanced the SCS efficiency of IFRS, while K3N and T56P reduced the SCS efficiency, indicating a positive sign epistasis among the mutations (Supplementary Table 2).

We then attempted to use machine learning to explore the combinatorial space of these 12 single-point mutations, a total of 4096 (2¹²) variants. The FFT-partial least squares regression (PLSR) approach was applied for predicting the fitness of combinatorial variants, which uses FFT for the protein sequence encoding and PLSR as the algorithm of the ML model. The FFT-PLSR model has demonstrated the ability to be trained on a small dataset of enzyme mutant activity data to accurately predict the activity across the entire combinatorial space. The variant sequences were first transferred into numerical features and then transformed into a two-dimensional energy versus frequency representation using FFT (Supplementary Fig. 3). Based on the transformed data, a PLSR model was trained to predict the activities of other mutants. We first trained the FFT-PLSR model using dataset 1, composed of 12 single variants datasets. By scoring with leave-one-out cross-validation (LOOCV), we screened 566 amino acid encodings from the AAindex database and selected the index with the best score to encode the amino acid sequences and build the PLS regression model (Supplementary Figs. 4 and 5). We constructed 6 double and 19 triple mutants, and measured their SCS efficiency to form dataset 2 as a test set (N = 25) (Fig. 1c). The trained model achieved an R² of 0.843 and an MSE of 2.887 on the test set, indicating that the model showed good prediction ability for high-activity combinatorial mutants (Fig. 1d). Interestingly, the best combinatorial mutant predicted by the model was D2N/R61K/H62Y, which was also the variant with highest activity in the test set (Supplementary Table 3).

In dataset 2, we observed epistasis between mutations. For example, T56P showed a decreased amber codon suppression efficiency, while H62Y improved the amber codon suppression efficiency compared to the IFRS. However, the T56P/H62Y showed a normalized fluorescence intensity significantly higher than that of H62Y, which led to a positive sign epistasis between T56P and H62Y. The positive sign epistasis effect was also found for R61K and H63Y, and D2N and K3N (Supplementary Table 2). By contrast, R19H and H29R showed a positive reciprocal sign epistasis (Supplementary Table 2). To help the model gain a more thorough understanding of epistasis between mutants, we added dataset 2 to the training set, which raised the total number of mutants in the training set to 38. Additionally, we rationally designed and constructed an additional test dataset 3 (N = 56) including 56 combinatorial mutants, mainly based on combination of improved variants (Fig. 1e). The retrained model achieved an R² of 0.835 on the test set, still showing high accuracy for predicting the activities of combinatorial mutants (Fig. 1f). We then constructed the top 8 mutants predicted by the model and tested their SCS efficiency (Fig. 1g and Supplementary Table 4). All eight variants showed a significant increase in the activity compared to the IFRS, with the mutant Z7 (D2N/V31I/T56P/R61K/H62Y/T122S/S193R, Com1-IFRS) exhibiting the highest SCS efficiency, which was 11-fold higher than the IFRS (Fig. 1g). We then attempted to apply the ML model to predict single-point variants in the TBD to improve the activity of Com1-IFRS. Fifteen single-point variants were predicted by the model to be more active than Com1-IFRS, with 12 of them located at previously unseen positions. However, all 15 variants failed to further improve the activity of Com1-IFRS, indicating the model’s limited ability to predict effects of mutations at unseen positions (Supplementary Fig. 6). The outcome is reasonable, given that the training set was small, only containing 38 variants across the 12 positions.

We utilized the Uniform Manifold Approximation and Projection algorithm, along with one-hot encoding for sequence representation, to reduce the dimensionality of the 12 single-point mutations combinational space and visualized it as a two-dimensional scatter plot. Our analysis revealed that mutations with similar activities clustered closely together (Fig. 1h). The mutants from the training set were distributed across the entire sequence space, which was crucial for developing an accurate machine learning model (Fig. 1h). Furthermore, our model’s predictions aligned well with the experimental data with high-activity variants in the outer clusters and low-activity variants in the inner clusters (Fig. 1h). This indicates that the model accurately mapped the fitness landscape, allowing us to identify mutants with significantly increased activity in the high-activity clusters.

Further improvement of IFRS activity by deep learning models

To further improve IFRS activity, we explored additional mutation sites in TBD with Com1-IFRS as the template, by employing three deep learning models that enable zero-shot prediction of high-fitness variants, including ESM-1v, MutCompute, and ProRefiner (Fig. 2a, Supplementary Fig. 7 and Table 5). EMS-1v is a protein language model that was trained on extensive datasets of protein sequences spanning the evolutionary tree of life, and hence learned the fundamental principles of protein structure and function²¹. Given the sequence of Com1-IFRS, each amino acid in the TBD region was individually masked and analyzed by the ESM-1v model to predict the impact of potential mutations at that site. We constructed and characterized 16 single-point mutants that were predicted to have a higher likelihood than wild type (Supplementary Figs. 7 and 8). MutCompute is a structure-based three-dimensional self-supervised convolutional neural network model that was trained to associate local protein micro-environments with their central amino acid²². Given the N-terminal and C-terminal structures of MmPylRS (5UD5.pdb & 4TQD.pdb), MutCompute predicted single-point mutations optimizing the protein structure, and we characterized the top 44 mutants based on the probability predicted (Supplementary Figs. 7 and 8). ProRefiner is a global graph attention model for inverse protein folding that designs sequences compatible with a given backbone structure²³. We restrict the candidate mutation sites to the TBD of PylRS and leverage sequence design models to compute a quality score for every candidate site. For each site to be examined, we masked this site in the sequence to get the input partial sequence, and the input backbone structure is from the Com1-IFRS structure predicted by Alphafold 3²⁴. The model then predicted the identity of the masked site in the form of a probability distribution over all amino acid types, and the top 42 single-point mutants were selected for characterization (Supplementary Figs. 7 and 8). Interestingly, 7 mutants were predicted to have improved activities by both Mutcompute and ProreRiner (Supplementary Fig. 7). Based on three methods, a total of 95 single-point mutants were constructed on 85 amino acids of COM1-IFRS and assayed for enzyme activity (Fig. 2b and Supplementary Fig. 7).

Since the enzyme activity is represented by the fluorescence intensity of sfGFP incorporated with 3BrF, the maximum activity that could be detected is the fluorescence intensity of wild-type sfGFP. The fluorescence intensity of sfGFPS2TAG in the presence of Com-1-IFRS reached 45% of wild-type sfGFP. Among all the single-point variants constructed, we did obtain several variants with higher activities than the Com1-IFRS (Fig. 2b), which were I176S predicted by ESM-1v, D76A, T68V, H28K, T20S predicted by MutCompute, K67S, N7S, V74I predicted by ProRefiner, and H63N predicted by both MutCompute and ProRefiner. The best variant D76A, designed by MutCompute, showed a 31% improvement in activity compared to Com-1-IFRS. However, most of the variants designed by the three approaches exhibited reduced or even lost activity, such as W16E, showing a 99.7% decrease in activity compared to IFRS.

We then wondered if these data are useful for training a supervised ML model to predict high-fitness single-point variants across the protein. With the above single-point mutants’ activity data as the training set (N = 96, including Com1-IFRS), the FFT-PLSR model was built. There are 566 amino acid encodings in the AAindex database, and we tested the performance of the models trained with different numbers of amino acid encodings. To optimize the amino acid encoding, we first screened the single AAindex encoding, which was then fixed for optimizing the second encoding. The third encoding was optimized with the first two encodings fixed. When screening two or three indices, the protein sequence was first subjected to FFT separately, and then the results were combined to train the model. When one index was used, the R² of the model was only 0.452, while when three amino acids encodings were used, the R² of the model increased to 0.926 (Supplementary Fig. 9). We then used the three-index model to predict fitness of all single-point mutations in the PylRS TBD region, and the top 20 variants were constructed and characterized for enzyme activity (Fig. 2b, Supplementary Fig. 10). The best variant, K67G, showed a 31.9% improvement in activity, which was even higher than D76A predicted by MutCompute and K67S predicted by ProRefiner, indicating that the FFT-PLSR model was effective in exploring sequence space (Fig. 2c).

To further explore the sequence space, we constructed the saturation mutagenesis on the nine sites where the improved variants were obtained, to build a mutability landscape (Fig. 2d). The mutability landscape is defined by the impact of all possible point mutations on protein function by substituting the native amino acid at each residue position with each of the 19 non-native amino acids, one at a time^25,26. The mutability landscape showed that most of the improved variants were found at positions D76, H63, and K67, and the mutations at sites T20, H28, and I176 were mostly detrimental. We screened the variants that showed over a 10% improvement in SCS efficiency compared to Com1-IFRS (Fig. 2e). There were 3, 7, 5, 1, 2, 9 mutations at positions N7, H63, K67, T68, V73, D75, respectively, meeting this criterion. We hence attempted to use the FFT-PLSR model to explore this sequence space containing 11,520 (4 × 8 × 6 × 2 × 3 × 10) mutations. To achieve epistasis information, the top 3 variants at each mutation site were combined to construct double variants, resulting in a total of 92 combined variants to be used for model training (Fig. 2f). Combination of improved single variants did generate further improved variants, with the best double mutant of N7E/T68F showing a SCS 70% higher than the Com1-IFRS (Fig. 2f). However, the variants with decreased activities were also observed, indicating a strong epistasis among several mutations. For example, both V74W and K67G improved the SCS efficiency, while V74W/K67G showed a significant decrease compared to the Com1-IFRS, which resulted in a strong negative reciprocal sign epistasis effect between V74W and K67G. The antagonistic effect was also observed for V74W and single variants including K67L, D76F, D76L, and D76Y (Supplementary Table 6).

We then conducted another round of FFT-PLSR building for predicting high-activity combinatory mutations. There were 27 single-point mutations and 92 combinatorial mutations, plus Com1-IFRS, 120 data points in total in the training dataset. These were used to train the FFT-PLSR model. We first used 10-fold cross-validation to evaluate and identify the optimal single-index model. However, this model only achieved an R² score of 0.558 on the training set, indicating a poor performance. To improve the model performance, multiple indices were used for encoding amino acids, which improved the R² score to 0.668 for the two-index model and 0.677 for the three-index model, respectively (Supplementary Fig. 11 and Table 7). The retrained FFT-PLSR model achieved an R² of 0.729 on the test set, which consisted of 15 double-point combinatory mutants not included in the training set (Fig. 2g, Supplementary Fig. 12). Based on the three-index model, we predicted the activity data of 11520 mutations, and the top 20 combinatorial variants were selected for experimental verification (Supplementary Table 8). All the variants showed comparable or higher activities than the Com1-IFRS, with 15 of them possessing activities improved by more than 2-fold, suggesting the great effect of the ML model. The best Z7-3 variant (N7Y/H63L/K67N/V74W, Com2-IFRS) showed a 2.8-fold increase in SCS efficiency compared to Com1-IFRS, more than 30-fold higher than IFRS, reaching fluorescence intensity comparable to wild-type sfGFP (Fig. 2h). Biochemical characterization of IFRS, Com1-IFRS and Com2-IFRS using 3BrF confirmed that the catalytic efficiency (k_cat/K_m) of Com1-IFRS and Com2-IFRS for tRNA was improved by 1.4-fold and 5.6-fold, respectively, compared to IFRS (Supplementary Table 9). Additionally, we tested if the mutations H20S, H28K, I176N, and I176S, which were not selected for making combinatory mutations, could further improve the activity of Com2-IFRS. It was found that none of these single-point mutations improved the activity of Com2-IFRS (Supplementary Fig. 13).

Combination of tRNA-binding domain mutations with catalytic domain mutations to generally enhance the incorporation efficiency of diverse ncAAs

As the IFRS is polyspecific and could accept various phenylalanine derivatives²⁷, we hence tested if Com1-IFRS and Com2-IFRS improved the incorporation efficiency of other ncAA substrates. It was found that normalized fluorescence of sfGFPS2TAG incorporated with 12 diverse ncAAs was significantly increased by the two variants compared to IFRS, with the largest 101.9-fold improvement for 3FF enabled by Com2-IFRS (Fig. 3a). However, since IFRS is a promiscuous enzyme with a certain degree of misincorporation of canonical amino acids (cAAs), Com1-IFRS and Com2-IFRS also increased the incorporation efficiency of cAAs, and different extent of misincorporation was observed in presence of different ncAAs (Fig. 3a). Subtracting fluorescence intensity of sfGFP2TAG in presence of 3FF with that in absence of 3FF revealed a 3944.8-fold improvement in SCS efficiency for Com2-IFRS compared to IFRS (Supplementary Fig. 14a). Biochemical characterization of IFRS, Com1-IFRS and Com2-IFRS using 3FF confirmed that the catalytic efficiency (k_cat/K_m) of Com1-IFRS and Com2-IFRS for tRNA was improved by 1.8-fold and 8.8-fold, respectively, compared to IFRS (Supplementary Table 9).

CTD of PylRS containing the catalytic sites has evolved to accept various ncAAs. The mutations Com1 and Com2 were then combined with different catalytic domain (CD) mutations to enhance the incorporation efficiency of the corresponding ncAAs. To test the universality of Com1 and Com2 in enhancing the incorporation efficiency of various ncAAs, we selected CD mutations with large sequence diversity that could accept ncAAs with significantly different side chains, including Phe derivatives, Tyr derivatives, Trp derivatives, Cys derivatives, His derivatives, and Lys derivatives (Supplementary Table 10). CD mutant NACA was combined with Com1 and Com2 to test the incorporation of three Phe derivatives, including 3-bromo-L-phenylalanine (3BrF),2-chloro-L-phenylalanine (2ClF), and 3-L-phenyllactic acid (PLA). The two variants significantly improved efficiency of NACA to incorporate these three ncAAs, and the Com2 led to improvement of 38.9-fold, 29.2-fold and 7.7-fold, for 3BrF, 2ClF, and PLA respectively (Fig. 3b). When the fluorescence intensity was normalized to cAA incorporation, the fold change was 61.1-fold, 68.1-fold and 9.9-fold, for 3BrF, 2ClF, and PLA, respectively (Supplementary Fig. 14b). IFRS could also accept 3BrF and 2ClF, and Com1, Com2 mutations significantly improved the misincorporation of cAAs in presence of 3BrF and 2ClF (Fig. 3a). Interestingly, when combined with NACA, the two variants did not improve the incorporation efficiency of cAAs in presence of 3BrF and 2ClF, compared to IFRS (Fig. 3b). This could be due to that NACA exhibited lower activity against cAAs than IFRS, and Com1, Com2 did not influence substrate specificity and generally enhanced the activity of CD mutations against all substrates.

This was further confirmed by the performance of Com1 and Com2 introduced in wild-type MmPylRS. MmPylRS showed limited misincorporation of cAAs, and the two variants did not increase the misincorporation of cAAs either, compared to the wild type. On the other hand, the two variants significantly increased the incorporation efficiency of N6-(tert-butoxycarbonyl)-L-lysine (BocK) and N6-((allyloxy)carbonyl)-L-lysine (AlocK) and the best Com2 resulted in 40.2-fold and 32.7-fold improvement for BocK and AlocK, respectively (Fig. 3b, Supplementary Fig. 14b). This was further confirmed by the huge increase in expression level of sfGFP incorporated with BocK in presence of Com2-WT compared to WT (Supplementary Fig. 15). Mass-spectra analysis also confirmed the correct incorporation of BocK and no misincorporation of cAAs was observed (Supplementary Fig. 16 and Table 11). Interestingly, the incorporation efficiency of BocK directed by Com2-WT was 2.7-fold higher than chPylRS-IPYE obtained previously⁷. Additionally, the SCS efficiency of Com2-WT against BocK was nearly 17.9-fold higher than MaPylRS, the PylRS enzyme lacking an NTD, although a previous study showed that MaPylRS exhibited higher activity than wild-type MmPylRS¹⁰ (Supplementary Fig. 17a). Biochemical characterization of WT and Com2-WT using BocK confirmed that the k_cat of Com2-WT improved 2.6-fold, the K_m for BocK weakly increased, such that the catalytic efficiency (k_cat/K_m^BocK) of the evolved variant was enhanced by 2.2-fold compared to wild-type MmPylRS (Supplementary Fig. 18). The kinetic parameters were also measured for tRNA, and the catalytic efficiency (k_cat/K_m^tRNA) of Com2-WT was improved by 1.8-fold compared to wild type (Supplementary Table 9). Interestingly, K_m for tRNA was increased by 1.4-fold for Com2-WT than wild type. Similarly, the binding affinity of Com2-WT with tRNA was around 2.1-fold lower than that of WT, suggesting that the mutations on the TBD improved the enzyme activity by modifying the tRNA binding conformation instead of enhancing the binding affinity (Supplementary Fig. 19).

Com1 and Com2 also enhanced the activities of other CD mutations toward their corresponding ncAAs. The addition of Com2 increased the incorporation efficiency of N-epsilon-Acetyl-L-lysine (AcK) by 13.3-fold compared to the original CD mutations MLAF, and the fold change was 32.2-fold when the incorporation efficiency was normalized to cAA incorporation (Fig. 3b, Supplementary Fig. 14b). The extent of improvement was significantly higher than the IPYE variant tested previously⁷. The SDS-PAGE also revealed that the amount of sfGFP expressed in presence of Com2-MLAF was significantly higher than that in presence of MLAF (Supplementary Fig. 15). Mass-spectra analysis confirmed the correct incorporation of AcK in sfGFP, but a misincorporation of lysine was also observed (Supplementary Fig. 16). Additionally, addition of Com2 improved incorporation efficiency of GML mutant by 5.4-fold, 5.8-fold and 3.9-fold against 3BrY, 3ClY and 3IY, respectively, while the fold change reached 117.3-fold, 587.6-fold and 6.5-fold when the SCS efficiency was normalized to cAA incorporation (Fig. 3b, Supplementary Fig. 14b). We checked expression of sfGFP with 3IY incorporated, and found that Com2-GML indeed significantly increased amount of expressed protein compared to GML. Mass-spectra analysis also confirmed the correct incorporation of 3IY (Supplementary Fig. 16).

Com1 and Com2 mutations were also constructed in BtaRS to explore their effect on the incorporation of the Trp derivatives²⁸. The addition of Com2 increased the incorporation efficiency of 3-benzothienyl-L-alanine (Bta) and 3-(1-naphthyl)-L-alanine (1NaA) by 56.8-fold and 63.3-fold, respectively (Fig. 3b). SDS-PAGE revealed the significantly improved amount of sfGFP with Bta incorporated enabled by Com2-BtaRS compared to BtaRS. Kinetic parameters measurement using Bta revealed that the catalytic efficiency (k_cat/K_m) of Com2-BtaRS for tRNA was 4-fold higher than that of BtaRS (Supplementary Table 9). Mass-spectra analysis confirmed the correct incorporation of Bta and no misincorporation of cAAs was observed (Supplementary Table 11). However, misincorporation of cAAs was indeed observed when 1Na was incorporated by Com2-BtaRS (Fig. 3b). To explore the effect of Com1 and Com2 on the incorporation of Cys derivatives, Com1 and Com2 were combined with CD mutation WS. Com1 did not exhibit improved effect on fluorescence of sfGFP, while the addition of Com2 improved fluorescence intensity of sfGFP with Sac incorporated by 2.8-fold compared to WS (Fig. 3b), and the enhanced amount of sfGFP expressed was confirmed on SDS-PAGE (Supplementary Fig. 15). WS also showed a certain degree of misincorporation of cAAs, which was also observed for Com2-WS (Fig. 3b).

Com1 and Com2 were combined with two CD mutations, including QF and IFGFF, to explore their effect on incorporating His derivatives, including 3-(2-thienyl)-L-alanine (2ThiA), 2-(5-bromothienyl)-L-alanine (BrThiA), and 3-methyl-L-histidine (3MeH). In the results, the addition of Com2 increased the incorporation efficiency of 2ThiA, BrThiA, and 3MetH by 93.7-fold, 41.5-fold, and 40.9-fold, respectively, compared to their CD variants. The fold change reached to 223.1-fold, 61.4-fold and 201.5-fold, respectively, when the amber codon suppression activity was normalized to cAA incorporation (Fig. 3b, Supplementary Fig. 14b). SDS-PAGE revealed that Com2-QF and Com2-IFGFF indeed dramatically enhanced the amount of purified sfGFP incorporated with 2ThiA and 3MetH, respectively, compared to CD mutations alone (Supplementary Fig. 15). Kinetic parameters measurement using 3MetH confirmed that the catalytic efficiency of Com2-IFGFF for tRNA was improved by 4.3-fold compared to IFGFF (Supplementary Table 9). Moreover, according to mass-spectra analysis, no misincorporation of cAAs was observed for GFP-2MetH and GFP-2ThiA (Supplementary Fig. 16).

To explore if the variants obtained were useful in improving the expression of other proteins with ncAAs incorporated, we tested the expression of myoglobin containing 3MetH (Fig. 4a). 3MetH has been used as a heme ligand to enhance the activity of myoglobin or as a catalytic residue of an artificial esterase possessing a non-canonical organocatalytic mechanism²⁹. Here, we used MmPylRS variant Com2-IFGFF to incorporate 3MetH into the position His93 of myoglobin as a ligand of heme. The protein expression was explored by carrying out protein purification from the same amount of cells in the presence of Com2-IFGFF and IFGFF. It was found that the concentration of 3MetH-containing myoglobin was 28.3 mg/L for Com2-IFGFF, 6.3-fold higher than 4.5 mg/L of IFGFF (Fig. 4b). We then measured the activities of Mb-3MetH against guaiacol using the purified protein without dilution. The myoglobin catalyzes the oxidation of guaiacol by hydrogen peroxide to generate a stable tetrameric product whose formation can be readily monitored by absorbance at 470 nm. The yield of product was significantly higher for Com2-IFGFF compared to IFGFF (Supplementary Fig. 20). The ΔOD₄₇₀ reached 0.063 in the reaction system containing Com2-IFGFF after a 40-min reaction, 7.9-fold higher than that using IFGFF (Fig. 4b, Supplementary Fig. 20), confirming the higher expression of target protein aided by the Com2-IFGFF variant.

Suppression of multiple amber codons of PylRS variants

We then characterized the ability of Com1-IFRS and Com2-IFRS to suppress multiple amber codons in sfGFP, which is important for incorporating multiple unnatural amino acids into proteins. One to five consecutive amber codons were inserted second position of sfGFP (Fig. 4c). Both Com1-IFRS and Com2-IFRS exhibited higher fluorescence intensity than IFRS in all situations, while Com2-IFRS showed higher suppression ability than the Com1-IFRS, with 122.4-fold, 99.4-fold, 91.2-fold and 53.3-fold improvement compared to IFRS, for S2TAG × 2, S2TAG × 3, S2TAG × 4 and S2TAG × 5, respectively (Fig. 4d, Supplementary Fig. 21a). It was also found that Com1-IFRS and Com2-IFRS improved the incorporation efficiency against native amino acids compared to IFRS. We also tested the incorporation efficiency of multiple unnatural amino acids at different positions of sfGFP (Fig. 4e). When 3BrF was incorporated at the position of D36 of sfGFP, the fluorescence intensity was different from that of sfGFP with 3BrF at the second position, for both the wild-type IFRS and mutant Com1-IFRS and Com2-IFRS (Fig. 4f, Supplementary Fig. 21b). Similar site-dependent incorporation efficiency has previously been observed for other ncAAs⁷. Despite this, Com2-IFRS still showed higher amber codon suppression efficiency than Com1-IFRS and IFRS, with 3.8-fold, 7.9-fold, 27.3-fold, 4.7-fold and 5.2-fold improvement compared to IFRS, for 1TAG, 2TAG, 3TAG, 4TAG and 5TAG, respectively.

MD simulations to explore the molecular change of PylRS variants

The whole 3D structure of MmPylRS was not yet available as the full-length protein is insoluble. Hence, AlphaFold3 was used to predict the structures of MmPylRS (WT), Com1-WT, and Com2-WT in complex with tRNA^Pyl (Supplementary Fig. 22). The predicted MmPylRS structure aligned well with the separate NTD structure and CTD structure determined previously (Fig. 5a). 50-ns MD simulations were then conducted for these structures in complex with Pyl-AMP to understand how the mutations influenced the binding of tRNA and the enzyme activity (Supplementary Figs. 23–25). The root-mean square deviation (RMSD) values revealed that the trajectories were well equilibrated at the last 10 ns, which were used for further analysis (Supplementary Fig. 25a). The reaction distance between the 3′-OH of tRNA A76 and the carboxyl carbon atom of amino acid Pyl was first analyzed. It was generally shorter for Com2-WT than wild type and Com1-WT (Fig. 5b, Supplementary Fig. 26). A total of 1000 snapshots were analyzed, and the number of snapshots with a distance shorter than 4 Å was 462, which is 21-fold and 10-fold higher than that of wild type and Com1-WT, respectively (Fig. 5c). These indicated that Com2 mutations mediated the binding of tRNA^Pyl to make the aminoacylation reaction happen more easily. Interestingly, in the reaction conformations, we observed new hydrogen bonds formed between Pyl-AMP and tRNA in the two variants compared to WT. Com1-WT showed a new hydrogen bond formed between the main chain -NH₂ of Pyl and 2′-OH of tRNA A76, while Com2-WT exhibited a hydrogen bond formed between the main chain -NH₂ of Pyl and the 3′-OH of tRNA A76. No such hydrogen bonds were found in wild type (Supplementary Fig. 27). These new hydrogen bonds will contribute to the interaction between Pyl-AMP and tRNA, and hence accelerate the reaction.

The hydrogen bonds formed between tRNA and the protein were also analyzed. In last 10-ns MD simulations, the number of hydrogen bonds in PylRS TBD with occupancy over 60% was 22, 19, and 22 for WT, Com1-WT, and Com2-WT, respectively (Supplementary Fig. 28). Specifically, both Com1-WT and Com2-WT formed new hydrogen bonds including LYS3-A58, ARG19-A46, ARG52-G52, ARG193-A5, ARG193-C13 and ARG193-U12, while Com2 formed several extra hydrogen bonds such as ARG55-C45, ARG55-A46, Arg58-A20 (Fig. 5d). Additionally, several hydrogen bonds were disrupted in the variants, such as ASN49-G47, ARG55-A46, ARG55-G21, ARG58-A58, R66-G21 and so forth. Specifically, in the conserved Motif 2 loop that is responsible for tRNA recognition, two hydrogen bonds Lys336-C71 and Lys336-C72 were disrupted, and a new hydrogen bond Asp334-C71 was formed in both the Com1 and Com2 variants (Supplementary Fig. 29). Additionally, an extra H-bond GLU332-C74 was formed in Com1 variant but not in WT and Com2. This indicated that mutations reshaped the interactions between tRNA and protein, and the Com1 and Com2 improved the tRNA binding in a different way. As a result, the interaction energy between tRNA and protein was different for the wild type and variants.

The binding free energy between tRNA and different domains of the protein was analyzed. Com1-WT and Com2-WT exhibited lower binding free energy compared to the wild type, which was mainly attributed to the decreased binding free energy of the tRNA binding domain (Supplementary Fig. 30). Interestingly, the binding free energy of full-length Com1-WT was lower than Com2-WT, while that of the tRNA binding domain of Com1-WT was higher than Com2-WT. The binding free energy determined in MD simulations seemed to contradict with binding affinity and K_m values measured for WT and Com2-WT. This could be attributed to the different ncAA substrate used. In the MD simulations, the substrate was native substrate pyrrolysine, while the substrate used in biochemical characterization was BocK. Kinetic parameters measurement using different ncAA did reveal that Com2 influenced the K_m for tRNA in a different way in the presence of different ncAA (Supplementary Table 9).

Residue-level binding energy contribution analysis for both protein and tRNA was carried out. Several amino acids and tRNA bases were indeed found to impact the binding free energy. For example, GLU332 in Com1-WT exhibited reduced binding energy, while no significant changes were observed in WT or Com2-WT, which might be attributed to the newly formed H-bond Glu332-C74 in Com1-WT. S193R mutation formed new salt bridges with tRNA, including Arg193-A5, Arg193-C13, Arg193-U12, which led to a significant decrease in the binding free energy (Fig. 5e). Also, Arg55 of Com2-WT showed a significantly low binding energy due to the T56P mutation, although it is not the case for Com1-WT. As for the analysis of tRNA binding energy, it was found that the mutations significantly increased the binding free energy for the last three bases C74, C75, and A76, which might facilitate the aminoacylation reaction (Fig. 5e).

Dynamics cross-correlation matrices (DCCMs) were also computed for the WT and two variants to understand how the mutations impact protein dynamics. Since the MD simulation systems are large and complex, a robust coupling analysis of dynamic cross-coupling correlation might need simulations of a longer timescale than 50 ns. We hence carried out 200-ns simulations for WT, Com1-WT, and Com2-WT (Supplementary Figs. 31–33, Supplementary Note 1), and the last 150 ns trajectories were used to compute DCCM of the protein C_α atom pairs and tRNA P atom pairs. Generally, the variants showed more dynamics cross-correlations between residue pairs compared to the wild type, while Com1-WT and Com2-WT exhibited similar C_ij values in most of the regions (Fig. 5f). Specifically, it was found that the dynamics correlation between NTD and CTD was more significant in the Com1-WT and Com2-WT, compared to WT. Although NTD and CTD are distant from each other, they are connected together with the tRNA. The mutations on tRNA binding domain hence impact the dynamics of CTD through modification of interactions with tRNA. The increased correlated dynamical network would help to maintain more reaction conformations, thereby enhancing the enzyme activity of Com1-WT and Com2-WT.