Biochemical and Computational Insights on a Novel Acid-Resistant and Thermal-Stable Glucose 1-Dehydrogenase

Due to the dual cofactor specificity, glucose 1-dehydrogenase (GDH) has been considered as a promising alternative for coenzyme regeneration in biocatalysis. To mine for potential GDHs for practical applications, several genes encoding for GDH had been heterogeneously expressed in Escherichia coli BL21 (DE3) for primary screening. Of all the candidates, GDH from Bacillus sp. ZJ (BzGDH) was one of the most robust enzymes. BzGDH was then purified to homogeneity by immobilized metal affinity chromatography and characterized biochemically. It displayed maximum activity at 45 °C and pH 9.0, and was stable at temperatures below 50 °C. BzGDH also exhibited a broad pH stability, especially in the acidic region, which could maintain around 80% of its initial activity at the pH range of 4.0–8.5 after incubating for 1 hour. Molecular dynamics simulation was conducted for better understanding the stability feature of BzGDH against the structural context. The in-silico simulation shows that BzGDH is stable and can maintain its overall structure against heat during the simulation at 323 K, which is consistent with the biochemical studies. In brief, the robust stability of BzGDH made it an attractive participant for cofactor regeneration on practical applications, especially for the catalysis implemented in acidic pH and high temperature.


Heterologous Expression and Purification
The specific activity of the purified BzGDH was 194 ± 2 U·mg −1 at 25 °C using NAD (nicotinamide adenine dinucleotide) as a cofactor. SDS-PAGE (sodium dodecyl sulfate polyacrylamide gel electrophoresis) analysis showed a homogeneous band corresponding to 30 kDa (Figure 3). By using gel filtration chromatography through a Zorbax Bio-series GF-450 column, the molecular weight of Figure 2. Unrooted phylogenetic tree of GDHs. The phylogenetic tree was constructed using the neighbor joining method [30] in MEGA7 software [31], with a bootstrap test of 1000 replicates. The evolutionary distances were computed using the Poisson correction method [32] and are in the units of the number of amino acid substitutions per site.

Figure 2.
Unrooted phylogenetic tree of GDHs. The phylogenetic tree was constructed using the neighbor joining method [30] in MEGA7 software [31], with a bootstrap test of 1000 replicates. The evolutionary distances were computed using the Poisson correction method [32] and are in the units of the number of amino acid substitutions per site.

