Probing the Effect of Glucose on the Activity and Stability of β-Glucosidase: An All-Atom Molecular Dynamics Simulation Investigation

β-Glucosidase (EC 3.2.1.21) plays an essential role in the removal of glycosyl residues from disaccharide cellobiose to produce glucose during the hydrolysis of lignocellulosic biomass. Although there exist a few β-glucosidase that are tolerant to large concentrations of glucose, these enzymes are typically prone to glucose inhibition. Understanding the basis of this inhibition is important for the production of cheaper biofuels from lignocellulose. In this study, all-atom molecular dynamics simulation at different temperatures and glucose concentrations was used to understand the molecular basis of glucose inhibition of GH1 β-glucosidase (B8CYA8) from Halothermothrix orenii. Our results show that glucose induces a broadening of the active site tunnel through residues lining the tunnel and facilitates the accumulation of glucose. In particular, we observed that glucose accumulates at the tunnel entrance and near the catalytic sites to block substrate accessibility and inhibit enzyme activity. The reduction of enzyme activity was also confirmed experimentally through specific activity measurements in the presence of 0–2.5 M glucose. We also show that the increase in glucose concentrations leads to a decrease in the number of water molecules inside the tunnel to affect substrate hydrolysis. Overall, the results help in understanding the role of residues along the active site tunnel for the engineering of glucose-tolerant β-glucosidase.


INTRODUCTION
The production of biofuels from lignocellulosic biomass could be broadly divided into three stepspretreatment of the biomass to reduce cellulose recalcitrance, enzymatic conversion of the polysaccharides in pretreated biomass to sugars, and the conversion of the sugars to fuels by the same or other microbes through fermentative processes. A group of enzymes, known as cellulase, break down the polysaccharides into sugars. 1,2 These enzymes work synergistically to hydrolyze the biomass and are made up of the minimum set of endoglucanase, cellobiohydrolase, and β-glucosidase. Endoglucanase (EC 3.2.1.4) randomly cleave the β-1,4 glycosidic linkages of cellulose, cellobiohydrolase (EC 3.2.1.91) attack the cellulose chain ends to produce a dimer of glucose called cellobiose, which is linked by a β-1,4 glycosidic bond. β-Glucosidase (EC 3.2.1.21) hydrolyze the β-1,4 linkage in cellobiose to produce glucose. 3−6 The final hydrolysis step is generally recognized as a limiting step in the conversion of lignocellulosic biomass to sugars as β-glucosidase are inhibited by its reaction product glucose. 4,6,7 The efficient and economic hydrolysis by cellulase requires enzymes that are also active in high glucose concentrations. Understanding the basis of glucose inhibition and engineering toward relieving this product inhibition is, therefore, a major challenge.
Whereas the typical inhibition constant (K i,app ) of glucose on the β-glucosidase chromogenic substrate, p-nitrophenyl-β-Dglucoside (pNPGlc), is less than 50 mM, recently there have been reports of bacterial β-glucosidase that are either not inhibited by glucose or are stimulated in the presence of glucose. 6,8,9 Thus, there exists a large difference in glucose inhibition amongst β-glucosidase. Though the molecular mechanism of glucose inhibition or tolerance is not well understood, experimental studies of glucose-tolerant βglucosidase have attributed tolerance to narrower and deeper active site tunnels, specific residues in the tunnel, and the presence of allosteric sites on these enzymes. 10−12 Molecular dynamics (MD) simulation is an important tool to characterize structural dynamics, flexibility, and mechanism of enzymes and help gain deeper insights into the structure−function relationship. 13−16 Computational modeling and MD simulation have been reported toward understanding catalytic activity and thermostability of cellulases. 15,17−23 However, there are no reports of theoretical studies for a molecular-level understanding of β-glucosidase inhibition or the effect of glucose on the structural stability of such enzymes.
As part of a program to understand how the β-glucosidase reaction product, glucose, interacts with the enzyme, we have been studying the effects of glucose on β-glucosidase across mesophilic and thermophilic organisms. 6,9,24 One such enzyme is the highly active and thermostable GH1 β-glucosidase (B8CYA8) from Halothermothrix orenii. 24,25 The higher stability of this enzyme has been shown to be a promising target toward industrial applications. 26 In the present work, we report the effect of glucose on B8CYA8 activity and employed MD simulations to investigate the influence of glucose on B8CYA8 activity and stability. In order to understand the effect of different concentrations of glucose on protein structure and activity, we computed the root mean square deviation (rmsd) of amino acids, root mean square fluctuations (RMSFs) of C α atoms, radial distribution functions, and distance distribution between the selected pair at different temperatures. The MD studies provide us a way to understand the basis of B8CYA8 glucose inhibition of catalytic activity at high glucose concentrations.

