The gene for HvG6PDH (Uniprot code D4GS48), with the optimized sequence for expression in Escherichia coli, was synthesized by GenScript USA Inc., (Piscataway, NJ, United States) and subsequently cloned into the vector pET-TEV-28a, which includes a His-tag in the N-terminal and resistance to kanamycin (KAN). The cells of the E. coli BL21(DE3) strain were transformed with the expression vector and 50 ml of Luria–Bertani (LB), pH 7.0, supplemented with 30 μg ml^–1 of KAN and grown overnight with shaking at 200 rpm. Then, 500 ml of the LB/KAN medium was inoculated and supplemented with 1 mM IPTG when the culture reached 0.8–1.2 units of OD₆₀₀ and grown overnight at 37°C. The cells were centrifuged for 10 min at 6,130 g, and the pellet was resuspended in 30 ml of lysis buffer (2 M KCl, 50 mM Tris–HCl, pH 7.8, 1 mM EDTA, 5 mM 2-mercaptoethanol, and 0.1% v/v Triton x-114). The resuspended pellet was sonicated with 15 pulses of 20 s at 40% of amplitude, and pauses of 40 s between each pulse were performed. The lysate was centrifuged at 38,360 g for 10 min to recover the insoluble fraction corresponding to the inclusion bodies. Afterward, the inclusion bodies were resuspended in 12.5 ml of solubilization buffer (8 M urea, 50 mM Tris–HCl, pH 7.8, 1 mM EDTA, and 50 mM dithiothreitol) and incubated for 1 h at 37°C. After solubilization, the supernatant (solubilized inclusion bodies) was centrifuged at 38,360 g for 30 min to recover the soluble fraction. HvG6PDH-solubilized inclusion bodies were filtered through a 0.2-μm syringe filter, refolded by rapid dilution in 237.5 ml of renaturation buffer (2 M KCl, 50 mM Tris–HCl, pH 7.8, and 5 mM 2-mercaptoethanol), and incubated at 4°C overnight with agitation. The refolded inclusion bodies were centrifuged at 38,360 g for 20 min, and the supernatant was filtered through a 0.2-μm syringe filter and loaded on a 5-ml HisTrap-HP column (GE Healthcare) equilibrated with renaturation buffer (2 M KCl, 50 mM Tris–HCl, pH 7.8, and 5 mM 2-mercaptoethanol). The column was washed with 10 volumes of buffer A, and protein elution was performed with a linear gradient between 50 and 250 mM imidazole in 50 ml. Fractions with G6PDH activity were pooled, and protein purity was evaluated by sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE). The samples were precipitated with 20% trichloroacetic acid to remove the excess of KCl.

Enzyme activity was determined at 25°C by following NADH formation at 340 nm in a buffer of 100 mM Tris–HCl, pH 8.5, 2 M KCl, 3 mM NAD⁺, and 50 mM G6P. A UV/visible Synergy 2 spectrophotometer (BioTek, Winooski, VT, United States) and 96-well plates (model 269620; Nunc, Rochester, NY, United States) were used for measurements. An extinction coefficient of 6.22 mM^–1 × cm^–1 for NADH, and a path length correction at 1 cm was used for specific activity (U × mg^–1) determination. Enzyme unit (U) was defined as micromole of product per minute. To ensure initial velocity conditions, 4–6 nM of protein concentration was used. We determined the kinetic parameters for G6PDH by varying the substrate concentration at a fixed concentration of the co-substrate at pH 7.5 and NaCl and KCl concentrations between 0.5–3.0 and 0.5–2.5 M, respectively. The equation of Michaelis–Menten was adjusted to initial velocities by the non-linear method using the GraphPad Prism 8 software.

The ADP-dependent glucokinase activity was determined using the ADP-dependent GK form Thermococcus litoralis (TlGK) by a coupled assay at 25°C following NADH formation at 340 nm in a buffer of 100 mM Hepes, pH 7.5, 3 mM NAD⁺, 5 mM glucose, 1 mM ADP, and 6 mM MgCl₂. The enzyme concentration was 12 nM, and the KCl concentrations varied between 0.2 and 2.0 M.

