Engineered Stable 5-Hydroxymethylfurfural Oxidase (HMFO) from 8BxHMFO Variant of Methylovorus sp. MP688 through B-Factor Analysis

What is known as Furan-2,5-dicarboxylic acid (FDCA) is an attractive compound since it has similar properties to terephthalic acid. Further, 5-hydroxymethylfurfural oxidase (HMFO) is an enzyme, which could convert HMF to FDCA directly. Most wild types of HMFO have low activity on the oxidation of HMF to FDCA. The variant of 8BxHFMO from Methylovorus sp. MP688 was the only reported enzyme that was able to perform FDCA production. However, the stabilization of 8BxHMFO is still not that satisfactory, and further improvement is necessary for the industrial application of the enzyme. In this work, stability-enhanced HMFO from 8BxHFMO was engineered through employing B-factor analysis. The mutation libraries were created based on the NNK degeneracy of residues with the top ten highest B-factor value, and two of the effective mutants were screened out through the high throughput selection with the horseradish peroxidase (HRP)-Tyr assay. The mutants Q319K and N44G show a significantly increased yield of FDCA in the reaction temperature range of 30 to 40 °C. The mutant Q319K shows the best performance at 35 °C with a FDCA yield of 98% (the original 8BxHMFO was only 85%), and a half-life exceeding 72 h. Moreover, molecular dynamic simulation indicates that more hydrogen bonds are formed in the mutants, which improves the stability of the protein structure. The method could enhance the design of more stable biocatalysts; and provides potential for the further optimization and utilization of HMFO in biotechnological processes.


Introduction
Since the fossil resource depletion and degradation of the environment, renewable materials and green technologies attract more attention [1]. Therefore, using biomass as raw material to produce chemicals and biomaterials has become the future development trend.
Polyethylene terephthalate (PET) is one of the widely utilized polyesters in our life. However, as the basic material of its synthesis, terephthalic acid (TPA) is mainly produced from the nonrenewable petroleum industry [2]. It was found that furan-2,5dicarboxylic acid (FDCA) has similar properties to TPA, which could be used as an ideal substitute for TPA to produce PET [3,4]. Since FDCA could be produced through 5hydroxymethylfurfural (HMF), and the latter could be converted from lignocellulose, FDCA was listed as a bioderived C6 platform compound by the International Energy Agency (IEA) in 2012 [5]. Therefore, based on the topic of the green and low-carbon circular economics [6], the development of FDCA production technology has gained a lot of interest in the past decade.
Among all kinds of production methods of FDCA, the oxidation of HMF to produce FDCA was considered as the most promising application prospect since HMF is renewable, plentiful, and cheap [7]. At present, the oxidation process is mostly catalyzed by homogeneous or supported metal catalysts [8]. Nevertheless, disadvantages such as high costs, critical reaction conditions, low yields, uncontrolled side reactions, and environmentally unfriendly requirement of toxic organic solvents, make the method less applicable [9]. FDCA could also be synthesized from HMF through 5-hydroxymethylfurfural oxidase (HMFO) [10]. However, the conversion rate of FDCA from HMF through HMFO catalysis is not satisfactory. Since most of the HMFO can only oxidize HMF to 5-formyl-2furancarboxylic acid (FFA), other enzymes are needed for a complete oxidation of HMF [11]. Among all HMFOs, a specific one from Methylovorus sp. However, MP688 is an exception which can independently convert HMF into FDCA, although the reaction was inefficient [12]. Evolutionary works were also carried out on this special HMFO and it was found that a variant (8BxHMFO) including eight mutations: I73V, H74Y, G356H, V367R, T414K, A419Y, A435E, and W466F, could transfer over 90% of HMF to FDCA in 24 h with a substrate concentration of 5 mM, and enzyme loading of 2 µM. This work promoted an alternative way to assess numerous previous investigations on the utilization of 8BxHMFO with enzymatic catalysis for the FDCA production [13]. However, the stabilization of 8BxHMFO is still not that satisfactory, the enzyme only displayed good performance at low temperature, further improvement was necessary.
In this work, we applied homology modeling and molecular dynamic (MD) simulations to identity the potential unstable residues, and further mutation libraries were build based on the NNK degeneracy (N: adenine/cytosine/guanine/thymine; K: guanine/thymine) to encode all 20 proteinogenic amino acids [14,15]. The beneficial mutants were selected based on a quick and efficient high-throughput screening method. The performance of the mutants, and the improvement mechanism were also investigated.