Effects of pH and Temperature on the Activity and Stability
BzGDH exhibited activity at a wide pH range from 4.0 to 10.5, and displayed maximum activity at pH 9.0 in Tris-HCl buffer among all buffers. Actually, the chemical composition of the sodium citrate, sodium phosphate, and Tris-HCl buffers, showed no significant influence on the specific activity of the enzyme, and the differences in the specific activity are mainly caused by the change of the pH of the solution. However, a significant decrease of specific activity was observed in the glycine-NaOH buffer at pH values of 8.5 and 9.0 when compared to those of the same pH values of the Tris-HCl buffer, indicating that glycine might inhibit the activity of BzGDH. Surprisingly, the optimum pH was determined as 9.5 in Glycine-NaOH buffer (Figure 4a), which is inconsistent with the maximum activity pH of 9.0. A reasonable explanation for this discrepancy is that the observed activity of the enzyme is not only affected by the pKa of its catalytic residues which played critical roles on the activity, but is influenced by the stability of the enzyme which might be unstable at its optimum pH (Figure 4b), and is even sometimes affected by the chemicals in the buffer such as glycine in this case. In regards to its pH stability, BzGDH was stable over a broad pH range, especially in the acidic region, which could maintain around 80% of its initial activity in the pH range of 4.0-8.5 after incubating for 1 hour (Figure 4b).
the maximum activity pH of 9.0. A reasonable explanation for this discrepancy is that the observed activity of the enzyme is not only affected by the pKa of its catalytic residues which played critical roles on the activity, but is influenced by the stability of the enzyme which might be unstable at its optimum pH (Figure 4b), and is even sometimes affected by the chemicals in the buffer such as glycine in this case. In regards to its pH stability, BzGDH was stable over a broad pH range, especially in the acidic region, which could maintain around 80% of its initial activity in the pH range of 4.0-8.5 after incubating for 1 hour (Figure 4b). As demonstrated in Figure 4c, the optimum catalytic temperature of BzGDH was determined as 45 °C. The activity of BzGDH decreased linearly from 45 to 65 °C and could not be measurable at 75 °C. In consistent with its higher optimum reaction temperature, the recombinant enzyme also possessed good thermal stability, which was stable after incubation at temperatures below 50 °C for 30 min and still maintained 50% of its initial activity after incubation at 65 °C for 30 min (Figure 4d). BzGDH As demonstrated in Figure 4c, the optimum catalytic temperature of BzGDH was determined as 45 • C. The activity of BzGDH decreased linearly from 45 to 65 • C and could not be measurable at 75 • C. In consistent with its higher optimum reaction temperature, the recombinant enzyme also possessed good thermal stability, which was stable after incubation at temperatures below 50 • C for 30 min and still maintained 50% of its initial activity after incubation at 65 • C for 30 min (Figure 4d). BzGDH exhibited superior thermal stability to its homologous counterparts, BgGDH [23] and BcGDH [24], which were almost completely inactivated after incubation at 50 • C without any protective agent.
Since stability is an indispensable characteristic for the utilization of enzymes in real life, the considerable stability of BzGDH against both heat and acid made it a very promising candidate in practical application in harsh conditions.

Substrate Specificity and Steady-State Kinetics
As shown in Table 1, the substrate spectrum of BzGDH was similar to that of BcGDH. However, both BzGDH and BcGDH displayed stricter substrate specificity toward various sugars than that of BgGDH, especially for galactose and mannose, indicating that BzGDH could be a potential diagnostic reagent for blood glucose measurement as well as BcGDH.
The steady-state kinetic constants of BzGDH were determined by using a nonlinear fitting plot ( Table 2). Although BzGDH had similar k cat values for both NAD and NADP, the K m value for NADP was 5.6-fold higher than that for NAD, indicating that BzGDH preferred NAD rather than NADP as the cofactor. The cofactor preference of BzGDH resembled that of BmGDHIII, BmGDHIV [4], and BtGDH [3], while BmGDH, BmGDHI, BmGDHII [5], and BgGDH [23] preferred NADP.