A homology model of HvG6PDH was constructed using MODELLER, v9.24 (Webb and Sali, 2017). The reported structure of the SDR uronate dehydrogenase from Agrobacterium fabrum (PDB: 3RFV) was used as a template since it is the phylogenetically closer protein to HvG6PDH (29.2% identity) with determined three-dimensional structure in the presence of ligands. The 10 models obtained were evaluated based on DOPE score. The model with the best score was used to remodel the non-conserved loop between residues 69 and 80 by ab initio modeling using the DaReUS-Loop web tool (Karami et al., 2019), generating 10 loop variants. These models were evaluated using the ProSA-web tools (Wiederstein and Sippl, 2007), Procheck (Laskowski et al., 1993), and Verify3D (Eisenberg et al., 1997) to choose the best one. NAD⁺ and the sugar were added to the best resulting model using MODELLER, v9.24 (Webb and Sali, 2016) software and based on the ligands present in PDB 3RFV. The G6P was built based on the galactaro-1,5-lactone present in the 3RFV PDB structure, changing manually the atomic differences in PyMOL. Finally, the loop β6-α10 (residues from 178 to 190) was re-modeled by the ab initio function of MODELLER to solve clashes with G6P.

From the homology model of HvG6PDH, we prepared two systems of molecular simulations, one with 1.5 M of KCl and another without salt. The protonation states were assigned using Propka Software (Olsson et al., 2011) implemented in Maestro, version 9.7 (Schrodinger, LLC, New York, NY, United States). The topology and coordinate input files were prepared using the tLEaP software, and protein parameters were obtained using the ff14SB force field. The parameters for G6P were estimated with Antechamber/Chimera using the GAFF force field, whereas for NAD⁺ the parameters were taken from the Bryce group database (Ryde, 1995). The protein/complex was immersed in a periodic TIP3P water box composed of 25,987 water molecules with 12-Å extension and 58 K⁺ ions to neutralize the charge of the model. For the 1.5 M salt condition, 729 K⁺ ions were added along with the same amount of Cl^– counterions.

Once the systems were prepared, a restricted minimization protocol was applied [adapted from Shaw et al. (2010)] with a gradient of spatial restrictions, both to the protein and ligands, to minimize separately the solvent and the protein. These restrictions were gradually released (300, 200, and 10 kcal/mol Ų) in three minimization steps, ending with generally 40,000 minimization cycles of the entire system without restrictions. All minimizations were performed in AMBER using pmemd.MPI on a 64-core server. The initial steps utilized the steepest descent method and then switched to a conjugate gradient.

The minimized systems were equilibrated in two steps under the canonical ensemble (NVT: moles, volume, and temperature conserved) conditions; first, for solvent arrangement, the systems were heated to 298 K during 200 ps using a Langevian protocol and 500 kcal/mol Ų restrictions over the protein and ligand atoms. Afterward, the last coordinates of this system were used as a starting point for heating the whole system to 298 K, with no restriction to protein or ligands. This heating was maintained for 1,000 ps, being 75% of this time a ramp from 0 to 298 K.

Then, both systems (low and high salt) were equilibrated under isothermal–isobaric (NPT) ensemble (moles, pressure, and temperature conserved) conditions at a constant temperature of 298 K for 35 ns. From this step, productive trajectories were performed over 100 ns for both conditions, using pmemd.cuda – available in the Amber18 suite, with a time step of 2 fs and periodic border conditions. For each system, three replicas were performed. The energy interactions between the protein and ligands as well as root mean square deviation (RMSD), root mean square fluctuation (RMSF), and radius of gyration (R_g) were calculated using CPPTRAJ V4.14.0 (AmberTools V19.12) (Roe and Cheatham, 2013).

To preserve a competent catalytic conformation, a second batch of molecular dynamics, using a protocol like the one described above, was applied. In this case, an extra minimization step was included before the final minimization, and a new set of restrictions (10 kcal/mol Ų) to only G6P, NAD⁺, and to the catalytic residue (tyrosine 152) was applied.