METHODS
2.1. Protein Expression, Purification, and Assay. The synthetic gene corresponding to the β-glucosidase from H. orenii was constructed (BankIt1930137BG_Halotherm KU867899) as reported previously and expressed in Escherichia coli Top 10F′ cells (Life Technologies, La Jolla, CA). 24 B8CYA8 was purified using a protocol detailed earlier and purity was confirmed by 10% sodium dodecyl sulfate polyacrylamide gel electrophoresis. 24 The specific activity of enzymes in the presence of glucose was determined on the substrate p-nitrophenyl β-D-glucopyranoside (pNPGlc) after incubation in the buffer containing different concentrations of glucose under optimum conditions.
2.2. Computational Methodology and Simulation Details. The X-ray crystallographic structure of the βglucosidase was obtained from the protein data bank (PDB ID: 4PTX). 25 At first, the energy-minimized protein molecule (B8CYA8) was kept at the center of a cubic simulation box 120 Å long and then solvated with water. We used the TIP3P water model in all of our simulations. 27 We added glucose molecules by using PACKMOL 28 to obtain the following desired glucose concentration, X glu of 0, 0.2, 0.5, 1.5, and 2.0 M. All the potential parameters were obtained from the CHARMM36 force field. The simulations were performed by NAMD-2.9 simulation tools. 29−32 By taking initial configurations at different X glu , the conjugate gradient method was applied for 300 000 steps to remove all the energetically unfavorable contacts. Three different temperatures were chosen to represent the room temperature (300 K), high catalytic activity (335 K), and maximum temperatures at which any catalytic activity could be measured (375 K). We have also performed simulation at 450 K at which the protein is denatured (450 K). The starting configurations were equilibrated for 3.0 ns in the NPT ensemble to fix the simulation box length. In NPT simulations, Nose−Hoover thermostat and barostat coupling constants were taken to be 0.5 and 2.0 ps, respectively. 33,34 The pressure was kept constant at 1.0 atm. After the simulation box length was fixed, we equilibrated the system for 5.0 ns in the NVT ensemble. To analyze different properties of the system, a production run of 50 ns was performed in the NVT ensemble for all cases. At 335 K, a production run of 100 ns was also performed for all glucose concentrations and at 300 and 375 K, a longer production run was performed for 0.5, 1.5, and 2.0 M glucose concentrations to check the validity of the production run up to 50 ns. We observed that there were no qualitative differences in the results for 50 and 100 ns production runs. The temperature of the systems was kept constant at the respective temperature by using the damping coefficient (γ) of 1.0 ps −1 by Langevin dynamics. Long-range interactions are handled by the particle mesh Ewald method, 35−37 with real space cut-off of 16 and 2 Å pair-list cut-off. We used a scaling factor of 1−4 in all our simulations. The time step was of 1.0 fs and all the properties were computed from the trajectories stored at an interval of 4.0 ps during the production run.

