Mechanism-Based Strategy for Optimizing HaloTag Protein Labeling
Highlights• A comprehensive kinetic and computational study provided detailed insights into the mechanism of HaloTag ligand incorporation.• Both the ligand binding and the subsequent chemical step can be decisive for the efficiency of covalent ligand incorporation.• Natural dehalogenases provide high efficiency in ligand incorporation with no need for intensive optimization using directed evolution.• A simple strategy is proposed to identify the optimal tag from a wide range of natural dehalogenases.
Introduction<!>A B<!>Results<!>Protein expression and purification<!>Kinetic analysis<!>Computational studies<!>Discussion<!>Conclusions<!>Materials and Methods
<p>Genetically encoded protein labeling methods are widely employed in protein chemistry, molecular and cellular biology. More recently, self-labeling protein tags designed for covalent conjugation to small-molecule ligands functionalized with biorthogonal linkers have gained widespread attention. One of the most popular self-labeling methods, HaloTag®, uses engineered haloalkane dehalogenase (HLD) genetically fused to the proteins of interest, which covalently binds synthetic ligands bearing various functionalities, such as a strong light-up fluorescence response. The original concept of bifunctional linkers developed by Dick Janssen and co-workers 1 for covalent capturing and ribosomal/phage display of HLDs was translated to in vivo and in vitro analysis of mammalian proteins by Wood and co-workers [2][3][4] .</p><p>Since its development and commercialization by Promega, HaloTag® has become a valuable research tool for a broad range of applications (Figure 1A) including protein purification 5 and immobilization 6 , enhancement of the soluble expression of recombinant proteins 7 , cellular protein imaging 8,9 , imaging in vivo 10 , and single-molecule studies [11][12][13] . The technology is applicable for analyses of protein-protein and protein-nucleic acid interactions 14,15 , proteome stress 16,17 , protein folding and aggregation 18,19 , dynamics and hydration [20][21][22] , or cellpermeability 23 . HaloTag fusions enable protein control in vivo [24][25][26] , including degradation 27,28 or dimerization 29 of proteins of interest. Further applications include high-throughput screening methods, microarrays, and chip technology [30][31][32] , intracellular detection of pH 33 or biologically important ions 34,35 , mechanochemistry 36,37 , functionalization of nanoparticles 38 , and quantum dots 39 . Recently, the potential of the HaloTag technology in cell therapy was discovered, as it has been used for cell surface modification permitting angiogenesis, increased motility, and immune shielding 40 . A wide range of diverse HaloTag ligands have been designed and synthesized offering a variety of properties (Figure 1B), e.g., improved photostability and brightness 12 , high biocompatibility or fluorogenicity allowing "no-wash" labeling protocols [41][42][43] or providing specific affinity handles 4 . Despite the great diversity of ligands used, the majority of their applications always utilize the same tag protein DhaAHT, without considering the choice of another protein partner for better recognition of a specific ligand. The labeling reaction proceeds via a twostep kinetic pathway, the binding of the ligand, and the following chemical conversion, leading to the stable covalent alkyl-enzyme complex. The latter unimolecular step cannot be easily optimized by modifying the labeling protocol and it depends solely on the optimal reactive orientation of the bound ligand. The architecture and dynamics of the protein access tunnel and active site thus plays a key role. The 10,000-fold improvement in binding efficiency of the DhaAHT tag required for the successful protein imaging was achieved using a focused directed evolution on the access tunnel residues 44 . This study has already shown that the efficiency obtained by molecular evolution differs significantly among individual ligands, despite sharing the same reactive linker. Similar effects were observed in our recent studies focused on the engineering of access tunnels in HLDs 45,46 . The binding efficiency of HaloTag ligands varied across seven orders of magnitude for HLDs with different architectures of their access tunnels. Interestingly, the change in the functional reporter strongly affected the labeling efficiency even for the ligands with the same reactive linker 45,46 . The results collectively suggest that the broadly used DhaAHT tag may not be the optimal tag for the incorporation of various available ligands. Since DhaAHT was introduced, the portfolio of available dehalogenases has significantly expanded 47,48 and currently offers an interesting range of variants (Figure 2, Table 1).  The main tunnel calculated with CAVER 3.02 50 in the respective structures after adding hydrogens atoms using PyMOL 2.3.2 51 ; the tunnel origin is defined at the carboxylic-O atoms of Asp106; a probe radius is 0.6 Å; shell depth 4 Å; and shell radius 4 Å; b the mutations affecting the p1 tunnel are marked with asterisk; NA means not applicable; c all tunnels are shown from the same viewpoint, being the active site located at the lower end of the tunnels and the protein surface at the top.</p><!><p>In this study, we present a comprehensive kinetic and computational study of a mechanism of the HaloTag ligand incorporation. We investigated the effect of two different functional reporters, tetramethylrhodamine (TMR) and 4-stilbazolium, and different lengths of the reactive linker. We compared DhaAHT optimized by the directed evolution with three natural dehalogenases DhaA, LinB, and DmmA. The results have shown striking differences in the kinetics of both of the involved reaction steps, the ligand binding and the consecutive chemical step. The efficiency of each step is different for individual protein-ligand pairs. Strikingly, the most efficient reaction was not determined for DhaAHT and TMR, albeit systematically optimized by directed evolution. The wild-type enzymes LinB and DmmA showed the highest efficiency of the probe incorporation with the 4-stilbazolium probes. Our study proposes a new concept of selecting the protein tag matching the specific ligands.</p><!><p>We performed comprehensive kinetic and theoretical study for the incorporation of two different HaloTag ligands, commercial tetramethylrhodamine (TMR) and three 4stilbazolium-based ligands (1B, 1D, and 1E) with different lengths of the reactive linker (Figure 3). The 4-stilbazolium-based dyes have shown a stronger fluorogenic response upon labeling and easier synthetic routes than any of the previous HaloTag labels, which could be highly beneficial for their applicability 42 . We compared DhaAHT optimized by directed evolution with analogs of three natural dehalogenases DhaA, LinB, and DmmA. DhaAH272F, LinBH272F, and DmmAH315F, correspond to the natural enzymes with a single additional mutation in the catalytic base (H272F for DhaA and LinB, H315F in DmmA). This mutation in the catalytic histidine leads to the interruption of the catalytic cycle by preventing the hydrolytic step, and to the formation of the covalent alkyl-enzyme intermediate as the final complex.</p><!><p>The haloalkane dehalogenases genes linBH272F, dhaAH272F, dhaAHT, and dmmAH315F were cloned into pAQN, pET21b or pET24a vectors and transformed into Escherichia coli BL21 (Supplementary Table S1). The enzymes were overexpressed and purified by metal-affinity chromatography. The purity of proteins was analyzed by SDS-PAGE (Supplementary Figure S1). A detailed description of the gene cloning, expression, and purification are provided in the Supplementary Information (Section I).</p><!><p>Fluorescence intensity and anisotropy were used to systematically analyze the concentration dependence of DhaAHT, DhaAH272F, LinBH272F, and DmmAH315F reactions with TMR and 1B, 1D, and 1E. The conventional fitting and numerical integration methods were applied to obtain detailed information about the individual rate and equilibrium constants related to the two-step model of the HaloTag reaction (Figure 3).</p><p>Initially, the kinetic data were analyzed by conventional exponential fitting using nonlinear regression. To compare the consistency of the data with earlier published results, the apparent second-order rate constants were calculated following the procedure used originally by Los and co-workers 4 . The value of the apparent rate constant obtained for the incorporation of TMR into DhaAHT, 2.3 × 10 6 M -1 •s -1 , corresponds well with the value reported by Los and co-workers of 2.7 × 10 6 M -1 •s -1 . 4 Next, the concentration dependence of the kinetic data (Figure 3A) was explored to provide detailed information about the kinetic pathway and the estimates of the true rate and equilibrium constants. Although the single-exponential fit of DhaAHT traces obtained with TMR provided satisfying statistics (χ2/DoF = 2.29; p-value = 0.28), the use of a double-exponential function showed significantly improved goodness of fit (χ2/DoF = 1.21; p-value = 0.43) and distinguished two separate phases, well consistent with expected two-step kinetic model for the reaction (Figure 3C).  52 The kinetic model of HaloTag reaction: E is the enzyme, L is the ligand, E.L is enzyme-ligand complex, E-L is the covalent alkyl-enzyme complex, k1 and k-1 is the rate of association and dissociation of enzyme-ligand complex, respectively, and k2 represents the rate of chemical step (nucleophilic substitution SN2). (D) The kinetic parameters obtained by numerical analysis of anisotropy data. Error bars represents the standard error of the fitted parameters. The rigorous confidence contour analysis of variance of fitted parameters is presented in Supplementary Information (Supplementary Table S5). The kinetic experiments were performed in two to three independent replicates.</p><p>The concentration dependence of the obtained rates was analyzed analytically by a secondary fitting to the approximate rate equations derived for the two-step model (Supplementary Equations S4 and S5). The analysis provided initial estimates of the rate constant for association (k1 = 41 ± 4 μM -1 •min -1 ) and dissociation (k-1 = 0.08 ± 0.04 min -1 ) for TMR. DhaAHT bound complex and the rate constant for subsequent chemical step resulted in the formation of the covalent alkyl-enzyme complex (k2 = 0.06 ± 0.04 min -1 ). In the case of DhaAHT reaction with 4-stilbazolium-based ligands (Figure 3B), the binding phase gradually disappears in the dead-time of the measurement with increasing concentration of the enzyme, and only singleexponential fit provided reasonable estimates of the rates and amplitudes. The concentration dependence of the observed rate (Supplementary Equation S6) allowed to define the initial estimates of the equilibrium dissociation constant for the enzyme-ligand bound complex (KD = k-1/k1) ranging from 0.6 to 1.3 μM and the rate of consecutive chemical step (k2) ranging from 0.03 to 0.33 min -1 for individual 4-stilbazolium ligands.</p><p>The conventional analysis of the kinetics of DhaAHT indicated substantial differences in the reaction with the two types of tested ligands. The TMR probe shows a rapid binding to a tight enzyme-ligand complex followed by relatively slow chemical conversion leading to the final covalently bound complex. The accumulation of the reversible enzyme-ligand bound complex (E.L) is well visible from the anisotropic data (Figure 3A middle), which shows a significant concentration dependence of the equilibrium signal with amplitudes defined by the equilibrium dissociation constant of enzyme-ligand complex. The numerical simulation of the fraction of individual reaction species (Figure 3A right) illustrates the course of the reaction involving rapid binding of TMR into DhaAHT associated with the accumulation of enzymeligand bound complex (E.L, red), which is slowly transformed into the final covalent alkylenzyme complex (E-L, green) (Figure 3A right). In contrast, the kinetics of the 4-stilbazolium ligand reaction is dominated by the chemical step leading to the dominant accumulation of the final covalent complex (E-L). Anisotropic traces thus reach the same level of the signal in equilibrium, which is defined by the total ligand concentration (Figure 3B).</p><p>Although the conventional approach is currently the most widely used method of kinetic data analysis, it has several limitations, such as the loss of an important relationship between velocity and amplitude, or the accumulation of errors associated with successive calculation of a large number of temporary parameters to estimate a small number of relevant kinetic constants. To overcome these limitations, we performed a global data analysis based on numerical methods. The parameter estimates obtained by conventional analysis (Supplementary Figures S2 and Table S2) were used as initial starting values for the numerical fitting. A detailed description of the conventional and numerical analysis of the kinetics data including rigorous statistical assessment is provided in the Supplementary Information (Section II).</p><p>The global data fitting used numerical integration of the rate equations from an input model (Figure 3C) searching a unique set of kinetic parameters (Supplementary Figure S3 and Table S3) that explain the original raw data and produce a minimum χ2 value 53 . The observable anisotropy signal was defined as the sum of the contributions of each species to the total signal with scaling factors for each species (Supplementary Table S4). In addition to monitoring the standard errors and residuals, the global fitting of kinetic data allowed to perform a rigorous analysis of variance referred to as a confidence contour analysis 54 . This analysis confirmed the high quality of the global fit, with all obtained kinetic parameters being well constrained by the experimental data (Supplementary Table S5). In the same way, the complex kinetic analysis was performed systematically comparing DhaAHT with three nonoptimized natural variants DhaAH272F, LinBH272F, and DmmAH315F in the reaction with TMR and 1B, 1D, and 1E (Supplementary Figure S3 and Table S3). The specific rate constants defining velocity of individual reaction steps, ligand binding (k1) and chemical conversion (k2), as well as the overall labeling efficiency defined by K1.k2, the product of the equilibrium constant for the ground-state binding K1 = 1/KD = k1/k-1 and the rate of the consecutive chemical step k2, are summarized for each enzyme variant in Figure 3D.</p><p>Unlike TMR, which provides only the possibility of instrumentally more complex anisotropy/polarization measurements, the 4-stilbazolium-based ligands provide the additional advantage in the tracking of fluorescence intensity signal (Supplementary Figure S4), commonly available in most laboratories. The increase in fluorescence intensity observed upon the incorporation of ligands into the enzymes was 5-, 2-and 10-fold for 1B, 1D, and 1E ligand, respectively. Such increases provide sufficient signal for in vitro enzymology studies. Additionally, the strong signal change observed especially for the 1E ligand provides a promising alternative to TMR in no-wash applications for cell labeling experiments.</p><!><p>DhaAHT, DhaAH272F, LinBH272F, and DmmAH315F were modeled from the crystal structures by in silico mutagenesis with Rosetta 55 , and the respective access tunnels from the active site to the surface were calculated using CAVER 3.02 50 . The tunnels found in DhaAHT, DhaAH272F, LinBH272F, and DmmAH315F showed considerably different geometric properties, especially the main tunnel, p1 (Figure 4A and Supplementary Figure S5). P1 is considerably wider in DhaAHT (bottleneck radius of 1.68 Å) than in DhaAH272F (1.29 Å) or LinBH272F (1.35 Å), and slightly wider than in DmmAH315F (1.62 Å). This fact suggests a higher accessibility of the DhaAHT active site as compared to DhaAH272F and LinBH272F, which provides a first explanation for the generally higher binding rates of the probes to DhaAHT and DmmAH315F than to the other two proteins. Moreover, the orientation of p1 is very different in LinBH272F and DmmAH315F compared to that of the DhaA variants. This suggests that LinBH272F and DmmAH315F may have different chemical and geometric preferences for the ligands that they can bind, compared to DhaAHT and DhaAH272F.</p><p>To understand the large differences found in the experimental kinetic measurements, we selected two representative probes (TMR and 1E) and two proteins (DhaAHT and LinBH272F) to study their molecular binding in more detail. These systems were chosen because of the high binding specificity found among two of the corresponding pairs, i.e., DhaAHT with TMR and LinBH272F with 1E. The TMR and 1E ligands were modeled and then refined with the Density Function Theory, which provided the energy-minimized structures and the partial atomic charges (see Supporting Information, Section III, for details). The binding of both probes to DhaAHT and LinBH272F was studied by molecular dynamics (MD), using the adaptive sampling approach 56 . The simulations started with the probes located in the bulk solvent, and consecutive rounds (epochs) of multiple MD simulations were performed. According to the adaptive sampling method, the starting points for the new MDs in each epoch were chosen from the previously sampled states based on the distance between the reacting groups in the probe and the enzyme. Each system was simulated for a total time of 20 µs. Markov state model (MSM) analysis was performed to obtain the relevant kinetic ensembles describing the binding of the molecular probes to the active sites of the proteins. Four Markov states could describe well the binding process, consisting of one fully bound state, two intermediates, and a fully unbound state (Figure 4B and Supplementary Figures S9 and S10). The kinetic parameters were calculated for the transitions between the most unbound state and the fully bound state. The results showed that the estimated binding rates (DhaAHT: k1 = 2.99 ± 0.45  10 8 M -1 •s -1 for TMR, 1.09 ± 0.23  10 8 M -1 •s -1 for 1E; LinBH272F: k1 = 7.6 ± 2.5  10 7 M -1 •s -1 for TMR, 1.22 ± 0.21  10 8 M -1 •s -1 for 1E; see Supplementary Table S7) followed the exact same order as the experimental ones, and k1 was highest for DhaAHT with TMR, and lowest for LinBH272F with TMR. The computational and experimental rates, however, differed by several orders of magnitude, being higher for the theoretical values. Such discrepancy is not unprecedented 56,57 and will be discussed below. Regarding the unbinding rates (k-1; Supplementary Table S7), they were all slower than the binding, which is consistent with the majority of the experimental results obtained here, although the order was not strictly observed. The slowest unbinding was predicted for TMR with DhaAHT, while the experiments showed the lowest unbinding rate for TMR with LinBH272F. The predicted binding affinity, given by K1 = 1/KD, also partially followed the experimental trends, where DhaAHT and TMR were correctly predicted with the highest affinity. However, the predicted order among the other pairs was incorrect, where LinBH272F and 1E showed experimentally the second strongest affinity. Interestingly, the highest probability of the bound state (Pbound) was obtained for TMR with DhaAHT, with 0.260 ± 0.043, and the lowest was found for 1E with LinBH272F, with 0.047 ± 0.008 (Supplementary Table S7).</p><p>Inspecting visually the bound states, we found that, in every system, the probes used exclusively the p1 tunnel (Figure 4B and Supplementary Figure S10). We also observed that the 4-stilbazolium aromatic system of 1E was partially inserted in the tunnel, both in DhaAHT and LinBH272F. Conversely, due to its longer linker, TMR presented its aromatic moiety completely outside of the protein. Moreover, the bound conformations of TMR followed the natural orientation of the p1 tunnel in DhaAHT, while 1E followed the tilted orientation of the p1 tunnel in LinBH272F (Figure 4). This interesting finding is in line with the fact that 1E is a better binder with LinBH272F and DmmAH315F than TMR. It also suggests a higher complementarity of the 4-stilbazolium-based probes with the LinB and DmmA variants in comparison with TMR. Analyzing in more detail the interactions found in the LinBH272F-1E complex, we found that in the bound state the aromatic system formed close hydrophobic contacts with L150, V173, L177, and L248, located in the tunnel, and with L179, P245, A247, A271, located on the extension of the tunnel mouth. Interestingly, the negatively charged D147 residue also formed electrostatic interactions with the 4-stilbazolium system due to its delocalized positive charge. We expect that the many interactions and constraints of the aromatic system of 1B, 1D, and 1E within the access tunnel contribute to strong fluorescence effects upon binding, as Clark and co-workers 42 have previously suggested. between DhaAHT and TMR, obtained from an adiabatic mapping of the distance between the reacting atoms of the protein (D106-COO -) and the TMR probe (CH2Cl). ∆G ‡ is the activation barrier of the reaction, where the ground state (GS), transition state (TS), and ligand-enzyme covalent complex (CC) are depicted. TMR is shown as the magenta sticks, the chloride ion as the green ball, the nucleophile D106 and the halide-stabilizing residues N41 and W107 as grey sticks.</p><p>Next, we predicted the reactivity of the TMR and 1E probes towards DhaAHT and LinBH272F and compared them with the experimental results. We started by analyzing the pre-reactive complexes found during the respective MDs, hereafter termed as near-attack conformation (NAC). We estimated the constant of formation of this pre-reactive complex, KNAC, based on the total number of NACs found and the probability of the bound state, Pbound. Surprisingly, the highest KNAC was obtained for 1E with LinBH272F, and the lowest for TMR with DhaAHT (Supplementary Table S7), with a difference of nearly two orders of magnitude. This suggests that, despite the binding of TMR with DhaAHT being extremely fast, the probability of the system adopting a potentially reactive conformation is very low. In contrast, once 1E and LinBH272F reach the bound state, the reactive conformation is achieved much faster than for any of the other systems. This perfectly follows the trends of the experimental kinetic rates obtained of the chemical step (k2 in Figure 3 and Supplementary Table S7). We then applied a hybrid QM/MM adiabatic mapping to estimate the energy barriers of the SN2 reaction, ∆G ‡ SN2 (Figure 4C). The predicted ΔG ‡ SN2 values (Supplementary Table S7) showed the lowest activation barrier for LinBH272F-TMR (12.1 ± 1.9 kcal•mol -1 ), followed by LinBH272F-1E (13.8 ± 1.8 kcal•mol -1 ), and the highest barrier for DhaAHT-TMR (15.5 ± 1.3 kcal•mol -1 ). This indicates that once the NAC has been achieved, LinBH272F provides a better environment for performing the SN2 step with both of the probes than DhaAHT.</p><p>Finally, the KNAC and ΔG ‡ SN2 were combined (Supplementary Information, Section III) to estimate the overall activation energy of the second kinetic step, ΔG ‡ 2. This step 2 (Figure 3C) is a direct measure of reactivity, and the estimated and experimental values of ΔG ‡ 2 can be directly compared. As a result, the highest calculated ΔG ‡ 2 value was obtained for DhaAHT with TMR, 2.5 kcal•mol -1 above the second highest energy barrier of DhaAHT-1E and 4.4 kcal•mol -1 above the ΔG ‡ 2 value for LinBH272F-1E (Supplementary Table S7). This is in reasonably good agreement with the experimental values, where DhaAHT with TMR also showed the highest ΔG ‡ 2 value, 2.7 kcal•mol -1 above that of LinBH272F with 1E.</p><!><p>We conducted a detailed kinetic and computational study of the reaction of a tetramethylrhodamine-based (TMR) and three 4-stilbazolium-based ligands (1B, 1D, and 1E) with reacting linkers of different lengths, and the haloalkane dehalogenase optimized by directed evolution (DhaAHT) and the analogs of three natural dehalogenases (DhaAH272F, LinBH272F, and DmmAH315F). The kinetic study showed substantial differences in the reaction kinetics between the individual enzymes as well as among the different ligands. The TMR probe showed a very high rate of binding towards the engineered DhaAHT (k1 = 39.7 ± 0.6 μM -1 •min -1 ), which is in high contrast with the other studied enzymes, which showed rates of binding ranging between 0.001 and 0.4 μM -1 •min -1 . In comparison to the non-optimized DhaAH272F, the engineering of the access tunnel of DhaA 44 led to a significant improvement of the ligand binding, while it did not compromise the catalytic step. The engineered DhaAHT showed even a slightly decreased activation barrier of the chemical step (ΔΔG ‡ 2 = -0.6 kcal•mol -1 ), although the ground state energy (ΔΔG0) of the enzyme-ligand complex was significantly lower (-3.3 kcal•mol -1 ) in comparison to DhaAH272F (Supplementary Table S6). Interestingly, the non-optimized LinBH272F showed a slow binding towards TMR, but it displayed the highest rate of the consecutive chemical conversion leading to the formation of the covalent enzyme-TMR complex.</p><p>The importance of the chemical step for the efficiency of the HaloTag labeling reactions is more pronounced with the 4-stilbazolium-based ligands. Even though the binding of 1E into LinBH272F is orders of magnitude slower and weaker, the elevated velocity of the following chemical step ensures the fully comparable labeling efficiency of 1E with LinBH272F (K1.k2 = 3.0 μM -1 •min -1 ) to that observed for the reaction of TMR with the engineered DhaAHT (K1.k2 = 3.1 μM -1 •min -1 ). The weaker binding of 1E to LinBH272F (ΔΔG0 = 2.8 kcal•mol -1 ) is compensated by a lower activation energy, resulting in a fast consecutive chemical conversion (ΔΔG ‡ = -2.8 kcal•mol -1 ). The reaction of 1E with LinBH272F illustrates that the desired labeling efficiency can be achieved not only by an improved ligand binding, but also by an acceleration of the chemical reaction. The reaction mechanism observed for LinBH272F can be explained by the specific architecture its access tunnel. The narrower tunnel bottleneck compromises to some extent the ligand transport, but at the same time reduces the active site solvation and makes the productive binding more probable. Moreover, the tunnel lining residues lower the initial entropy and promote the contact of the reacting atoms, possibly through specific interactions with the ligand. All these effects may have a positive effect on increasing the rate of the carbon-halogen bond cleavage (SN2), and have been described in previous studies focused on the engineering of access tunnels, for both DhaA and LinB [58][59][60][61] . However, the narrow architecture of the access tunnel in LinB makes this enzyme more sensitive to the length of the reactive linker. This is important to allow the formation of favorable interactions between the aromatic system of the probes and the residues lining the tunnel. The eightcarbon linker of 1E was the only one providing optimal length for LinB labeling since shorter linkers have not been able to achieve even remotely the efficiency of 1E.</p><p>Wider and more accessible access tunnels seem to be more universal, as evidenced by the reaction of the non-optimized DmmAH315F with 4-stilbazolium ligands. DmmA has the most open main tunnel ensuring easy access of the substrates to its active site. The reaction of DmmAH315F with all 4-stilbazolium-based ligands showed rapid binding, but also rapid chemical steps. The resulting labeling efficiency thus surpasses the commercial reaction of TMR with DhaAHT. The reaction of LinBH272F with 1E and DmmAH315F with all the 4stilbazolium ligands showed that natural variants can provide high efficiency useful for HaloTag applications without time-demanding protein optimization by directed evolution. It is also interesting that, just as DhaAHT is highly specific for the reaction with TMR, DmmAH315F showed high efficiency only in reaction with the 4-stilbazolium-based ligands, but not with TMR. Clearly, it is important to select an appropriate protein for the binding of a specific ligand.</p><p>Using several computational methods, we simulated and predicted the kinetics and thermodynamics of the two-step binding process of four representative fluorescent probe/protein systems, namely for the TMR and 1E probes with DhaAHT and LinBH272F. We found disparities in the absolute values of the calculated kinetic rates. Such differences have been reported previously 56,57 , and can be attributed to the bias intrinsic to the simulation method (adaptive sampling) or the conditions used in our MD simulations, namely the force field and the solvent model. The ligand transport in proteins is highly influenced by the and its respective bulk properties, such as diffusivity. In spite of being one of the most widely used water models in molecular simulations, the TIP3P model has a higher diffusivity than pure water. It is also known to overestimate the diffusion properties of amino acids 62 , and we can presume that that same holds for many other solvated molecules. Importantly, our results showed significant correlations with some of the experimental parameters and revealed important clues for different aspects of the molecular binding on focus here. We could qualitatively replicate the order in the k1 binding rates, with DhaAHT-TMR showing the highest value, followed by LinBH272F-1E, and partially the order of affinities. The unbinding rates, however, were less consistent with the experimental results. The simulation of the binding process revealed how the probes interacted differently with the proteins. Both used solely p1 tunnel to bind the proteins, but the preferred orientation of TMR was more compatible with the geometry of the p1 tunnel in the DhaA variants, while 1E adopted an orientation more similar to the p1 tunnel found in the native LinBH272F and DmmAH315F. This reveals a complementarity intrinsic to those two pairs that seem to explain the labeling efficiencies described above. The large number of interactions formed between 1E and residues lining the tunnel of LinBH272F support this hypothesis.</p><p>The overall chemical step was dissected into the pre-organization of the bound state to form the pre-reactive complex (NAC) and the SN2 reaction. We estimated these two partial steps from our computational approach and calculated the total activation energy of the second kinetic step (∆G ‡ 2), which can be compared with the parameter determined experimentally. We found that DhaAHT-TMR displayed the lowest overall reactivity (with the highest ∆G ‡ 2 value), which is in good agreement with the experimental data. It showed not only the worst efficiency in achieving a productive binding mode (lowest KNAC) but also presented the highest activation barrier to the SN2 reaction. Conversely, LinBH272F-1E was the most efficient system in adopting the pre-reactive conformation after the binding (highest KNAC). A low ∆G ‡ SN2 also resulted in LinBH272F-1E having a rather low activation barrier to the overall chemical step, with an estimated ∆G ‡ 2 value below that of DhaAHT-TMR by 4.4 kcal•mol -1 . Some of the discrepancies between the theoretical and experimental values were likely due to a poor sampling of the fully bound states, which may have not been sufficient to provide an accurate ensemble distribution of the pre-reactive state. However, our results provided sufficient clues to explain why, although DhaAHT-TMR presented the highest binding rate, its reactivity is very far from ideal. In contrast, although the LinBH272F-1E system had poorer binding rates, it is much more efficient on the chemical step. Overall, the binding of the 1E probe to the non-optimized LinBH272F protein revealed a reasonable binding/reactivity trade-off, which resulted in a labeling efficiency very close to that of DhaAHT-TMR. Some of the effects discussed above can be extrapolated to the DmmAH315F-1E system, which presented the best binding efficiency among all the tested pairs. Hence, we hypothesize that the binding of 1E to DmmAH315F is fast due to the combination of a sufficiently wide access tunnel and a good complementarity of its architecture with the 1E probe, which lead to a high number of favorable interactions. Secondly, strong probe-enzyme interactions can contribute to a stable and highly reactive DmmAH315F-1E complex, thus leading to a fast chemical step. The combination of a fast binging rate and a fast chemical step resulted in a system with the highest labeling efficiency.</p><!><p>Here we have demonstrated that not only the ligand accessibility is important for the binding of diverse probes to HaloTag proteins, but also the subsequent chemical step can significantly affect the ligand specificity and labeling efficiency. We have identified substantial differences in the kinetics of binding and chemical reaction between individual enzymes with different ligands. The TMR probe showed a rapid binding to DhaAHT, which was followed by a slow chemical conversion to the alkyl-enzyme complex. In contrast, the binding of the 4stilbazolium-based ligands to DhaAHT and other tag proteins was much slower than with TMR, but the chemical step was greatly improved in most cases. Interestingly, we found that the best efficiencies for the incorporation of several 4-stilbazolium-based probes (namely 1D and 1E) were achieved with the analogs of natural non-optimized dehalogenases, LinBH272F and DmmAH315F, which provided high kinetic rates for both binding and chemistry. This demonstrates that different natural proteins can be effective for the incorporation of specific probes without the need for demanding protein engineering procedures. Moreover, the 4stilbazolium-based ligands, due to a better light-up response upon binding, may provide better detection limits and thus could be preferable to the traditional probes, e.g. for simple fluorescence assays, analysis of binding interactions, or microscopy imaging.</p><p>We propose that, before conducting laborious optimization rounds by directed evolution, a rapid screening of the available natural dehalogenases could lead to the identification of potential candidates for optimal tag proteins. Thus, one could benefit from the very diverse pool of tunnel architectures already available among the known haloalkane dehalogenases. Calculation of the respective access tunnels with CAVER and molecular docking could provide good first filters for this selection. The subsequent utilization of more robust computational methods, like molecular dynamics and quantum mechanics, can help identifying the ideal enzyme-probe pairs.</p><!><p>All materials and methods are described in detail in Supplementary Information.</p>