A restriction of 10 kcal/mol Ų to ligands and to the catalytic residue Y152 was applied during the second step of the NVT equilibration in both systems. Finally, during the 35 ns NPT equilibration and the 100 ns productive trajectories, a spatial restriction of 1 kcal/mol Ų to G6P and to the hydroxyl group of the Y152 catalytic residue was included.

The amino acid distribution in the protein structure was evaluated through the GetArea (Fraczkiewicz and Braun, 1998), and residues were clustered based on solvent accessible surface (SAS). If the SAS is ≥20 Ų, the residue is considered exposed to the surface (external shell); if not, it is considered as part of the protein core (internal shell). The electrostatic surface potential was calculated using the APBS tools implemented in PyMOL 2.0 (Baker et al., 2001).

We searched for homologous sequences of SDR-G6PDH from H. volcanii by a PSI-blast search, with five iterations, in the non-redundant NCBI database and in the UniProtKB database. A non-redundant set with 80% identity was chosen using Multiseq (Roberts et al., 2006), and sequences of characterized SDR members as glucose dehydrogenase (GlcDH) from E. coli, human GDP-mannose dehydratase, UDP-glucose-4,6-epimerase from Arabidopsis thaliana and E. coli, L-arabinose from H. volcanii, and uronate dehydrogenase from A. fabrum were manually included to identified groups of homolog sequences with an unknown function. A multiple sequence alignment, using ClustalX, was constructed with the initial set of sequences. All these sequences were aligned based on the structural profile of the following structures (PDB ID 3RFV, 3AY3, and 6D9Y). The misaligned positions were manually corrected, and a maximum-likelihood tree was constructed using the PhyML server (Guindon et al., 2010), using the VT evolutive model (Müller and Vingron, 2001) according to AIC-SMS (Lefort et al., 2017).

As was expected, the expression of HvG6PDH in E. coli results in inclusion bodies formation, probably due to the requirement of high ionic strength conditions to be properly folded. Therefore, HvG6PDH was purified from inclusion bodies using a refolding protocol as described in the section “Materials and Methods.” The procedure was accomplished with one chromatographic step in a Ni-NTA column, where 5 mg of highly pure HvG6PDH with a specific activity of 134 U/mg was obtained from 500 ml of culture (Figure 1).

We assessed the effect of salt concentration on the optimum pH and kinetic parameters. A significant increase in enzyme activity was observed over 0.5 M of salt, although the optimum pH (between 8 and 10) was not affected by salt (Supplementary Figure 1). The saturation curves for both substrates, G6P and NAD⁺, were determined in the presence of different KCl and NaCl concentrations. In all the conditions tested, the enzyme presents a hyperbolic behavior for both substrates (Supplementary Figure 2). The determination of the kinetics parameters show that an increase in KCl concentrations produced an increment in the V_max value, while the K_M for either NAD⁺ or G6P was not significantly affected (Table 1). A similar effect on the kinetic parameters was observed in the presence of NaCl (Supplementary Figures 2C,D) except for 0.5 M NaCl, where saturation with G6P was not achieved even at 30 mM, and then a high error in the kinetic parameters is expected. The increase in the k_cat/K_M values at increasing salt concentrations (either KCl or NaCl) is mainly due to the effect of salt on k_cat.

To explore the biotechnological potential of this halophilic enzyme, we assess the effect of different salts on the enzyme activity under saturating substrate conditions (considering the parameters obtained for NaCl and KCl) (Figure 2). In agreement with the halophilic character of the enzyme, maximum activity was achieved in the presence of KCl or NaCl (over 1.0 M of KCl or NaCl), although the activity is 1.9-fold higher in the presence of 2 M KCl compared to the same NaCl concentration. Other salts, like CsCl, NH₄Cl, and RbCl, also produce an increment in enzyme activity, although the maximum activity is lower than the one obtained in the presence of KCl. On the other hand, LiCl and, to a lesser extent, (NH₄)₂SO₄ produced an inhibitory effect.