Potential Mutation Site Prediction
Normally, the protein stabilization is affected by the strength of interactions between residues that form microstructures of protein, which include intramolecular hydrogen bonds, hydrophobic forces, salt bridges, and disulfide bonds [16]. Highly flexible residues with less interactions with other residues could trigger protein unfolding. Therefore, the intensification of these kinds of flexible residues could increase the stabilization of the whole protein [17]. The B-factor (or B value) of residues represents their flexibility, and indicates the conformation of residue is vulnerable. Thus, the residues with higher B-factor value might be the un-stabilizing point of the protein [18]. We targeted the residues with the top ten highest B-factor values in 8BxHMFO ( Figure 1a) and analyzed their distribution in the 3D structure of the protein (Figure 1b). Most of those residues are distributed on the surface of the protein and form as a loop structure, which indicates that the microstructure might be unstable, and the mutation of them may have less effect on the activity of the enzyme. Thus, those ten residues were targeted for randomization by using standard saturation mutagenesis [19].

Variant Screening
In this study, NNK degeneracy to encode all 20 proteinogenic amino acids was used for saturation mutagenesis of each candidate. Since at least 96 clones are required statistically for 95% coverage of residue randomization, a high through-put screening method is necessary. Normally, mutation screening of oxidase could be achieved by coupling the enzymatic oxidization product of H 2 O 2 to a horseradish peroxidase (HRP) and 2,2 -azinobis (3-ethyl-benzthiazoline)-6-sulfonic acid (ABTS) reaction. In the coupling reaction, ABTS forms a product colored in green, which could be measured spectrophotometrically at 405 nm [20]. However, the instability and low sensitivity makes the HRP-ABTS assay unsuitable for this investigation. Tyrosine (Tyr) could also be oxidized with H 2 O 2 catalyzed by HRP (Scheme 1) [21], and a fluorescent product is formed with the excitation and emission wavelengths of 330 and 415 nm, respectively (Figure 2a). Compared to the ABTS oxidation product, this fluorescent product is more stable, and the linearity range of the assay confirmed its proportionality with the concentration of H 2 O 2 varying from 0.02 to 10 µM (Figure 2b). Therefore, the HRP-Tyr assay was chosen as a high through-put screening method.
Catalysts 2021, 11, x FOR PEER REVIEW3 of 13 (a) (b) Figure 1. Mutation residues selection based on the B-factor. (a) the information of the B-factor for residues of 8BxHMFO, and the potential mutational sites, the red squares are selected residues. (b) the distribution of the potential mutational residues (green), residues of active centers were colored in cyan [19], FAD was colored in yellow, the orange zone is the active center, and the blue zone is the substrate channel.

Variant Screening
In this study, NNK degeneracy to encode all 20 proteinogenic amino acids was used for saturation mutagenesis of each candidate. Since at least 96 clones are required statistically for 95% coverage of residue randomization, a high through-put screening method is necessary. Normally, mutation screening of oxidase could be achieved by coupling the enzymatic oxidization product of H2O2 to a horseradish peroxidase (HRP) and 2,2′-azinobis (3-ethylbenzthiazoline)-6-sulfonic acid (ABTS) reaction. In the coupling reaction, ABTS forms a product colored in green, which could be measured spectrophotometrically at 405 nm [20]. However, the instability and low sensitivity makes the HRP-ABTS assay unsuitable for this investigation. Tyrosine (Tyr) could also be oxidized with H2O2 catalyzed by HRP (Scheme 1) [21], and a fluorescent product is formed with the excitation and emission wavelengths of 330 and 415 nm, respectively ( Figure 2a). Compared to the ABTS oxidation product, this fluorescent product is more stable, and the linearity range of the assay confirmed its proportionality with the concentration of H2O2 varying from 0.02 to 10 µM (Figure 2b). Therefore, the HRP-Tyr assay was chosen as a high through-put screening method. Figure 1. Mutation residues selection based on the B-factor. (a) the information of the B-factor for residues of 8BxHMFO, and the potential mutational sites, the red squares are selected residues. (b) the distribution of the potential mutational residues (green), residues of active centers were colored in cyan [19], FAD was colored in yellow, the orange zone is the active center, and the blue zone is the substrate channel. utation residues selection based on the B-factor. (a) the information of the B-factor for residues of 8BxHMFO, and l mutational sites, the red squares are selected residues. (b) the distribution of the potential mutational residues idues of active centers were colored in cyan [19], FAD was colored in yellow, the orange zone is the active center, e zone is the substrate channel.