RESULTS AND DISCUSSION
3.1. Validation of the Forcefield. We determined the rmsd of the energy-minimized B8CYA8 backbone at 300, 335, 375, and 475 K in the absence of glucose. Figure 1a shows that up to 50 ns, rmsd is less than 2.5 Å for temperatures up to 375 K. However, above 375 K, rmsd increased drastically, suggesting a temperature-induced denaturation of B8CYA8. This agreed well with the experimentally determined unfolding temperature of this protein 24 and proved that the force field used in this study was capable of reproducing the experimentally observed stability of B8CYA8. Figure 1b shows the root-mean-square fluctuation (RMSF) of each amino acid residue at different temperatures. Up to 375 K, the residues have a maximum RMSF of less than 3.5 Å. However, at 475 K the RMSF of many B8CYA8 residues are more than 5.0 Å and indicates protein denaturation at that temperature. 38 As the experimental data indicated that B8CYA8 displays maximum activity at 335 K and denatures at temperature close to 375 K, we will hereafter restrict our discussions to enzyme behavior at 300, 335, and 375 K.
3.2. Structural Changes of B8CYA8 in the Presence of Glucose. The B8CYA8 crystal structure indicates the presence

ACS Omega
Article of a TIM barrel with a (α/β) 8 fold and the active site located deep inside a tunnel. 27 This enzyme active site tunnel is made up of catalytic residues E354 and E166 located at the glyconebinding site at the bottom of the tunnel, aglycone-binding site located toward the middle of the tunnel, and gatekeeper residues at the edge of the active site tunnel near the tunnel opening. The aglycone-binding site contains a diverse combination of residues and generally forms hydrophobic interactions as well as water-mediated H-bonds with the substrate. 39 The glycone-binding sites contain a highly conserved combination of residues that can also form Hbonds with the substrate. Thus, the amino acids in the active site tunnel are thought to play an important role in βglucosidase activity and glucose tolerance. 39−44 The B8CYA8 glycone, aglycone, and gatekeeper residues are listed in Table  S1.
In order to understand any structural changes caused by the reaction product glucose at different regions of the enzyme including the active site tunnel, we started by computing the rmsd of the B8CYA8 backbone in the presence of glucose. Figure 2a−c shows rmsd at 300, 335, and 375 K, respectively, in the presence of different concentrations of glucose. At 300 K, the rmsd is nearly independent of glucose concentrations. However, at 335 and 375 K subtle rmsd variations were observed across different glucose concentrations. To better understand the gross structural changes, we computed average rmsd (rmsd avg ) from the last 20 ns production runs at different temperatures and glucose concentrations ( Table 1). The columns under temperatures in Table 1 represent average rmsd. At 335 and 375 K, the rmsd avg of B8CYA8 increases at high glucose concentrations with a maximum rmsd avg of 3.51 Å. This suggests structural changes in B8CYA8 up to 375 K, when the protein remains folded. To visualize these structural changes, we took the last configuration of the enzyme from each simulation and aligned the backbone C α atoms. Figure  S1A−E shows the last B8CYA8 structure from each MD run at 335 K and Figure 3 shows their aligned structures.
Significant structural changes at residues 324, 325, 410, and 411 in the gatekeeper region and residue 415 at the glyconebinding site in the active site tunnel were observed and are shown in the region within black dotted circles in Figure 3. It has been previously reported that residues in the glyconebinding and gatekeeper region may play an important role in the activity of B8CYA8. 27 Thus, the residues causing structural changes in the presence of glucose may also have a significant effect on B8CYA8 activity.
To further pinpoint the structural changes of all the residues, we first computed RMSF as a function of glucose concentration (X glu ). Figure 4a−c shows RMSF at T = 300, 335, and 375 K, respectively, in the presence of 0−2.0 M glucose. RMSF is mostly invariant with glucose concentration, except for the residues in the active site tunnel region. As structural changes around the catalytic sites (E354 and E166) are expected to play an important role in regulating protein activity, we zoomed Figure 4a−c to highlight the RMSF of residues between 160−170 and residues 350−360 containing the catalytic residues E166 and E354, respectively. These regions are shown in Figure 4d−f,g−i, respectively, with the RMSF measured at three different temperatures. We observed small changes in RMSF near the catalytic site residues in the presence of glucose. Previous reports in the literature implicate specific regions around the active site tunnel like a loop region, gatekeeper region, glycone-binding site residues, and aglyconebinding site to play an important role in regulating protein activity. 27