Since G6PDH is an enzyme usually employed in coupled assays to measure glucose phosphorylation reactions, we explore the ability of HvG6PDH to be used as an auxiliary enzyme for halophilic hexokinase/glucokinase reactions. As a control, it is necessary to determine the ability of HvG6PDH to oxidize glucose in the presence of NAD⁺ (GlcDH activity). This was assessed by performing glucose saturation curves in the presence of KCl or NaCl (Figure 2). In the presence of NAD⁺, HvG6PDH oxidizes glucose with a linear behavior even up to 500 mM. The catalytic efficiency (k_cat/K_M) of the GlcDH activity in the presence of KCl and NaCl is 2.3 × 10 and 1.7 × 10¹ M^–1s^–1, respectively, values that are more than two orders of magnitude lower than the ones obtained for the G6PDH activity (Table 1).

Based on the mathematical approximation reported by Storer and Cornish-Bowden (Storer and Cornish Bowden, 1974), we design a glucokinase (GK)-coupled assay by employing HvG6PDH as an auxiliary enzyme, considering its GlcDH activity. We calculated the minimum amount of HvG6PDH required to measure an initial velocity of 0.01 mM min^–1 for the GK activity, with a lag time of 1–5 min before the steady state is reached and with a v₂/v₁ ratio of 0.95–0.99 [the v₂/v₁ ratio corresponds to the quotient between the auxiliary enzyme velocity (v₂) and the interest enzyme velocity (v₁) at the steady state] (Supplementary Table 1). Considering these parameters, 2 U/ml of HvG6PDH should be used. This was experimentally tested with the hyperthermophilic ADP-dependent glucokinase from T. litoralis. First, we measured the GlcDH activity in the presence of increasing concentrations of KCl and saturating concentrations of NAD⁺ and glucose. Then, TlGK addition causes a significant increase in the initial velocity, which allows us to calculate the actual glucokinase activity by subtracting the GlcDH velocity (Figure 3). However, the GlcDH activity displayed by HvG6PDH restricts its use as an auxiliary enzyme for GK-coupled assays to conditions where low glucose concentrations are employed (below 1 mM); otherwise, a careful determination of the GlcDH activity must be performed. Despite these limitations, the TlGK hyperthermophilic enzyme remains active up to 2.0 M of KCl, indicating that HvG6PDH can be used as an auxiliary enzyme with the appropriate controls.

We modeled the dimeric structure of HvG6PDH (Figure 4A) using the non-halophilic uronate dehydrogenase from Agrobacterium fabrum (AfUDH, PDB ID: 3RFV) as a template (Figure 4A). The AfUDH is the protein phylogenetically closer to HvG6PDH (29.2% identity) with a determined structure in the presence of ligands. The region between residues 69–80 of HvG6PDH was ab initio modeled using DaReUs-loop (Karami et al., 2019) since it lacks any homology with the template. Moreover, when ligands (G6P and NAD⁺) were added to the model, significant clashes between them and the residues located in the 178–193 region (β6-α10 loop) were found. For this reason, this loop was also ab initio modeled using the loop refine tool of MODELLER. After modeling, we selected the model with the best DOPE score and the higher number of interactions between the model and G6P. The selected model was minimized with ff14SB force field in Amber18. For this minimized structure, the surface electrostatic potential and the solvent accessibility surface of both monomers were calculated. For comparison, these parameters were also calculated for the non-halophilic AfUDH. For HvG6PD, a markedly negative surface potential was observed over the entire protein (Figure 4B), in contrast with the non-halophilic AfUDH (Supplementary Figure 3B). The results obtained agree with those described for haloarchaeal proteins (Britton et al., 2006; Baker et al., 2009) where a high presence of negatively charged residues, such as aspartic and glutamic, are observed (Supplementary Figure 4). However, in both the halophilic as well as the non-halophilic protein, a region of positive electrostatic potential is observed, which corresponds to the substrate binding site, according to the chemical nature of NAD⁺ and G6P, which have bulky negatively charged groups.