Variant Screening
In this study, NNK degeneracy to encode all 20 proteinogenic amino acids was used for saturation mutagenesis of each candidate. Since at least 96 clones are required statistically for 95% coverage of residue randomization, a high through-put screening method is necessary. Normally, mutation screening of oxidase could be achieved by coupling the enzymatic oxidization product of H2O2 to a horseradish peroxidase (HRP) and 2,2′-azinobis (3-ethylbenzthiazoline)-6-sulfonic acid (ABTS) reaction. In the coupling reaction, ABTS forms a product colored in green, which could be measured spectrophotometrically at 405 nm [20]. However, the instability and low sensitivity makes the HRP-ABTS assay unsuitable for this investigation. Tyrosine (Tyr) could also be oxidized with H2O2 catalyzed by HRP (Scheme 1) [21], and a fluorescent product is formed with the excitation and emission wavelengths of 330 and 415 nm, respectively ( Figure 2a). Compared to the ABTS oxidation product, this fluorescent product is more stable, and the linearity range of the assay confirmed its proportionality with the concentration of H2O2 varying from 0.02 to 10 µM ( Figure 2b). Therefore, the HRP-Tyr assay was chosen as a high through-put screening method. Scheme 1. The mechanism of a coupled HRP-Tyr spectrofluorimetric assay; (Step 1): the oxidation of HMF to FDCA [12,13]; (Step 2): the polymerization of Tyr with H 2 O 2 catalyzed by HRP, and the product is a fluorescent compound [21].
In order to screen all the stable enhanced variants, 150 clones (statistically for 99% coverage) were screened from the mutation library of each candidate. Finally, five mutations (Q319K, N44G, N44K, K328G, and R252L) with higher fluorescence intensity than 8BxHMFO were selected ( Figure 3a) to further investigate the effect of mutagenesis on the FDCA conversion at different temperatures. However, only Q319K and N44G show better conversion rates than the original 8BxHMFO ( Figure 3b). Variant Q319K displays the best performance at temperature in excess of 30 • C. The B-factor analysis reveals that the Q319K and N44G show high values besides Arg 252 (Figure 1a), which indicates that they are flexible, and have a major impact on the stability of the protein [22]. Interestingly, the mutagenesis of residue Arg 252 (with the highest B-factor value) shows less effect. As illustrated in Figure 1b, Arg 252 is located on a separated loop which has less interaction with other secondary structures, and the residue is far from the active center and substrate channel. Thus, the mutation of the residue might have less effect on the stabilization of the protein. The aim of this work is to improve the stability of HMFO for FDCA production, thus mutants of Q319K and N44G were selected, and their combinatorial mutations were utilized for further catalytic performance investigation. In order to screen all the stable enhanced variants, 150 clones (statistically for 9 coverage) were screened from the mutation library of each candidate. Finally, five mutatio (Q319K, N44G, N44K, K328G, and R252L) with higher fluorescence intensity than 8BxHM were selected (Figure 3a) to further investigate the effect of mutagenesis on the FDC conversion at different temperatures. However, only Q319K and N44G show bet conversion rates than the original 8BxHMFO ( Figure 3b). Variant Q319K displays the b performance at temperature in excess of 30 °C. The B-factor analysis reveals that the Q31 and N44G show high values besides Arg 252 (Figure 1a), which indicates that they are flexib and have a major impact on the stability of the protein [22]. Interestingly, the mutagenesis residue Arg 252 (with the highest B-factor value) shows less effect. As illustrated in Figure  Arg 252 is located on a separated loop which has less interaction with other seconda structures, and the residue is far from the active center and substrate channel. Thus, mutation of the residue might have less effect on the stabilization of the protein. The aim this work is to improve the stability of HMFO for FDCA production, thus mutants of Q31 and N44G were selected, and their combinatorial mutations were utilized for further cataly performance investigation.

Catalytic Performance of Engineered HMFO
To determine the catalytic performance changes of the enhanced stable variants and the combinatorial mutation, the catalytic process of each variant at temperatures of 25 to 40 °C were investigated (Figure 4), and detailed FDCA yields after 24 h reaction were listed in Table  1. The results indicated that 8BxHMFO and all the mutants obtained similar yield of FDCA after 24 h reaction at 25 °C. However, when the temperature was increased to 30 and 35 °C,