ACS Omega
Article and lower enzyme activity. Indeed, the reduction of enzyme activity was confirmed upon experimental determination of B8CYA8 specific activity in the presence of 0−2.5 M glucose ( Figure 6). Ongoing experimental studies have confirmed some of these results. For example, W327R mutant activity decreases sixfold compared to the wild-type enzyme (data not shown).
To further confirm that the fluctuations of the abovementioned amino acid residues are due to glucose, we computed the radial distribution function [g(r)], where r represents the distance between glucose and the solventexposed amino acid residue of the protein that can interact with glucose. We chose those solvent-exposed residues based on the maximum change in RMSF with X glu and compared with the residues where there were no changes in the presence of glucose. At 300 K, we chose D307 and at 335 and 375 K, we chose E318 as these residues showed the maximum RMSF. We chose L309 as RMSF for this residue did not vary with glucose concentration. We measured g(r) between OH 1 (see Figure  S2a) of glucose, which can form hydrogen bonds with the protein, and Oδ 1 , Oε 1 , and Hδ 12 of D307, E318, and L309, respectively (atom notations are defined in Figure S2b−d). Figure 7a−d,e−h,i−l shows g(r) at 300, 335, and 375 K, respectively, in the presence of different concentrations of glucose. At large distances, g(r) at different glucose concentrations do not tend to 1 (one) because the bulk density of glucose is different. At all temperatures when r ≤ 5.0 Å, we observed multiple distinct peaks in Oδ 1 −OH 1 and Oε 1 − OH 1 g(r), whereas no such peak was observed for Hδ 12 −OH 1 g(r). The peak heights indicate that glucose accumulation near D307 and E318 is large as compared to L309. Other such chosen pairs also show similar results and confirm that the fluctuations of amino acid residues are indeed due to the presence of glucose. The relatively larger peak heights suggest that the amino acids with large RMSF interacts favorably with

ACS Omega
Article glucose. Moreover, glucose approaches closer to D307 and E318 (r ≤ 2.4 Å) than L309 at all temperatures and concentrations studied, suggesting the possibility of hydrogen bond interaction with D307 and E318. As the glucose concentration increases from 0.2 to 2.0 M, the peak height in g(r) mostly decreases because of the finite number of glucose molecules that can be accommodated within the steric limits imposed by D307 and E318 in the active site tunnel.
3.2.1. Broadening of the Active Site Tunnel Because of Glucose. Enzyme activity typically depends on the concentration of solvent and substrate molecules near the active site residues. As the catalytic site in B8CYA8 is located at the bottom of the active site tunnel, the sterics of the tunnel entrance may affect the accessibility of substrate and solvent at the catalytic site. To measure any changes in the diameter of the active site tunnel, we measured the distances between different pairs of residues such as N308−A410, D307−A410, Y411−K316, and A410−V314 located nearly opposite to each other at the gatekeeper region near the tunnel entrance [ Figure  8a]. We computed the distance distribution [P(R)], where R is the distance between the center of mass (c. o. m.) of the chosen residue pairs. The P(R) at 300, 335, and 375 K for the four pairs of residues are shown in Figures S3−S5, respectively. R varied up to 7.0 Å at 300, 335, and 375 K. Though R varies randomly with glucose concentrations at 300 K, the R increases with glucose concentrations at 335 and 375 K. In order to obtain the tunnel diameter, we computed the average distribution of the four chosen pairs of residues, as shown in Figures