The hydrophobic core of HvG6PDH and AfUDH is mainly composed of small hydrophobic residues, like alanine, glycine, leucine, and valine (Figure 4C and Supplementary Figure 3C). Nevertheless, HvG6PDH shows an important diminution of leucine (7.9% HvG6PDH and 14.9% AfUDH) and a slight increase in alanine content (13.4% HvG6PDH and 10.8% AfUDH) when compared to the non-halophilic AfUDH. Besides this, most of the tryptophan residues of the halophilic HvG6PDH are exposed to the solvent (five out of six are exposed), unlike the non-halophilic protein where most of them are located at the hydrophobic core (only 2 out of 10 are exposed). The halophilic protein also presents a reduction in the hydrophobic core, which represents 39.4% of the total volume of the protein, in comparison with the non-halophilic protein where the hydrophobic core represents 45.9%. These results agree with the structural adaptations reported for haloarchaeal proteins (Paul et al., 2008; Graziano and Merlino, 2014).

The minimized model of HvG6PDH shows that the binding of G6P and NAD⁺ at the active site involves residues from the β6-α10 loop in the first coordination shell. Residues H186 and R193 establish polar interactions with G6P, whereas N181 interacts with the di-phosphate of the NAD⁺ molecule and residue H185 with the nicotinamide moiety. The other residues participating in G6P stabilization are N116 and H117 (Figure 4D). Residues located at 5Å from NAD⁺ are mainly glycine and aspartic that interact with the adenosine moiety and residues R13 and R79 that interact with the di-phosphate group. The catalytic residues reported for SDR superfamily members are also found at the modeled active site of HvG6PDH, and correspond to S115, Y152, and K156. These results indicate that HvG6PDH would perform its catalytic function through the conserved mechanism described for SDR superfamily members, and this allows us to propose a role of the aforementioned residues in substrate binding.

To assess the effect of molar concentrations of salt on the structural dynamics of the dimeric structure of HvG6PDH, we performed classical molecular dynamics simulation (cMDS) at a constant temperature of 300 K for 100 ns at low salt concentration (58 atoms or 125 mM of K⁺), only providing the necessary counter-ions to maintain electrolytic neutrality in the system. A second batch of cMDS was also performed, where K⁺ and Cl^– ions (729 atoms of each) were added to simulate a 1.5-molar concentration of salt. The RMSD values are shown to evaluate the conformational integrity and the global changes of the HvG6PDH model at low and high salt concentrations. The average (±SD between replicates) RMSD values of three replicates at low KCl reached 4.08 ± 0.19 Å in the last 10 ns, whereas at 1.5 M KCl, this value is only 3.30 ± 0.14 Å (Figure 5A). To determine the basis of this difference, we separately analyzed the R_g of each monomer of the protein in both conditions (low and high salt). The average value during the trajectory is different for both monomers; in the case of monomer 1, R_g is 17.8 Å at high and 18.8 Å at low KCl concentration. Nonetheless, monomer 2 maintains a very similar R_g average in both conditions over all the trajectories (18.3 Å at high KCl and 18.2 Å at low KCl) (Figure 5B). This result suggests a role for salt concentration in the conformational dynamics of HvG6PDH; at low salt concentration, one of the monomers is more swollen, which could be the result of partial unfolding of the protein or a higher conformational dynamic. Then, we analyzed the atomic fluctuation of our model, evaluating each monomer as an independent chain. At low KCl concentration, the overall mean RMSF is 1.36 ± 0.30 Å (monomer 1) and 1.32 ± 0.28 Å (monomer 2). In contrast, at high KCl concentration, the RMSF values are lower than those obtained at low salt, 1.03 ± 0.19 and 1.08 ± 0.19 Å for monomers 1 and 2, respectively. A more detailed analysis of each monomer in both conditions revealed regions with high RMSF variations in monomer 1 (Figures 5C,D). At low salt concentration, regions such as the α2 helix, the β6-α10 loop, and the segment from α11 helix to β7 sheet of monomer 1 present very high RMSF average values, with these regions being stabilized at high salt concentration (Figures 5E,F).