Catalytic Performance of Engineered HMFO
To determine the catalytic performance changes of the enhanced stable variants and the combinatorial mutation, the catalytic process of each variant at temperatures of 25 to 40 • C were investigated (Figure 4), and detailed FDCA yields after 24 h reaction were listed in Table 1. The results indicated that 8BxHMFO and all the mutants obtained similar yield of FDCA after 24 h reaction at 25 • C. However, when the temperature was increased to 30 and 35 • C, two of the single site mutants showed better FDCA yields than others. When further increasing the temperature to 40 • C, all three kinds of mutants show better catalytic performance than 8BxHMFO. Compared to others, Q319K shows better performance under each of the reaction temperatures, and from the tendency of FDCA yields, it could be concluded that the improved HMFO mutants have better yields with shorter reaction time than the original 8BxHMFO within a temperature range of 30 to 40 • C. It is interesting that the catalytic result with combinatorial mutation N44G-Q319K were worse than the two of single site mutants; the reason for this may be that the combinatorial mutation influenced the structure of active site and thus decreased the catalytic performance. Similar results were also proved in other works of enzyme evolution [23].     The performance of each mutant with higher substrate concentration was also investigated. When HMF concentration reached 10 mM, 0.02 µM of HMFO could not provide the satisfactory yields of FDCA in the same reaction time (Figure 5a). Meanwhile, the pH stability of each candidate was also investigated (Figure 5b). The results indicate that the mutation has no effect on the pH stability improvement.
Catalysts 2021, 11, x FOR PEER REVIEW7 of 13 8BxHMFO. These parameters indicated that the FDCA yields improving with mutants Q319K and N44G (Figure 4) are due to the improvement of operational stability. Besides, fro the results of Figure 4, the combinatorial mutation of N44G-Q319K obtained less FDCA yiel than others. From the kinetic results of Table 2, it could be found that the catalyst performan of N44G-Q319K mutation was quite low, which could not provide satisfactory FDCA yield   [13]. In this work, with the increasing catalytic stability, only 0.02 µM HMFO mutan could reached FDCA yield higher than 95% with the reaction operated at 35 °C. All the resu indicate the operation stability of the mutants N44G and Q319K are improved significantly The kinetic parameters of each mutant were furthermore determined and compared with the original 8BxHMFO (Table 2). Except for the combinatorial mutation N44G-Q319K, the catalytic efficiencies of the two single site mutations are slightly lower than the original 8BxHMFO. These parameters indicated that the FDCA yields improving with mutants of Q319K and N44G (Figure 4) are due to the improvement of operational stability. Besides, from the results of Figure 4, the combinatorial mutation of N44G-Q319K obtained less FDCA yields than others. From the kinetic results of Table 2, it could be found that the catalyst performance of N44G-Q319K mutation was quite low, which could not provide satisfactory FDCA yields. The influence of incubation time on the enzyme activity at 35 • C was also investigated. The half-life of mutant N44G and Q319K are all longer than 72 h, with the activity of N44G reduced from 27.95 to 13.79 U/mg, and Q319K from 30.74 to 17.40 U/mg, respectively ( Figure 6). In comparison, the half-life of the original 8BxHMFO is only about 40 h. From the original publications of Methylovorus sp. MP688 8BxHMFO, the work was carried out with an enzyme concentration of 2 µM, and the yield of FDCA was around 98% at 25 • C, and greater than 95% at 40 • C [13]. In this work, with the increasing catalytic stability, only 0.02 µM HMFO mutants could reached FDCA yield higher than 95% with the reaction operated at 35 • C. All the results indicate the operation stability of the mutants N44G and Q319K are improved significantly.

Stability Improvement Mechanism Analysis
Compared to the original, the residue of Lys 319 forms more h residues around it in the mutant Q319K (Figure 7a), which microstructure [24], and similar results were also observed in the stru of N44G (Figure 7b). Root mean square fluctuations (RMSF) va flexibility and stability of individual protein residues. Normall corresponds to high flexibility and instability of the residue [25]. Tak flexibility index, the Q319K (Figure 8a,b) and N44G (Figure 8c,d) mu results compared to the 8BxHMFO. In addition, the RMSF values of 42 to 46 were decreased at 313K, indicating that the mutation conformation in these peptide fragments [17], and the simulation res the results of Figure 7.

