Rational Design of a Thermostable 2′-Deoxyribosyltransferase for Nelarabine Production by Prediction of Disulfide Bond Engineering Sites

One of the major drawbacks of the industrial implementation of enzymatic processes is the low operational stability of the enzymes under tough industrial conditions. In this respect, the use of thermostable enzymes in the industry is gaining ground during the last decades. Herein, we report a structure-guided approach for the development of novel and thermostable 2′-deoxyribosyltransferases (NDTs) based on the computational design of disulfide bonds on hot spot positions. To this end, a small library of NDT variants from Lactobacillus delbrueckii (LdNDT) with introduced cysteine pairs was created. Among them, LdNDTS104C (100% retained activity) was chosen as the most thermostable variant, displaying a six- and two-fold enhanced long-term stability when stored at 55 °C (t1/255 °C ≈ 24 h) and 60 °C (t1/260 °C ≈ 4 h), respectively. Moreover, the biochemical characterization revealed that LdNDTS104C showed >60% relative activity across a broad range of temperature (30–90 °C) and pH (5–7). Finally, to study the potential application of LdNDTS104C as an industrial catalyst, the enzymatic synthesis of nelarabine was successfully carried out under different substrate conditions (1:1 and 3:1) at different reaction times. Under these experimental conditions, the production of nelarabine was increased up to 2.8-fold (72% conversion) compared with wild-type LdNDT.


Introduction
Nucleosides are essential molecules in many different biochemical pathways, such as nucleic acid synthesis, regulation and metabolism, and signal transduction. Consequently, nucleoside analogues (NAs) have attracted considerable attention as antiviral [1], anticancer [2], or antibacterial [3] active pharmaceutical ingredients (APIs) in the pharma industry.
Despite the undeniable potential of biocatalysis as a valuable tool for the industrial production of NAs, several operational drawbacks hamper the application of enzymes as industrial biocatalysts [20,21]. Among these, achieving enzyme stabilization under harsh Despite the undeniable potential of biocatalysis as a valuable tool for the industrial production of NAs, several operational drawbacks hamper the application of enzymes as industrial biocatalysts [20,21]. Among these, achieving enzyme stabilization under harsh conditions used in industry is the most challenging from a practical standpoint. Modern biocatalysis demands new biocatalysts active and stable under aggressive reaction conditions, e.g., extreme pH values, high temperatures, or the presence of organic solvents and hazardous reagents. In this respect, the use of different strategies to get robust biocatalysts is a habitual practice in the pharma industry. These include the immobilization of the enzymes onto different matrixes [22,23], the use of enzymes from extremophilic organisms [24,25], or the tailor-made engineering of more stable enzymes [26,27].
The rational engineering of more stable variants is strongly dependent on structurebased computational design, which can identify hotspots on the protein sequence for mutagenesis to reach the optimization of flexible regions, the rigidification of the structure, or the introduction of specific stabilizing interactions and metal-chelating sites [28,29]. In this context, several computational tools are commonly used to predict stabilizing mutations, either based on sequence (consensus mutagenesis, proline substitutions, and disulfide bond formation) [30][31][32] or on structure knowledge (molecular dynamics, B-FITTER, constrain network analysis, and FoldX) [33][34][35][36].
Since LdNDT displays a homohexameric assembly, similar to other reported bacterial homologs from L. leichmannii, L. helveticus, or Bacillus psychrosaccharolyticus [11,37,38], this oligomeric character opened up new design possibilities for the incorporation of intersubunit disulfide bridges by following a structure-guided approach. In this context, based on the in silico predictions, four variants (LdNDTV63C, LdNDTV93C, LdNDTS104C, and the double mutant LdNDTV93C/S104C) were produced and characterized. Among them, LdNDTS104C displayed similar activity to that shown by wild-type LdNDT, but also an enhanced thermostability (up to six-and two-fold t1/2 enhancement when stored at 55 and 60 °C, respectively), which demonstrated the effect of this mutation on thermal stability. Finally, LdNDTS104C was successfully applied as a biocatalyst for the synthesis of nelarabine ( Figure  1), improving the conversion displayed by the wild-type counterpart.