ACS Omega
Article other glucose concentrations because of large fluctuations of the gatekeeper residues.
3.2.2. Accumulation of Glucose at the Entrance and inside the Active Site Tunnel. The experimentally observed inhibition of B8CYA8 activity by glucose [shown in Figure 6] may be caused by the accumulation of glucose at the gatekeeper region near the tunnel entrance and/or near the active site. To delineate between these two possibilities, we measured the number of glucose molecules (N glucose ) as a function of distance from the catalytic sites [r d ] at all the glucose concentrations. To compute N glucose , we considered three hydrogen atoms of glucose, namely HO 3 , HO 4 , and HO 6 and oxygen atom (Oε 1 ) of the catalytic site residues as shown in Figure 9. Figure 10 shows N glucose as a function of r d at 335 K at different glucose concentrations. As the approximate length of the tunnel is 19.0 Å, we restricted r d up to 20 Å. r d = 0 Å represents the location of Oε 1 atom of the two catalytic residues. We observed that at the gatekeeper region, N glucose increases significantly with glucose concentration, in particular at X glu = 1.5 and 2.0 M. Thus, the accumulation of glucose at the gatekeeper region probably hinders substrate entry into the tunnel. We also observed that the distance of closest approach of glucose to the catalytic residues decreases at high glucose concentrations. For example, at X glu = 2.0 M the distance between Oε 1 and hydrogen atoms of glucose is nearly 2.5 Å whereas at 0.2 M, the distance increases to around 7.5 Å. The increased presence of glucose near the active sites at high glucose concentration competes with substrate binding, hinders substrate accessibility, and inhibits B8CYA8 activity. At 300 and 375 K (shown in Figures S6 and S7), a similar trend was observed. Our studies therefore show that the presence of a high concentration of glucose induces tunnel broadening and glucose accumulation both at the catalytic site and the tunnel entrance.
3.3. Proton-Donating Capacity of the Active Site Residues in the Presence of Glucose. It is known that the catalytic site E166 acts as a proton donor and another catalytic residue E354 acts as a nucleophile. 27 GH1 β glucosidase are typically hydrolyzed through the well-known retaining mechanism where E166 is first protonated by neighboring water and then acts as a proton donor. 39 The proton-donating capacity of the catalytic residues should depend on the aqueous environment adjacent to the catalytic site. We, therefore, computed the average number of water molecules (N water ) in the tunnel within 3.5 Å from the Oε 1 atom of E166. Figure 11 shows N water at 335 K for different glucose concentrations. It is seen that the number of water molecules adjacent to E166 decreases significantly. The reduction of N water near E166 with glucose concentration may be due to increased accumulation of glucose near the catalytic site (see Figure 10) at a high glucose concentration. The fewer number of water molecules at high glucose concentrations might decrease the proton-donating ability of the catalytic site (E166) and/or affect the release of glucose from the intermediate state by the nucleophilic water and lead to decrease in specific activity of B8CYA8.

CONCLUSIONS
MD simulations at different temperatures and glucose concentrations assist us in understanding the effect of glucose on B8CYA8 activity. At high glucose concentrations, we observed significant fluctuations at residues corresponding to the loop region (303−308), gatekeeper residues (314, 316, 324, 325, 326, 410, and 411), aglycone-binding site residue (327), and glycone-binding site residues (401, 406, 408, 409, 415, and 417). The experimentally observed inhibition of B8CYA8 activity in the presence of increasing concentrations of glucose may be due to the effect of these fluctuations on substrate binding and enzyme activity. We found an increase in glucose accumulation at the tunnel entrance and near the active site which may compete with substrate accessibility and binding to the catalytic site. At 335 K when the enzyme shows

ACS Omega
Article high activity, glucose-induced tunnel broadening was observed that may further facilitate the accumulation of glucose at the tunnel entrance. Furthermore, in the vicinity of catalytic residues, the decrease in the number of water molecules in the presence of glucose could reduce proton donation to the substrate and inhibit enzyme hydrolytic activity. Our results implicate the role of the active site tunnel in maintaining glucose tolerance of B8CYA8 and argues against the role of any single residue. These results may guide a closer examination of the β-glucosidase active site tunnel toward understanding the role of protein−solvent interaction in the regulation of polysaccharide hydrolysis and engineering of improved enzymes.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsomega.9b00509.
Glycone, aglycone, and gate-keeper residues of B8CYA8; B8CYA8 structures for different glucose concentrations at 335 K; atom notation for glucose, D307, E318, and L309 amino acids; distance distribution [P(R)] for chosen pair of amino acids at 300, 335, and 375 K; and numbers of glucose (N glucose ) as a function of the distance (r d ) between different hydroxyl groups of glucose (OH 3 , OH 4 , OH 6 ) and the Oε 1 of E166 at