Stability Improvement Mechanism Analysis
Compared to the original, the residue of Lys 319 forms more hydrogen bonds with the residues around it in the mutant Q319K (Figure 7a), which means a more stable microstructure [24], and similar results were also observed in the structure analysis of mutant of N44G (Figure 7b). Root mean square fluctuations (RMSF) value is a measure of the flexibility and stability of individual protein residues. Normally, a high RMSF value corresponds to high flexibility and instability of the residue [25]. Taking the RMSF value as a flexibility index, the Q319K (Figure 8a,b) and N44G (Figure 8c,d) mutants show lower RMSF results compared to the 8BxHMFO. In addition, the RMSF values of residues 314 to 325, and 42 to 46 were decreased at 313 K, indicating that the mutation adopted a more stable conformation in these peptide fragments [17], and the simulation results were consistent with the results of Figure 7.

In Silico Simulation
Firstly, a homology model of 8BxHMFO was obtained by the Automated Modeling Tool of Swiss Model Web Service (http://swissmodel.expasy.org/, accessed on 8 September 2021),

In Silico Simulation
Firstly, a homology model of 8BxHMFO was obtained by the Automated Modeling Tool of Swiss Model Web Service (http://swissmodel.expasy.org/, accessed on 8 September 2021), and the crystal structure of HMFO from Methylovorus sp. MP688 (PDB ID: 4UDQ) were applied as template [18]. Then, molecular dynamic simulation was carried out by YASARA software package (version 21.8.26, CMBI, Radboud University, Nijmegen, Netherlands, 8 September 2021) using the AMBER03 force field [26]. The simulation was carried out with water molecules added as solvent, at pH of 7.5 and temperature of 313 K for 10 ns. Then, the root-mean-square displacement (RMSD) of backbone atoms were calculated to evaluate the overall motion of 8BxHMFO [27]. Following this, the RMSF values of the residues were obtained through the simulation, and B-factors could be calculated according to Equation (1). Finally, the residues with top ten highest B-factor values in 8BxHMFO were selected as potential mutational sites for further investigation.

Saturation Mutagenesis Procedure
Site-saturation mutagenesis experiments were carried out using the pET-22b(+) vector harboring the gene encoding 8BxHMFO from Methylovorus sp. MP68 (NCBI reference sequence: WP_013440946.1) as template [13]. The primers containing the NNK degenerate codons (Table S1) were purchased from GenScript Corp (Nanjing, China). PCR reactions were performed in a total volume of 50 µL, including 0.2 µL template DNA, NNK degenerate forward and reverse primers (10 µM, 1 µL each), 1µL dNTP Mix (10 mM each), 1 µL Phanta Max Super-Fidelity DNA Polymerase. The PCR protocol consisted of initial denaturation for 30 s at 95 • C, followed by 35 cycles of 15 s at 95 • C, 64 • C for 10 s, and 72 • C for 1 min, and a final extension at 72 • C for 5 min.
A total of 50 µL PCR products were treated with DpnI restriction enzyme for 2 h at 37 • C for the removal of the parent plasmid. Then, the product was purified with a PCR purification kit (CWBIO, Co. Ltd., Beijing, China), and transferred to competent E. coli BL-21(DE3) by heat shock. The transformed bacteria were spread on a Luria-Bertani (LB) plate containing 50 µg/ mL ampicillin and incubated at 37 • C overnight, then obtained the corresponding mutation libraries.

Variation Screening of Mutation Libraries
For each transformant library, 120 colonies were picked and transferred to 150 µL LB-medium containing 50 µg/mL ampicillin in the wells of a 96-deep-well polypropylene plate. The plate was incubated at 37 • C, 200 rpm, overnight. After that, 10 µL of each preculture were used to inoculate into 500 µL of LB-medium distributed within a 96-deepwell plate, followed by cell growth, and induction was carried out at OD600 0.6-0.8, with 0.1 mM IPTG. Then, the colonies were incubated at 17 • C, 200 rpm, for 12 h. Finally, the cultures were pelleted by centrifugation (5000 rpm, 4 • C, 15 min).
For the variation screening, the cells were disrupted by the method of multigelation between −80 and 20 • C for three times, and followed by centrifugation with 5000 rpm at 4 • C for 5 min. Then, 100 µL of supernatants were transferred into a new black 96-wellplate. The plate was incubated at 45 • C for 1.5 h, then 5 µL of L-tyrosine (50 mM), 10 µL of HRP (1 mg/mL), 40 µL of HMF (50 mM), was added. The enzymatic catalysis was carried out with a total volume of 200 µL (pH 7.5, 50 mM phosphate buffer), at the reactions were carried out at 30 • C. The fluorescence spectra of the HRP-catalyzed reaction of tyrosine with H 2 O 2 were used to determine the reaction performance of HMFO (Scheme 1), the excitation and emission wavelengths of oxidation product were 330 nm and 415 nm, respectively. The fluorescence intensity was utilized to express the thermal stability changes of variations (Figure 2b), and the variations with high fluorescence intensity were selected for the further sequencing work.

HMFO Expression and Purification
For enzyme expression, an overnight culture of E. coli BL21(DE3) with the HMFO coding plasmid was diluted 1:100 in 1 L of Luria-Bertani broth containing 50 µg/mL ampicillin, and incubated at 37 • C until it reached an OD600 value of 0.8. Then, cells were induced with 0.1 mM IPTG and grown for 24 h at 18 • C. After that, E. coli cells were harvested by centrifugation at 8000× g for 10 min at 4 • C. The sediments of cultures were resuspended in potassium phosphate buffer (pH 8.0, 50 mM), and then disintegrated by ultrasound in an ice bath. After centrifugation (20,000× g for 15 min), the supernatant was purified following the method of His-tag purification described by Ueda et al. [28]. After the purification, the purity of the HMFO protein was >95%, and the activity of each HMFO was: 30.10 U/mg of 8BxHMFO, 30.74 U/mg of Q319K, 27.95 U/mg of N44G, and 26.85 U/mg of N44G-Q319K, respectively.
The HRP-Tyr spectrofluorimetric method was applied for the enzyme activity assay. The assay conditions were: 10 µL Tyr solution (50 mM), 10 µL HRP solution (1 mg/mL), 40 µL HMF (50 mM), HFMO solution 100 µL, the total volume was 200 µL in 50 mM PBS buffer. The definition of HFMO activity is the amount of enzyme that produces one micromole of H 2 O 2 per minute.

Enzyme Kinetics
The kinetic measurements were performed by monitoring the production of H 2 O 2 concentration changes with the method of HRP-Tyr assay by using substrate concentrations of 1.5 to 10 mM and purified HMFO variants (by the method of His-tag purification [28]) at an enzyme concentration of 0.02 µM. The K m , and k cat values were obtained from the non-linear regression fitting of the Michaelis-Menten curves (Figures S1-S4).

General Procedure for the HMFO Enzymatic Oxidation
Unless otherwise stated, the enzymatic reaction contains 5 mM of HMF (dissolved in 50 mM PBS buffer pH 7.5), 1% v/v dimethylsulfoxide (DMSO), 10 µL of catalase from Micrococcus lysodeikticus (100,000 U/mL) and HMFO (0.02 µM). Reactions were carried out at 25 • C to 40 • C, with shaking at 220 rpm for 24 h. The samples were taken at specific times and analyzed by high-performance liquid chromatography (HPLC) with the method described by Dijkman et al. [12].

Conclusions
The present work demonstrates the effects of direct evolution for improving the stability of HMFO. In many cases, flexible residues may cause protein unfolding leading to denaturation. The B-factor indicates flexible residues in a protein. The evolution strategy based on the B-factor analysis is simple and effective. After screening saturation mutation libraries of the sites with top ten higher B-factors, the mutants of Q319K and N44G show better catalytic stability in the temperature range of 25 to 40 • C. Especially with mutant Q319K as a catalyst, the FDCA yields were significantly increased to 98% from 85% (the original 8BxHMFO) at incubation temperature of 35 • C. Interestingly, residue Arg 252 shows the highest B-factor value. However, its mutagenesis shows less effect; this might be attributed to the fact that the location of Arg 252 was on a separated loop that is far from the active center and substrate channel, with less interaction with other secondary structures. Thus, the mutation of Arg 252 might have less effect on the stabilization of the protein. Meanwhile, the half-life of both mutants N44G and Q319K are around 72 h, which is about two folds longer than 8BxHMFO. This improvement may mainly be due to a higher degree of hydrogen bonds formed in the mutants, which increase the stability of the protein. Our results might be facilitating the design of more efficient thermostable biocatalysts, and are important for the further optimization and utilization of HMFO in the biotechnological processes.  Table S1: primers for saturation mutagenesis.
Author Contributions: Q.W. and S.J. performed the experiments. D.L. and J.L. contributed to molecular dynamic simulation. F.W. conceived, designed the experiments. L.L. and K.N. supervised all data and wrote the paper together with all authors. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by National Natural Science Foundation of China (21978017, 21978020). The authors would like to thank for the above funding.

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