Computational Design of LdNDT Variants
The use of NDTs as biocatalysts for the synthesis of NAs has been extensively described [6,9,[37][38][39][40][41][42][43][44][45]. However, even though NDT-mediated transglycosylation emerges as a versatile tool for industrial manufacturing of NAs, the long reaction times (24-48 h) required for the synthesis of some NAs (2 and 3 modified nucleosides) usually compromise the enzyme lifespan. Recently, we have described the presence of inter-and intra-subunit disulfide bridges (DSB) in NDT from Desulfotalea psychrophila (DpNDT, PDB ID 7O62) [41], which contributes to stabilizing the protein folding. Inspired by these findings, we envisioned that the introduction of DSBs into Lactobacillus delbrueckii NDT structure may improve thermal stability.
To this end, based on the high sequence identity (98%) with hexameric NDT from L. leichmannii (PDB ID 1F8Y) [46], a 3D model was constructed [47]. Based on the results of calculations with the SSBONDPredict [48] several inter-and intra-subunit DSBs were proposed. Moreover, the effect of these mutations on protein structure was analyzed using HotSpot Wizard [30]. Finally, four variants: LdNDT V63C , LdNDT V93C , LdNDT S104C , and the double mutant LdNDT V93C/S104C were proposed as tentative candidates for in vitro experiments ( Figure 2 and Table 1). 3D models for each LdNDT variant were constructed, and MD simulations were run. Distances and angles for introduced Cys residues were consistent with disulfide bond formation [49,50]. which contributes to stabilizing the protein folding. Inspired by these findings, we envisioned that the introduction of DSBs into Lactobacillus delbrueckii NDT structure may improve thermal stability.
To this end, based on the high sequence identity (98%) with hexameric NDT from L. leichmannii (PDB ID 1F8Y) [46], a 3D model was constructed [47]. Based on the results of calculations with the SSBONDPredict [48] several inter-and intra-subunit DSBs were proposed. Moreover, the effect of these mutations on protein structure was analyzed using HotSpot Wizard [30]. Finally, four variants: LdNDTV63C, LdNDTV93C, LdNDTS104C, and the double mutant LdNDTV93C/S104C were proposed as tentative candidates for in vitro experiments ( Figure 2 and Table 1). 3D models for each LdNDT variant were constructed, and MD simulations were run. Distances and angles for introduced Cys residues were consistent with disulfide bond formation [49,50].    The location of the disulfide bridges within the hexameric assembly of LdNDT is indicated in the upper part of Figure 2, which shows LdNDT along its three-fold molecular symmetry axis. It should be remarked here that introduction of just one cysteine in the amino acid sequence of the enzyme results in the introduction of three disulfide bridges in the assembly. Therefore, the global effects observed in thermal stability of the LdNDT variants (see below) would result from a joint combination of the internal symmetry of the enzyme and the introduction of the Cys residue, namely, the effect of the mutation is amplified by the molecular symmetry. Conversely, the lower panels of Figure 2 show close-up views of the amino acid changes incorporated: V63C, V93C, and S104C (from left to right). The distances between sulfur atoms are within the range observed for disulfides, as well as those for the Cα-Cα distances (5.6 Å for V63C, 4.5 Å for V93C, and 6.5 Å for S104C; typical range: 3.0-7.5 Å).

Production of LdNDT Variants
Taking into account these previous features, the engineered LdNDT variants were produced and purified (see Section 3). SDS-PAGE analysis of the purified variants revealed only one protein band with an apparent molecular mass of 18 kDa ( Figure S1). Once the LdNDT variants were purified, we assayed the N-2 -deoxyribosyltransferase activity in the wild-type and engineered variants under standard assay conditions (Table 1 and Figure S2).
As shown in Table 1, LdNDT V63C and LdNDT V93C/S104C variants exhibited a drastic loss of activity compared to the wild-type counterpart, whereas LdNDT V93C did not show a remarkable loss of activity under standard assay conditions. Interestingly, a similar activity was observed for the LdNDT S104C variant.
In an attempt to explain the activity difference between LdNDT variants, distances from OE1 of catalytic Glu98 to C1 of dUrd (d1), as well as distances between H from the C-terminal (HTX) of catalytic Tyr157# and O4 of dUrd (d2), were measured ( Figure 3). The evaluation of these distances through MD simulations may explain the distortion of the active site. As shown in Figure 3, the presence of cystine mutations leads to significant changes in the d1 distance ( Figure 3A, average d1 distances) compared to the wild-type counterpart. However, no significant differences are observed in the d2 distance, and more importantly, the average values for all LdNDT variants fall below 2.5 Å ( Figure 3B, average d2 distances). To that effect, as previously reported, the productive complex is formed when d1 and d2 distances are less than 3.45 Å and 2.5 Å, respectively [51]. Interestingly, the MD simulation of LdNDT S104C reveals an average d1 distance located under 3.45 Å, which is consistent with the relative activity values reported in Table 1. For the rest of the mutants, the average d1 distance stands above 3.45 Å, which could be related to the observed lower activity values. More interestingly, the distribution curve for LdNDT S104C d1 distances falls within that shown for productive complex for LdNDT (d1 < 3.45 Å) [51], which is not observed for the other variants. All in all, the computational analysis of enzymeligand complexes revealed that the presence of a cystine could distort the enzyme-ligand interactions and therefore interfere with substrate binding and catalysis.

Thermal Stability
The effect of temperature on the stability of LdNDT, LdNDTV93C, and LdNDTS104C was assayed. As shown in Figure 4, significant differences in thermal stability were observed for wild-type LdNDT and mutant variants when incubated at 55 and 60 °C ( Figure 4A,B).  Although LdNDT losses~50% of the initial activity in the first 4 h of incubation, almost no further drop in enzyme activity is observed until reaching 24 h of incubation at 55 • C (23% relative activity) ( Figure 4A). LdNDT V93C displayed a similar tendency, losing 50% of the initial activity in the first 4 h of incubation (t 1/2 55 • C ≈ 4 h). However, the activity drops significantly in the following hours, observing 11% relative activity after 15 h of incubation (data not shown). In contrast, LdNDT S104C displays up to six-fold thermal stability enhancement (t 1/2 55 • C ≈ 24 h) compared to the wild-type enzyme ( Figure 4A). As a result of this first test, LdNDT S104C was selected for further thermostability experiments (60 • C) ( Figure 4B). LdNDT suffered a remarkable loss of activity at shorter incubation periods when stored at 60 • C (t 1/2 60 • C ≈ 2 h), whereas a higher half-time was observed for LdNDT S104C (t 1/2 60 • C ≈ 4 h, two-fold improvement). Additionally, LdNDT and LdNDT S104C display similar long-term stability when stored at 55 • C in the presence of the reducing agent 1 mM DTT ( Figure S3), probably due to the reduction of the disulfide bonds in LdNDT S104C .
To deepen into the molecular basis of experimental results, MD simulations at different temperatures were carried out. Throughout these simulations, variations of RMSD (root mean square deviation of atomic positions); RMSF (root mean square fluctuation of protein backbone); and B-factor (so called the Debye-Waller factor, temperature factor, or atomic displacement parameter) values were evaluated as a measurement of the effect of mutation on the overall conformation and flexibility shown by the enzyme [52][53][54]. For that matter, higher variations in these values are associated with superior conformational changes and increased flexibility in the protein structure. Since disulfide bond formation leads to a more rigid structure, nonsignificant variations in the RMSD, RMSF, and B-factor values are expected when the temperature rises.
As shown in Figure 5A, similar RMSD values are observed for both LdNDT and LdNDT S104C through MD simulations in the 25-50 • C temperature range. When simulations were performed at 95 • C, the RMSD values significantly increase for LdNDT after 50 ns, while a similar tendency is observed for LdNDT S104C after 150 ns. Furthermore, to explain these RMSD differences, the RMSF variations were also evaluated for both enzymes at different temperatures ( Figure 5B). Attending to the obtained results, a significant decrease in RMSF values was observed in the vicinity of the Cys104 position for LdNDT S104C , probably related to a structural rigidification by the presence of a disulfide bond. In addition, MD simulations of LdNDT S104C show lower RMSF values for catalytic residues Glu98 and Tyr157# at higher temperatures (compared to wild-type LdNDT), which is consistent with the long-term operational stability observed for LdNDT S104C .
Additionally, the overall flexibility of enzymes was analyzed by the measurement of average B-factor values through independent MD simulations at different temperatures (25-95 • C). The computational results revealed a 1.4-fold increment of the B-factor average value for LdNDT S104C and a 2.5-fold increase for LdNDT ( Figure 5C) when the temperature rises from 25 • C to 50 • C. Moreover, broader differences were observed in the 25-95 • C range, displaying 2.0-and 3.3-fold increments for LdNDT S104C and LdNDT, respectively. Finally, attending to their B-factor cartoon representation, significant differences are detected in the flexibility of LdNDT and LdNDT S104C ( Figure 5D). Finally, attending to their B-factor cartoon representation, significant differences are detected in the flexibility of LdNDT and LdNDTS104C ( Figure 5D).

Biochemical Characterization
Taking into account the abovementioned results, LdNDTS104C emerges as the most valuable candidate from an industrial perspective. Therefore, the next step was to establish the optimal conditions to reach maximum activity. To this end, the effect of temperature and pH on LdNDTS104C activity was studied ( Figure 4C,D) and compared to those results previously reported for LdNDT [11]. When compared to its wild-type counterpart, LdNDTS104C displayed high activity values (more than 60%) across a broader temperature

Biochemical Characterization
Taking into account the abovementioned results, LdNDT S104C emerges as the most valuable candidate from an industrial perspective. Therefore, the next step was to establish the optimal conditions to reach maximum activity. To this end, the effect of temperature and pH on LdNDT S104C activity was studied ( Figure 4C,D) and compared to those results previously reported for LdNDT [11]. When compared to its wild-type counterpart, LdNDT S104C displayed high activity values (more than 60%) across a broader temperature range (from 30 • C to 90 • C) ( Figure 4C). It is particularly interesting to note that a negligible loss of activity was observed at 80 • C and 90 • C for the mutant enzyme, while just a~20% relative activity was reported for wild-type LdNDT at the same temperatures [11]. Thus, besides the enhanced thermal stability, the proposed mutation could also improve the enzyme activity at extremely high temperatures. As reported for LdNDT, the maximum activity for LdNDT S104C was observed at 60 • C. Regarding the pH dependence, LdNDT S104C displays a similar tendency to that shown by LdNDT [11] (Figure 4D), with a maximum in the pH range of 6-7.

Enzymatic Synthesis of Nelarabine
Nelarabine, a water-soluble prodrug of ara G licensed by SmithKline Beecham in 2005, received approval for the treatment of T-cell acute lymphoblastic leukemia [1]. Chemical methods for nelarabine synthesis involved condensation between an arabinosyl derivative and 6-methoxy guanosine, which usually leads to anomeric mixtures and, therefore, to a significant decrease in the efficiency of the process [55]. In addition, these strategies usually employ expensive and environmentally unfriendly chemical reagents and organic solvents. Therefore, the underlying need for waste reduction and high (enantio)-selectivities fostered the introduction of biocatalysis as a sustainable technology for the industrial synthesis of nelarabine.
In contrast, enzymatic transglycosylation emerges as a cost-effective and eco-friendly alternative to chemical methods. In recent years, nucleoside phosphorylases (NPs) have been widely employed as powerful catalysts for nelarabine synthesis [6,56]. However, the high cost of α-D-arabinofuranose 1-phosphate (sugar donor) is a serious constraint for their industrial application. In this context, different alternatives are now under study by pharma industries. Among them, the NDT-mediated synthesis of nelarabine from arabinosyl nucleosides and 6-methoxyguanine (6-MeOGua) could offer an eco-friendly but also cheaper process alternative [10,39,[57][58][59][60].
Herein, the authors propose an NDT-mediated one-pot synthesis of nelarabine from ara U and 6-MeOGua. Since long reaction times are required for the NDT-mediated synthesis of arabinosyl nucleosides [10,39,[57][58][59][60], it is expected that a more stable NDT (LdNDT S104C ) would allow an increase in conversion rates.
As shown in Figure 6A, the variant LdNDT S104C displayed around a 3.7-fold higher production yield (compared to wild-type LdNDT) when using a 1:1 nucleoside/base ratio. The observed improvement could be due to the enhanced thermostability displayed by LdNDT S104C , which allows the enzyme to be more active during the long reaction periods employed for the production of nelarabine. Furthermore, to improve the nelarabine production yields, a 3:1 substrate ratio was employed under the same reaction conditions ( Figure 6B). Interestingly, the conversion rates were significantly improved for both LdNDT (from 10 to 26%) and LdNDT S104C (from 39 to 72%). Compared to previous results (1:1 ratio experiments), the differences in the conversion rates for both enzymes were diminished; however, a 2.8-fold increment was still observed for LdNDT S104C . Once again, when compared to LdNDT, the observed improvement could be due to the enhanced thermostability displayed by LdNDT S104C , which allows the enzyme to be more active during the long reaction periods. Finally, as a proof of concept of the potential of LdNDTS104C as an industrial biocatalyst, we performed the synthesis of nelarabine at high concentrations of ara U (15-30 mM) and 6-MeOGua (5-10 mM) over different reaction times (ranging from 12 to 72 h) at 50 °C since thermal stability of the enzyme was guaranteed (Figure 7). Finally, as a proof of concept of the potential of LdNDT S104C as an industrial biocatalyst, we performed the synthesis of nelarabine at high concentrations of ara U (15-30 mM) and 6-MeOGua (5-10 mM) over different reaction times (ranging from 12 to 72 h) at 50 • C since thermal stability of the enzyme was guaranteed (Figure 7).  1). (B) ara U:6-MeOGua (3:1). All determinations were carried out in triplicate, and the standard error of the mean is calculated using the standard deviation.
Finally, as a proof of concept of the potential of LdNDTS104C as an industrial biocatalyst, we performed the synthesis of nelarabine at high concentrations of ara U (15-30 mM) and 6-MeOGua (5-10 mM) over different reaction times (ranging from 12 to 72 h) at 50 °C since thermal stability of the enzyme was guaranteed (Figure 7).

Figure 7.
Time-course of the enzymatic production of nelarabine catalyzed by LdNDTS104C at different substrate concentrations: ara U (3-30 mM) and 6-MeOGua (1-10 mM). All determinations were carried out in triplicate, and the standard error of the mean is calculated using the standard deviation.
Based on the previous literature for wild-type LdNDT [11], high promiscuity on nucleobase recognition is expected. Together with the enhanced thermostability displayed by LdNDTS104C, suggest that other different APIs could be synthesized through this methodology, such as vidarabine (adenine arabinoside, ara A), fludarabine (2-fluoroadenine arabinoside, ara FA), or clofarabine (2-chloradenine-2′-fluoroarabinoside), among others. However, despite this enzymatic approach vastly surpassing chemical procedures in terms of efficiency (one-pot regio-and enantioselective reactions) or sustainability (green conditions and aqueous buffers), from an industrial perspective, several operational constraints need to be addressed. One of the major drawbacks is the high cost of production and purification of recombinant enzymes, which makes necessary the capability of reusing the catalyst on successive reactions to reach a cost-effective process. In this regard, enzyme immobilization offers valuable strategies to overcome this operational handicap but also the possibility to integrate the immobilized catalyst in continuous flow reactors. Figure 7. Time-course of the enzymatic production of nelarabine catalyzed by LdNDT S104C at different substrate concentrations: ara U (3-30 mM) and 6-MeOGua (1-10 mM). All determinations were carried out in triplicate, and the standard error of the mean is calculated using the standard deviation.
Based on the previous literature for wild-type LdNDT [11], high promiscuity on nucleobase recognition is expected. Together with the enhanced thermostability displayed by LdNDT S104C , suggest that other different APIs could be synthesized through this methodology, such as vidarabine (adenine arabinoside, ara A), fludarabine (2-fluoroadenine arabinoside, ara FA), or clofarabine (2-chloradenine-2 -fluoroarabinoside), among others. However, despite this enzymatic approach vastly surpassing chemical procedures in terms of efficiency (one-pot regio-and enantioselective reactions) or sustainability (green conditions and aqueous buffers), from an industrial perspective, several operational constraints need to be addressed. One of the major drawbacks is the high cost of production and purification of recombinant enzymes, which makes necessary the capability of reusing the catalyst on successive reactions to reach a cost-effective process. In this regard, enzyme immobilization offers valuable strategies to overcome this operational handicap but also the possibility to integrate the immobilized catalyst in continuous flow reactors.

Materials
All the substrates and reagents for microbiology experiments were from Condalabs (Madrid, Spain). Moreover, HPLC solvents and buffers were purchased from Sigma-Aldrich (Madrid, Spain). All other compounds involved in enzyme reactions were from Carbosynth Ltd. (Compton, UK).

Gene Expression and Protein Purification
All the DNA constructions were purchased from Genscript (Piscataway, NJ, United States) as NdeI-EcoRI fragments subcloned into the expression vector pET22b(+). The recombinant plasmids were used to transform E. coli BL21 (DE3) cells under conditions described elsewhere [11]. The overexpression of wild-type Lactobacillus delbrueckii NDT (UniProtKB Q1GC13) and variants and the cell extract preparation for protein purification were performed following the previous protocol reported by Acosta et al. [11].
LdNDT and mutant variants were purified by ammonium sulfate precipitation. To this end, the supernatant obtained after centrifugation of the lysed cells was brought to 40% (NH 4 ) 2 SO 4 saturation for the selective precipitation of the enzymes. The prepared solutions were incubated for 1 h at 4 • C and then centrifugated at 16,500× g. The collected pellets were resuspended in 10 mM sodium phosphate buffer, pH 7.0, heated at 67 • C for 7 min, and further incubated at 4 • C. After a centrifugation step (15,800× g), the supernatant was collected. Then, the protein mixture was purified by size-exclusion chromatography (HiLoad 16/60 Superdex 200, GE Healthcare, Madrid, Spain) (50 mM sodium phosphate, pH 7.0). The fractions comprising pure enzymes were identified by an SDS-PAGE and further confirmed by an activity assay. Finally, the protein concentration was spectrophotometrically quantified as described before [11].

N-2 -Deoxyribosyltransferase Activity Assay
The standard assay was carried out using 0.3 µg of the pure enzyme with 10 mM substrates (2 -deoxyuridine, dUrd, and 10 mM adenine, Ade) in 50 mM MES buffer, pH 6.5 (reaction volume 40 µL) at 50 • C for 5 min (300 rpm). Then, the reaction was stopped, and the samples were half-diluted with water for HPLC analysis. One unit of activity (U) was established as the milligrams of the enzyme that produces 1 µmol of product per min (IU).

Thermal Inactivation
LdNDT and mutants were stored at different temperature conditions (4 and −80 • C) for around 200 days. Moreover, to determine the effect of temperature on the operational stability of LdNDT, LdNDT V93C , and LdNDT S104C , 0.3 µg of purified enzyme were incubated in 10 mM sodium phosphate buffer, pH 7.0, for 30 h at different temperatures (55-60 • C). During these periods, the enzymatic activity was systematically assayed. Moreover, LdNDT and LdNDT S104C were also incubated in the presence of 1 mM DTT (dithiothreitol) for 24 h.

Enzymatic Synthesis of Nelarabine
To evaluate the potential use of LdNDT S104C as a catalyst, the enzymatic production of nelarabine was carried out under different reaction conditions (reaction time and substrate ratio). To this end, 10 µg of LdNDT or LdNDT S104C were incubated with 1 mM arabinosyl uracil (ara U) and 6-methoxyguanine (6-MeOGua) in 50 mM MES buffer pH 6.5, at 50 • C in a final volume of 40 µL, with 300 rpm orbital shaking at different reaction times (12-72 h). Thereafter, the reactions were stopped following the standard assay described above. Further experiments were carried out at a 3:1 nucleoside/base substrate ratio (ara U/6-MeOGua; 3 mM/1 mM) to improve the conversion yield. The enzymatic activity was quantified by HPLC as described above.
Finally, to demonstrate the industrial applicability of LdNDT S104C , the enzymatic synthesis of nelarabine from high concentrations of ara U (15-30 mM) and 6-MeoGua (5-10 mM) was performed.

Molecular Modeling
The hexameric model for apo LdNDT was constructed based on the structure of type II 2 -deoxyribosyltransferase from L. leichmannii complexed with 5-methyl-2 -deoxypseudouridine (PDB ID 1F8Y) (98% sequence identity). Titratable residues were properly protonated at pH 6.5 using the H ++ 3.0 webserver [61]. Then, dUrd was manually docked into the active site through structural best-fit superposition onto the former ligand (5-methyl-2deoxypseudouridine). Previously, the Antechamber (AmberTools16) [62] program was utilized for the optimization of the ground state geometry of the dUrd molecule, as well as for computing its electrostatic potential using quantum mechanical HF/6-31G* wavefunction fitted to the atoms as RESP charges.
The consistency of our 3D model was validated through MD simulation using the leaprcff14SB force field. The MD simulation was run in the Single-Precision-Fixed-Precision (SPFP) mode implemented in the pmemd.cuda module of Amber16. The LdNDT:dUrd complex was simulated as a hexamer, and the complex was embedded in a cubic box of TIP3P water molecules [63], including 27 Na + ions, to achieve charge neutrality. As described before for this enzyme [47], the total energy of the system was minimized in three consecutive steps (3 × 20,000 cycles). In the first 5000 cycles, a steepest descendent method was employed, while the following minimization cycles were calculated by the conjugate gradient method. After a 200 ps heating (100-300 K) restricting the movement of all solute atoms [64], the simulation was run with a fixed volume (NVT ensemble) as previously described [47].
Once the consistency of the model was confirmed, the system was employed for the construction LdNDT variants to improve the thermostability. Those mutations were predicted based on both structural and sequence information [41,51,65]. Mutations leading to disulfide bond formation were predicted based on the structural knowledge using the Neural Network Model AI algorithm implemented in SSBONDPredict software [48]. Those mutations were further confirmed by the sequence analysis through the HotSpot Wizard server [30]. All in all, three potential mutable positions were identified and evaluated by MD simulations following the abovementioned protocol. In this respect, selected amino acids were mutated to cysteine residues with PyMOL [66] for the construction of four LdNDT variants: LdNDT V63C , LdNDT V93C , LdNDT S104C , and LdNDT V93C/S104C . Cystines were formed by tleap software available in AmberTools16, and the suitableness of the models was confirmed by MD simulations. Results for the MD simulations were analyzed with the cpptraj module available in AmberTools16 [62].

QM/MM MD Simulations
MD simulations for all mutant candidates were run as described before to refine the LdNDT:dUrd Michaelis complex for the first half-reaction. The distances between the reactant OE1 of Glu98 and the C1' anomeric carbon of dUrd, as well as H from the C-terminal of Tyr157# and O4 of dUrd, were monitored. The best snapshot was selected as previously described by Del Arco and co-workers [65] and further simulated for 50 ps through QM/MM MD simulation without any restraint. For the first half-reaction, the ligand (dUrd) and catalytic residues, Glu98 and Tyr157#, were included in the QM region. As previously described [65], the QM region was treated at the self-consistent-charge density-functional tight-binding method (SCC-DFTB) [67] using the DFTB3 formulation, whereas the rest of the system was treated with MM MD, as described above. The results for QM/MM MD simulations were analyzed with the cpptraj module available in AmberTools16 [62].

MD Simulations for Thermal Stability Predictions
Once the best mutant candidate was selected (LdNDT S104C ), MD simulations at different temperatures were run to compare the thermal stabilities of LdNDT and LdNDT S104C . More specifically, these simulations give an insight into the flexibility of the enzymes at high temperatures, which can be related to thermal stability [52]. The flexibility was evaluated through the measurements of the RMSD, RMSF, and B-factor values [52][53][54] through the MD simulations. In this sense, these simulations were made following the above-described protocols but at different temperatures ranging from 300 K to 400 K. the RMSD, RMSF, and B-factor values from the MD simulations were analyzed with the cpptraj module available in AmberTools16 [62]. The MD simulations were carried out in triplicate.

Conclusions
Herein, we report the structure-guided design of a novel thermostable NDT based on the DSB design in LdNDT. After a first screening focused on enzyme activity, LdNDT V93C and LdNDT S104C were selected as candidates for a second screening strictly focused on operational stability. Since the experimental results demonstrated LdNDT S104C as the most stable catalyst, this variant was successfully employed as a catalyst in the enzyme production of nucleoside analogs at long reaction times. Interestingly, a remarkable effect on nelarabine production (2.8-fold increment, 72% conversion) was observed when compared with that shown for LdNDT, which reinforces our initial hypothesis.

Conflicts of Interest:
The authors declare no conflict of interest.