Computational Insights into the Effects of pH on β-glucosidase Stability and Catalysis

The catalytic activity of β-glucosidase from Neotermes koshunensis (BgluN) exhibits a very steep pH dependence, indicating the complexity of the active site, where ligand binding is intricately coupled to the protonation equilibria. To clarify this connection, we performed constant-pH molecular dynamics (CpHMD) simulations over a pH range of 4–7 for the apo, substrates-bound, and product-bound forms of BgluN. The protonation profiles were suggestive that the E193 (acid/base), E402 (nucleophile), and E451 (auxiliary) residues are the ones modulating the pH-dependent catalytic activity. Within the optimal catalytic range (pH 5.0–5.5), there is a proton shared between these residues that is needed for catalysis. At higher pH, the overall site becomes deprotonated, compromising the acid/base role in the catalytic cycle. Binding affinity analyses using MM/PBSA and Wyman–Tanford Linkage methods revealed that the ligand binding is most thermodynamically favorable at pH 5.0, in good agreement with the experimentally determined activity optimum. Overall, these findings offer a mechanistic framework that integrates the dynamics of protonation, binding energetics, and catalytic readiness. Additionally, it identifies the nearby negatively charged E451 residue as a potential target for mutagenesis to fine-tune BgluN variants with shifted pH optima, thereby making them more suitable for industrial biocatalysis.