Homology Modeling and Electrostatic Potential Analysis
The quaternary structure of BzGDH was constructed by SWISS-MODEL [33] and evaluated by ProSA-web [34] and PROCHECK [35]. Both of the Z-score and Ramachandran plot statistics indicated that the dimensional structure of BzGDH ( Figure 5a) had been modeled reasonably (Table 3). To investigate the electrostatic potential of BzGDH, the model of BzGDH was subjected to the software APBS [36] and PyMOL (The PyMOL Molecular Graphics System, Version 1.7 Schrödinger, LLC. available online: http://pymol.org/), to generate the electrostatic potential molecular surface. As shown in Figure 5, the contact surfaces of subunits AB, AC, and AD circled by black ellipses are mainly constituted by non-polar amino acid residues and are surrounded by acidic amino acid residues. The non-polar areas can maintain their electrically neutral state in either acidic or alkaline solutions, whereas the acidic areas would be negatively charged in alkaline solutions, leading to the mutual repulsion between subunits. Therefore, the acid-resistance of BzGDH could be explained by the electrostatic potential of contact surfaces between subunits, as well as BcGDH [24].  (c) Electrostatic potential surface representation of the interface between subunits A and C.
(d) Electrostatic potential surface representation of the interface between subunits A and D. Subunits ABCD were labeled using the corresponding capital letters nearby, respectively. Positive, negative, and neutral electrostatic potential surfaces are rendered by blue, red, and white, respectively. The non-polar regions of the contact surfaces of subunits AB, AC, and AD were circled by dashed ellipses. (c) Electrostatic potential surface representation of the interface between subunits A and C; (d) Electrostatic potential surface representation of the interface between subunits A and D. Subunits ABCD were labeled using the corresponding capital letters nearby, respectively. Positive, negative, and neutral electrostatic potential surfaces are rendered by blue, red, and white, respectively. The non-polar regions of the contact surfaces of subunits AB, AC, and AD were circled by dashed ellipses.

Global Structure Stability
To study the stability and mobility of BzGDH, the model was subjected to a 20-ns MD simulation at 323 K. The stability of BzGDH was analyzed by the all-atom and backbone-atom root mean square deviation (RMSD), respectively, both of which increased from the beginning of the simulation and reached an equilibrium state at about 10 ns (Figure 6a), suggesting no significant structural changes for BzGDH during the simulation. In addition, the radius of gyration, the hydrogen bonds of intra-protein, and the solvent accessible surface area (SASA) of BzGDH all displayed steadily dynamic changes against time (Figure 6b-d), further confirming the stable global behavior of BzGDH during the simulation at 323 K. mean square deviation (RMSD), respectively, both of which increased from the beginning of the simulation and reached an equilibrium state at about 10 ns (Figure 6a), suggesting no significant structural changes for BzGDH during the simulation. In addition, the radius of gyration, the hydrogen bonds of intra-protein, and the solvent accessible surface area (SASA) of BzGDH all displayed steadily dynamic changes against time (Figure 6b-d), further confirming the stable global behavior of BzGDH during the simulation at 323 K.

Structure Flexibility
The conformational flexibility of BzGDH was assessed using the root mean square fluctuation (RMSF) of C-alpha (Cα) atoms per residue. Generally, regular secondary structure regions display tiny fluctuations with small RMSF values during the simulation, whereas prominent fluctuations with large RMSF are observed for irregular secondary structure regions such as terminal or loop regions, which often bear certain function of proteins. As shown in Figure 6e, regions involved in coenzyme binding (39-55) and substrate binding (190-210) of each subunit are more flexible with large RMSF values than other regions. The RMSF values were converted to B-factors using the equation: to visualize global structure rigidity and flexibility of BzGDH. As shown in Figure 7, most regions of BzGDH are rigid, except for the aforementioned flexible regions, indicating that the enzyme is stable during the simulation at 323 K.

Structure Flexibility
The conformational flexibility of BzGDH was assessed using the root mean square fluctuation (RMSF) of C-alpha (Cα) atoms per residue. Generally, regular secondary structure regions display tiny fluctuations with small RMSF values during the simulation, whereas prominent fluctuations with large RMSF are observed for irregular secondary structure regions such as terminal or loop regions, which often bear certain function of proteins. As shown in Figure 6e, regions involved in coenzyme binding (39-55) and substrate binding (190-210) of each subunit are more flexible with large RMSF values than other regions. The RMSF values were converted to B-factors using the equation: to visualize global structure rigidity and flexibility of BzGDH. As shown in Figure 7, most regions of BzGDH are rigid, except for the aforementioned flexible regions, indicating that the enzyme is stable during the simulation at 323 K. In addition to the observation of RMSF, the bond-specific fluctuations in protein structure can further be captured by the Lipari-Szabo order parameter S 2 [37], which provide an intuitive description of the amplitude of spatial restriction of the internal motions of the bond vectors on a fast timescale from picosecond to nanosecond (ps-ns). More specifically, S 2 represents the component of the H-X bond vector autocorrelation function which is dissipated by global molecular tumbling, while (1 − S 2 ) characterizes the bond vector orientational disorder arising from internal motion occurring more rapidly than the molecular tumbling. The S 2 order parameter can range from 0 to 1, with 1 corresponding to a rigid bond vector (completely restricted) and 0 corresponding to the highest degree of disorder for a bond vector (completely isotropic). Higher order parameters (0.85) were observed in the regions of secondary structure, while unstructured regions showed lower order parameters (0.4-0.6).
The order parameter S 2 of the main chain N-H bonds of BzGDH has been calculated based on the equilibrium MD trajectories. The average value of the order parameter S 2 , over all residues, is 0.86 for BzGDH. The most flexible region that showed lower S 2 of each subunit is the substrate binding domain, with residues Lys 179, Gly 180, Arg 182, Asn 184, Asn 185, Ala 190, Asn 196, and Asp 202 involved, indicating that these residues exhibit considerable disorder on the ps-ns timescale. Similarly, residues Gln 257, Ala 258, and Gly 259 in the C-terminal region of the protein have low order parameters, also implying that this region is disordered on the ps-ns timescale. Indeed, the order parameter revealed that these regions are flexible on the ps-ns timescale, with the fluctuations functioning to allow substrate access to and release of products from the active site. The results of the computation of the order parameters are in considerable agreement with the RMSF profiles, with the greatest flexibility occurring in loop regions, while other secondary structural elements are more constrained.
description of the amplitude of spatial restriction of the internal motions of the bond vectors on a fast timescale from picosecond to nanosecond (ps-ns). More specifically, S 2 represents the component of the H-X bond vector autocorrelation function which is dissipated by global molecular tumbling, while (1−S 2 ) characterizes the bond vector orientational disorder arising from internal motion occurring more rapidly than the molecular tumbling. The S 2 order parameter can range from 0 to 1, with 1 corresponding to a rigid bond vector (completely restricted) and 0 corresponding to the highest degree of disorder for a bond vector (completely isotropic). Higher order parameters (0.85) were observed in the regions of secondary structure, while unstructured regions showed lower order parameters (0.4-0.6).
The order parameter S 2 of the main chain N-H bonds of BzGDH has been calculated based on the equilibrium MD trajectories. The average value of the order parameter S 2 , over all residues, is 0.86 for BzGDH. The most flexible region that showed lower S 2 of each subunit is the substrate binding domain, with residues Lys 179, Gly 180, Arg 182, Asn 184, Asn 185, Ala 190, Asn 196, and Asp 202 involved, indicating that these residues exhibit considerable disorder on the ps-ns timescale. Similarly, residues Gln 257, Ala 258, and Gly 259 in the C-terminal region of the protein have low order parameters, also implying that this region is disordered on the ps-ns timescale. Indeed, the order parameter revealed that these regions are flexible on the ps-ns timescale, with the fluctuations functioning to allow substrate access to and release of products from the active site. The results of the computation of the order parameters are in considerable agreement with the RMSF profiles, with the greatest flexibility occurring in loop regions, while other secondary structural elements are more constrained.

Essential Dynamics
To reveal the concerted fluctuations of BzGDH over time and spatial scales, essential dynamics (ED) is employed to extract information from sampled conformations over the molecular dynamics trajectory [38]. Practically, the essential dynamics of a protein is obtained by performing principal component analysis (PCA), which is a multivariate statistical technique involving diagonalization of the covariance matrix ( Figure 8) constructed from atomic displacements of Cα atoms, to reduce the number of dimensions required to describe protein dynamics and yield a set of eigenvectors that provide information about collective motions of the protein [39].
To reveal the concerted fluctuations of BzGDH over time and spatial scales, essential dynamics (ED) is employed to extract information from sampled conformations over the molecular dynamics trajectory [38]. Practically, the essential dynamics of a protein is obtained by performing principal component analysis (PCA), which is a multivariate statistical technique involving diagonalization of the covariance matrix ( Figure 8) constructed from atomic displacements of Cα atoms, to reduce the number of dimensions required to describe protein dynamics and yield a set of eigenvectors that provide information about collective motions of the protein [39]. The eigenvectors represent the directions of motion, and the corresponding eigenvalues represent the amount of motion along each eigenvector, where larger eigenvalues describe motions on larger spatial scales. Generally, the first 10-20 eigenvectors are enough to capture the principal motions of the protein and describe more than 90% of all cumulative protein fluctuations [40]. However, it can be seen that only 14.3% of the total Cα motion can be explained by the first two eigenvectors, even the first 20 eigenvectors merely contribute for 51.2% of the total Cα motion from Figure 9a. This shows that most of the internal motions of BzGDH are not confined within a subspace of small dimension, and no obvious collective motion of the backbone of BzGDH is observed from the MD simulation performed at 323 K, reflecting that the enzyme can maintain its overall structure against heat, which is in accordance with the biochemical experiments. Figure 9b shows the trajectory projected on the plane defined by the first two principal eigenvectors. The trajectories filled most of the expected ranges, suggesting the deficiency of a coupled force field, which leads to independent motions. The trajectories were projected onto the individual eigenvectors against time to further investigate the motion along the eigenvector directions. It is clear from Figure 10 that the fluctuations of the first six eigenvectors are relatively large, whereas those of the subsequent eigenvectors become successively flat, indicating that the The eigenvectors represent the directions of motion, and the corresponding eigenvalues represent the amount of motion along each eigenvector, where larger eigenvalues describe motions on larger spatial scales. Generally, the first 10-20 eigenvectors are enough to capture the principal motions of the protein and describe more than 90% of all cumulative protein fluctuations [40]. However, it can be seen that only 14.3% of the total Cα motion can be explained by the first two eigenvectors, even the first 20 eigenvectors merely contribute for 51.2% of the total Cα motion from Figure 9a. This shows that most of the internal motions of BzGDH are not confined within a subspace of small dimension, and no obvious collective motion of the backbone of BzGDH is observed from the MD simulation performed at 323 K, reflecting that the enzyme can maintain its overall structure against heat, which is in accordance with the biochemical experiments. Figure 9b shows the trajectory projected on the plane defined by the first two principal eigenvectors. The trajectories filled most of the expected ranges, suggesting the deficiency of a coupled force field, which leads to independent motions. The trajectories were projected onto the individual eigenvectors against time to further investigate the motion along the eigenvector directions. It is clear from Figure 10 that the fluctuations of the first six eigenvectors are relatively large, whereas those of the subsequent eigenvectors become successively flat, indicating that the motions belong to the last four eigenvectors have reached their equilibrium fluctuation, which cannot be used to describe the motions of the system. Due to the limitation of hardware, such simulations may not capture the essential motions related to function at much longer timescales. Improvements in computational power will fill the gap between reality and simulation. motions belong to the last four eigenvectors have reached their equilibrium fluctuation, which cannot be used to describe the motions of the system. Due to the limitation of hardware, such simulations may not capture the essential motions related to function at much longer timescales. Improvements in computational power will fill the gap between reality and simulation.   motions belong to the last four eigenvectors have reached their equilibrium fluctuation, which cannot be used to describe the motions of the system. Due to the limitation of hardware, such simulations may not capture the essential motions related to function at much longer timescales. Improvements in computational power will fill the gap between reality and simulation.

Strains, Plasmids, and Chemicals
Strain Bacillus sp. ZJ isolated from the soil near Yuhangtang River in Hangzhou, China, was used as a source for retrieving glucose 1-dehydrogenase. Escherichia coli (E. coli) DH5α and BL21 (DE3), expression vector pET28a (+), and Ni-NTA resin were purchased from Invitrogen. Taq DNA polymerase, PrimeSTAR HS DNA Polymerase, T4 DNA ligase, NdeI, and BamHI were purchased from TaKaRa. Genomic DNA, plasmid and gel extraction kits were purchased from Axygen. All other chemicals were of analytical grade.

Cloning of the Bzgdh Gene and Sequence Analysis
The gene bzgdh was amplified by using genomic DNA of Bacillus sp. ZJ as a template with the forward primer 5 -GGAATTCCATATGTATAGTGATTTAGCAGG-3 and the reverse primer 5 -CGGGATCCTATTACCCACGCCCAGC-3 , which carried cutting sites of NdeI and BamHI (underlined), respectively. The amplified fragments were digested with NdeI and BamHI simultaneously, then purified by using a gel extraction kit prior to ligate with the pre-digested vector pET-28a (+). The recombinant plasmid harboring gene bzgdh was transformed into competent cells of E. coli DH5α for sequencing.
Homologous searches in GenBank were performed using the BLAST server (available online: http://blast.ncbi.nlm.nih.gov). Alignment of multiple protein sequences was conducted by using the Clustal X 2.0 program [28] and rendered by ESPript [29]. The phylogenetic tree was constructed using the neighbor-joining method in MEGA7 [31], with a bootstrap test of 1000 replicates.
The nucleotide sequence for GDH of Bacillus sp. ZJ was deposited in GenBank under accession number KJ701281.

Expression and Purification of Recombinant BzGDH
The recombinant plasmid was transformed into competent cells of E. coli BL21 (DE3) for expression. The recombinant cells were cultivated in Luria-Bertani broth containing 50 µg kanamycin·mL −1 at 37 • C with a shaking speed of 250 rpm. The expression of recombinant protein was induced by adding 0.5 mM of IPTG to the medium when the OD 600 of the culture reached 0.5-0.8, followed by another 12 h incubation at 25 • C with a shaking speed of 200 rpm. The cells were harvested by centrifugation at 10,000× g for 10 min at 4 • C and were washed with the binding buffer (50 mM NaH 2 PO 4 , 500 mM NaCl, 20 mM imidazole, pH 8.0), and then lysed by ultrasonication. The cell debris was removed by centrifugation at 15,000× g for 30 min at 4 • C, and then the supernatant was loaded onto a column containing pre-equilibrated Ni-NTA resin. The column was washed with binding buffer and subsequently eluted with elution buffer (50 mM NaH 2 PO 4 , 500 mM NaCl, 250 mM imidazole, pH 8.0). The eluted enzyme was desalted and concentrated by ultrafiltration and stored at −80 • C in 25 mM sodium phosphate buffer (pH 6.5) with 30% of glycerol contained. The protein concentration was determined by Bradford's method using bovine serum albumin as the reference.
Denaturing discontinuous polyacrylamide gel electrophoresis was performed on a 5% stacking gel and a 12% separating gel. The native molecular weight of GDH was determined by size-exclusion chromatography according to the protocol of the manufacture (Zorbax Bio-series GF-450, Agilent, Santa Clara, CA, USA), using lysozyme (14.3 kDa), chicken ovalbumin (45 kDa), bovine serum albumin fraction V (67 kDa), and goat IgG (150 kDa) as standards.

Enzyme Activity Assay
Glucose dehydrogenase activity was determined by assaying the absorbance of NADH at 340 nm in 100 mM sodium phosphate (pH 8.0) containing 200 mM glucose and 1 mM NAD at 25 • C. All measurements were conducted in triplicate. One unit of enzyme activity was defined as the amount of the enzyme that catalyzed the formation of 1 µmol of NADH per minute.

Effects of pH and Temperature on the Activity and Stability of BzGDH
The optimum pH of BzGDH was measured at pH ranging from 4.0 to 10.5 at 25 • C. The effect of pH on the stability of BzGDH was determined by measuring the residual activity after incubating BzGDH in buffers with different pH values for one hour at 25 • C.
The optimal temperature of BzGDH was determined at different temperatures (25-75 • C) in phosphate buffer at pH 7.0. The thermal stability of BzGDH was assayed by measuring the residual activity after incubating BzGDH at different temperatures (25-75 • C) in phosphate buffer at pH 7.0 for 30 min.

Substrate Specificity of BzGDH
The substrate specificity of BzGDH was determined by the aforementioned enzyme activity assay, except that glucose was replaced by sucrose, lactose, maltose, xylose, galactose, mannose, fructose, and arabinose, respectively.

Steady-State Kinetics of BzGDH
In order to obtain the kinetic constants for the coenzyme, 200 µM of glucose was employed as the substrate and 0.01 to 0.2 mM NAD and NADP were used as the coenzymes, respectively. For analysis of the kinetics for glucose, 1 mM NAD was used as a cofactor, 1 to 200 mM glucose was used as the substrate. GDH activity was measured as described above. The kinetic constants were determined by using a nonlinear fitting of the Michaelis-Menten equation: where [S] is the concentration of the cofactor or substrate, K m is the Michaelis constants for the cofactor or substrate, v is the reaction velocity, and V max is the maximum reaction velocity. The turnover number k cat was calculated by the equation: where [E] is the concentration of the enzyme.

Homology Modeling and Electrostatic Potential of BzGDH
The crystal structure of glucose 1-dehydrogenase from Bacillus megaterium IWG3 (PDB code: 1GEE, 1.60 Å) [41], which shares 88.12% identity with BzGDH, was served as the template for homology modeling of BzGDH. The three-dimensional model of BzGDH was constructed by using the SWISS-MODEL [33]. Precise evaluation of the model quality was performed using ProSA-web [34] and PROCHECK [35]. The structure for electrostatics calculations was prepared by PDB2PQR [42] using the AMBER force field and assigned protonation states at pH 7.0. The electrostatic potential of BzGDH was calculated by APBS [36] using the linearized Poisson-Boltzmann equation (lpbe) at 298 K with the monovalent ion concentration of 0.1 M. The dielectric constants of protein and solvent were set as 2.0 and 78.0, respectively. The electrostatic potential molecular surface was represented by PyMOL.

Molecular Dynamic Simulations of BzGDH
The constructed model of BzGDH was subjected to the software package GROMACS 5.0.2 [43], with the AMBER99SB [44] force field adopted, for molecular dynamics simulations. The model was first placed into the center of a virtual cubic box with a side length of 11.049 nm and solvated with 39,486 TIP3P water molecules. The pH condition was 7.0 according to the ionization state of the protein with a charge of −20, and twenty Na + ions were added to the water box as counter ions to neutralize the negative charge of the entire system. Bond lengths were constrained by the LINCS algorithm to ensure covalent bonds to maintain their correct lengths during the simulation. Energy minimization of the system was conducted using the steepest descent algorithm for 5000 steps, followed by a 500-ps equilibration simulation with harmonic position restraints on the heavy protein atoms to equilibrate the solvent molecules around the protein. Subsequently, a 2-ns simulation without position restraints was conducted to equilibrate the entire system. Finally, the production simulation was performed for 20 ns at the target temperature. All simulations were performed under the NPT ensemble with periodic boundary conditions and a time step of 2 fs. The temperature of the system was kept at 323 K using the v-rescale method, and the pressure was kept at 1 bar using the Parrinello-Rahman method.
According to the RMSD profile of BzGDH, trajectories that reached the equilibrium state (10-20 ns) were used for further analysis. Principal component analysis was conducted to identify the direction and amplitude of the most prominent characteristics of the motions of BzGDH along the simulation trajectory. Generalized order parameters S 2 , employed as a measure of the degree of spatial restriction of motion, were also calculated for the N-H bonds of BzGDH.

Conclusions
In this study, a novel NAD(P)-dependent glucose 1-dehydrogenase from Bacillus sp. ZJ has been extensively characterized, with remarkable acidic tolerance and thermal stability. To better understand the stability feature of BzGDH against the structural context, molecular dynamics simulation was conducted to investigate the conformational flexibility and fluctuations of BzGDH over time and spatial scales. Analysis of the trajectory shows that BzGDH is stable and can maintain its overall structure against heat during the simulation at 323 K, which is in accordance with the biochemical studies. In brief, the robust stability of BzGDH made it a promising participant for cofactor regeneration in practical applications, especially for catalysis implemented in acidic pH and high temperature.