Moreover, a close inspection of the substrates located at the active site shows that they do not maintain a catalytic conformation in replicates of either condition. To analyze the enzyme–substrate interactions, a 1-kcal/mol Ų restriction to G6P, NAD⁺, and the catalytic tyrosine (Y152) was applied, and a new set (3) of molecular dynamics was performed at low and high salt concentrations. In these new systems, the substrates maintain a catalytic conformation, and the potential energy of non-bonded interactions as well as the electrostatic and the van der Waals component between G6P and the protein were evaluated. The average values for the electrostatic interactions at low KCl were −196.09 ± 21.84 and −230.77 ± 14.54 kcal/mol Ų for monomers 1 and 2, respectively. On the other hand, these values for the van der Waals interactions were −26.37 ± 3.45 kcal/mol Ų for monomer 1 and −15.95 ± 2.86 kcal/mol Ų for monomer 2 (Supplementary Figure 5A). These results indicate that attractive electrostatic interactions are the main energetic component in the interaction between G6P and the protein. At high KCl concentration, the average value for the electrostatic component in monomers 1 and 2 are −192.85 ± 51.73 and −114.50 ± 35.20 kcal/mol Ų, while for the van der Waals interactions these are −16.29 ± 3.13 and −19.45 ± 3.49 kcal/mol Ų, respectively (Supplementary Figure 5B). The variation observed in the mean values of the electrostatic component could be attributed to a shielding effect due to the high ion concentration, although the electrostatic interactions of the protein with G6P were still attractive.

Interestingly, the β6-α10 loop, which was modeled ab initio, shows a significant change in the molecular trajectories in the presence of low or high salt. At low salt concentration, the electrostatic component of the interactions between this loop and G6P is mainly repulsive, being +17.90 ± 5.05 and +25.74 ± 5.16 kcal/mol Ų for monomers 1 and 2, respectively (Figure 6). Surprisingly, in all the replicates performed at high salt concentrations, the electrostatic component became attractive, being −19.58 ± 17.19 kcal/mol Ų for the first monomer. However, for the second monomer, the average value diminished from 22.49 ± 4.92 in the first 25 ns of the simulation to −14.92 ± 7.66 in the last 25 ns, changing from a repulsive to an attractive interaction (Figure 6).

The evolutionary relationship of archaeal H. volcanii G6PDH with SDR superfamily members was analyzed by molecular phylogeny using a PSI-Blast search. Sequences of characterized SDR enzymes, like UDP-glucose 4,6 epimerase, GDP-L-fucose synthase, L-arabinose dehydrogenase, and uronate dehydrogenases from eukaryotes, bacterial, and archaeal lineages, were included in the alignment. We identified 11 monophyletic groups (Figure 7A) where HvG6PDH appears in a group of uncharacterized sequences from Halobacteria organisms, including the orders Natrialbales, Haloferacales, and Halobacteriales (Figure 7A, group 1, and Figure 7B). Members of this group have the conserved NLTX₂H motif of the β6-α10 loop (sequence logo in Figure 7B) identified as relevant for the interaction with G6P in molecular dynamic simulations. A second conserved motif of this group also corresponds to the YERG motif, where R193 could be important for phosphate binding (Figure 7C). Close to group 1, another uncharacterized group of Halobacteria sequences appear (Figure 7A, group 2), whose sequences show an insertion inside the β6-α10 loop, which could indicate a different substrate specificity (Figure 7C). A branch composed exclusively of bacterial sequences including Chloroflexi, Proteobacteria, and Nitrospira organisms (group 5) was also identified. The most populated group (group 6) comprises enzymes with activities like dehydrogenases, epimerase–dehydratase, and fucose synthase from the three domains of life. The other four groups are composed exclusively by sequences from different lineages of bacteria with unknown enzymatic function, except for group 10, where the sequence of uronate dehydrogenase from A. fabrum is located.