ORIGINAL RESEARCH published: 10 April 2019 doi: 10.3389/fchem.2019.00200 Comparing Hydrolysis and Transglycosylation Reactions Catalyzed by Thermus thermophilus β-Glycosidase. A Combined MD and QM/MM Study Sonia Romero-Téllez 1,2 , José M. Lluch 1,2 , Àngels González-Lafont 1,2* and Laura Masgrau 1,2* 1 Departament de Química, Universitat Autònoma de Barcelona, Barcelona, Spain, 2 Institut de Biotecnologia i de Biomedicina, Universitat Autònoma de Barcelona, Barcelona, Spain The synthesis of oligosaccharides and other carbohydrate derivatives is of relevance Edited by: for the advancement of glycosciences both at the fundamental and applied level. For Vicent Moliner, University of Jaume I, Spain many years, glycosyl hydrolases (GHs) have been explored to catalyze the synthesis Reviewed by: of glycosidic bonds. In particular, retaining GHs can catalyze a transglycosylation (T) Alessandro Silva Nascimento, reaction that competes with hydrolysis (H). This has been done either employing University of São Paulo, Brazil Xabier Lopez, controlled conditions in wild type GHs or by engineering new mutants. The goal, which University of the Basque Country, is to increase the T/H ratio, has been achieved with moderate success in several Spain cases despite the fact that the molecular basis for T/H modulation are unclear. Here *Correspondence: we have used QM(DFT)/MM calculations to compare the glycosylation, hydrolysis and Àngels González-Lafont angels.gonzalez@uab.cat transglycosylation steps catalyzed by wild type Thermus thermophilus β-glycosidase Laura Masgrau (family GH1), a retaining glycosyl hydrolase for which a transglycosylation yield of 36% has laura.masgrau@uab.cat been determined experimentally. The three transition states have a strong oxocarbenium Specialty section: character and ring conformations between 4 H3 and 4 E. The atomic charges at the This article was submitted to transition states for hydrolysis and transglycosylation are very similar, except for the more Theoretical and Computational negative charge of the oxygen atom of water when compared to that of the acceptor Glc. Chemistry, a section of the journal The glycosylation transition state has a stronger SN 2 character than the deglycosylation Frontiers in Chemistry ones and the proton transfer is less advanced. At the QM(PBE0/TZVP)/MM level, the Received: 20 December 2018 TS for transglycosylation has shorter O4GLC -C1FUC (forming bond) distance and longer Accepted: 15 March 2019 Published: 10 April 2019 OE2GLU338 -C1FUC (breaking) distance than the hydrolysis one, although the HACC proton Citation: is closer to the Glu164 base in the hydrolysis TS. The QM(SCC-DFTB)/MM free energy Romero-Téllez S, Lluch JM, maxima show the inverted situation, although the hydrolysis TS presents significant González-Lafont À and Masgrau L structural fluctuations. The 3-OHGLC group of the acceptor Glc (transglycosylation) (2019) Comparing Hydrolysis and Transglycosylation Reactions and WAT432 (neighbor water in hydrolysis) are identified to stabilize the oxocarbenium Catalyzed by Thermus thermophilus transition states through interaction with O5FUC and O4FUC . The analysis of interaction β-Glycosidase. a Combined MD and QM/MM Study. Front. Chem. 7:200. suggests that perturbing the Glu392-Fuc interaction could increase the T/H ratio, either doi: 10.3389/fchem.2019.00200 by direct mutation of this residue or indirectly as reported experimentally in the Asn390I Frontiers in Chemistry | www.frontiersin.org 1 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation and Phe401S cases. The molecular understanding of similarities and differences between hydrolysis and transglycosylation steps may be of help in the design of new biocatalysts for glycan synthesis. Keywords: glycosylhydrolase, transglycosylation, QM/MM, glycans, hydrolysis, glycosyl-enzyme complex, GH1 INTRODUCTION steps with formation of a covalent glycosyl-enzyme intermediate. Two catalytic carboxyl groups (separated by ∼5Å) act as general Carbohydrates and their derivatives (glycoconjugates and acid/base and nucleophile in the reaction (as exemplified in glycosides) play important biological functions, including energy Figure 1 for the process studied in this work). In a first step storage, structural roles and also the encoding of a molecular (referred as glycosylation, G), the nucleophile attacks the and cell recognition language that drives the specificity of anomeric carbon while the carboxylic acid protonates the such processes. Accordingly, glycans are vital for normal life glycosidic oxygen of the leaving group, and the glycosyl-enzyme development and are also key factors in the progression of intermediate is formed. In the final step (deglycosylation, D), many diseases (from pathogen infection to cancer) (Seeberger the residue previously acting as an acid now deprotonates and Cummings, 2017). In order to develop fundamental the incoming acceptor substrate that attacks the anomeric and applied research in the glycosciences field, access to carbon forming the final product with net retention of the pure glycans and in sufficient amounts is needed. Glycans configuration. If the acceptor is a water molecule, hydrolysis biosynthesis involves the action of a repertoire of enzymes, (H) occurs. However, with the presence of a different suitable amongst them glycosyltransferases (GTs, which catalyze the sugar acceptor, many retaining GHs are capable of catalyzing synthesis of a new glycosidic bond) and glycoside hydrolases transglycosylation (T). (GHs). These enzymes catalyze the corresponding reaction with Different approaches and experimental conditions have been two possible stereochemical outcomes, that is, retention or used to enhance the transglycosylation vs. hydrolysis (T/H) inversion of the configuration at the anomeric carbon of the ratio, e.g., high acceptor concentrations, the use of highly transferred/hydrolyzed sugar. reactive glycosyl donors (activated donors) like aryl glycosides Besides being the target of drug design investigations, or glycosyl fluorides, removal of the transglycosylation product carbohydrate-active enzymes are of interest in glycans processing or enzyme immobilization. With all, yields and regiospecificities and synthesis, an important area of fundamental research and are still low (Planas and Faijes, 2002). Engineering of GHs also for the preparation of commercially-valuable products. has also been investigated. A major breakthrough in the field Glycoside hydrolytic enzymes (e.g., amylases, cellulases) are was accomplished with the development of glycosynthases, used in the e.g., food, textile, detergents, cosmetics, pulp, and in which hydrolysis is abolished by replacing the catalytic paper industries, representing around one third of the global nucleophile by a non-nucleophilic residue and the use of industrial enzyme market (Plou et al., 2007; Sajith et al., 2016). activated substrates (Mackenzie et al., 1998; Malet and Planas, The enzymatic synthesis of glycosidic bonds by carbohydrate- 1998; Moracci et al., 1998; Faijes and Planas, 2007; Cobucci- active enzymes has also been a subject of study for over 60 years Ponzano and Moracci, 2012; Danby and Withers, 2016). Other and different strategies have been developed that circumvent studies have tried to understand the molecular determinants for or complement chemical approaches, which usually require having a transglycosylation vs. an hydrolytic enzyme activity multiple protection, and deprotection steps to obtain the desired by applying mutational studies (Bissaro et al., 2015b). Notably, oligosaccharide (Kiessling and Splain, 2010; Bissaro et al., 2015b; Dion and coworkers (Feng et al., 2005) successfully applied Danby and Withers, 2016). In Nature, GTs are the main catalysts directed evolution techniques to the β-glycosidase of Thermus for the synthesis of glycosidic bonds. However, their broad thermophilus (Ttβ-gly), belonging to the CAZy family GH1 application is hindered by the difficulties in their expression (Henrissat, 1991; Cantarel et al., 2009; Hart and Copeland, and purification and due to the economic costs of obtaining the 2010), to increase its ability to synthesize oligosaccharides nucleotide-phosphate sugars they use as donor substrates, despite by transglycosylation. progresses have been made in that direction (Field, 2011). As an Their findings were later “semi-rationalized” by realizing that alternative, GHs have been explored to catalyze glycosidic bond formation (Planas and Faijes, 2002; Cobucci-Ponzano et al., 2011; mutations leading to an increase in the T/H ratio in Ttβ-gly were Bissaro et al., 2015b). GHs are more abundant than GTs, are much located at highly conserved positions at the−1 subsite (the one easier to obtain, cover a wide range of substrate specificities and accommodating the hydrolyzed/transferred monosaccharide) their substrates are cheap. (Teze et al., 2014). In this way, they created Ttβ-gly mutants In particular, retaining GHs can operate in a synthetic (e.g., Tyr284Phe, Asn282Thr, Phe401Ser, Asn163Ala, Arg75Ala) mode if the equilibrium is displaced toward glycoside bond that catalyzed transglycosylation between 4-nitrophenyl β-D- formation (thermodynamically controlled processes) or fucopyranoside (pNP-Fuc) (donor substrate) and N-methyl- by using activated glycosyl donors (kinetically controlled O-benzyl-N-(β-D-glucopyranosyl)-hydroxylamine (BnON(Me)- transglycosylation) (Planas and Faijes, 2002). Retaining GHs Glc, acceptor), with increased yields of up to 82% against a 36% follow a double displacement mechanism in two subsequent in wild type enzyme, although in all cases the catalytic efficiency Frontiers in Chemistry | www.frontiersin.org 2 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation FIGURE 1 | Catalytic mechanism of the glycosylation and deglycosylation (hydrolysis or transglycosylation) steps in retaining β-glycoside hydrolases. In this work, nucleophile and acid/base residues are GLU 338 and GLU 164, respectively, the donor substrate is pNP-Fuc and the acceptor substrate for transglycosylation is BnON(Me)-Glc. is significantly reduced when compared to WT Ttβ-gly. The product release and noticed the importance of aromatic rings in approach has since then been applied to other GHs families the aglycone site to facilitate acceptor (other than water) binding. (glycosyl hydrolases have been classified into different families Compiling evidences support the idea that the properties of the according to sequence similarity) (Henrissat, 1991; Cantarel deglycosylation transition state, substrate-specific interactions et al., 2009; Teze et al., 2015; Saumonneau et al., 2016). involving the acceptor substrate and specific features to channel The authors postulated that mutation of these well-conserved and retain water molecules close to the catalytic center, are key residues around the−1 subsite may induce lower stabilization of determinants for the T/H partition (Kempton and Withers, 1992; the deglycosylation transition states (TS), being the effect larger Bissaro et al., 2015a,b; David et al., 2017). Still, the molecular for hydrolysis than for transglycosylation reaction, therefore details for the basis of T/H modulation are unclear. resulting in the higher T/H ratio measured experimentally (Teze Computational studies on retaining GHs have mainly been et al., 2014). This suggests that the H and T transition states focused on the glycosylation step, fewer to the hydrolysis one present different characteristics and stabilization interactions. and very few to transglycosylation. For enzymes belonging For another family GH1 β-glucosidase (NkBgl), mutation of to CAZy family GH1, the one investigated in this work, the catalytic acid/base glutamate to aspartate was also found to quantum mechanical molecular mechanical (QM/MM) studies generate new transglycosylation products not produced by the on glycosylation and hydrolysis steps catalyzed by Oryza sativa wild-type enzyme (Jeng et al., 2012). The authors also proposed (rice) β-glucosidase (Osβ-gly) acting on glucose disaccharides that this catalytic residue is important in substrate entry and have been reported (Wang et al., 2011, 2013; Badieyan et al., Frontiers in Chemistry | www.frontiersin.org 3 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation 2012). They corroborated the predicted oxocarbenium-like Here we present a QM(DFT)/MM study of the glycosylation, nature of both TSs, which present the corresponding breaking hydrolysis and transglycosylation steps catalyzed by wild- bond already broken and the forming one far from established. type Ttβ-gly using pNP-Fuc as donor and BnON(Me)-Glc as The interaction of the C2 hydroxyl group from the sugar with the acceptor. Ttβ-gly catalyzes the hydrolysis of β-D-galactoside, catalytic residues was also confirmed to play an important role in β-D-glucoside, and β-D-fucoside derivatives, showing the the catalytic reaction by facilitating the sugar ring conformational highest activity with the latter (Dion et al., 1999). As change toward the TS and by contributing to TS stabilization. mentioned, in the presence of a suitable acceptor, it can This is in agreement with experimental data that have for long also catalyze transglycosylation; BnON(Me)-Glc has been identified this interaction as important for TS stabilization, reported to produce a high transglycosylation yield in the with a contribution of up to 10 kcal/mol (Zechel and Withers, native enzyme and a sole regioisomer product (N-methyl-O- 2000). Recently, the effect of such interaction on lowering the benzyl-N-(β-D-fucopyranosyl(1→ 4)β-D-glucopyranosyl)-hydr free energy barrier for the glycosylation and transglycosylation oxylamine) (Teze et al., 2013). Our main goal is to gain reactions catalyzed by a transglycosidase of CAZy family 72, understanding and compare the molecular mechanisms was calculated to be of 11 and 16 kcal/mol, respectively (Raich underlying hydrolysis and transglycosylation. It has been et al., 2016). The effect of using different QM/MM partitions proposed that the main driving force for increasing the on the barrier heights of glycosylation and hydrolysis in Osβ- T/H ratio must be the relative destabilization of the gly has also been investigated (Badieyan et al., 2012). It was hydrolysis TS as compared to the transglycosylation one, concluded that a Tyr residue equivalent to Ttβ-gly Tyr284 (thus which implies that differential interactions are involved in interacting with the catalytic nucleophile and the sugar O5 atom), each case. has an important contribution to the energy profile, being the effect more significant for deglycosylation than for glycosylation. METHODOLGY Importantly, the study showed that this residue had to be included in the QM region to obtain reliable potential energy System Setup barriers and ring conformations. The effect of other residues Initial coordinates for the wild-type Ttβ-gly were taken from the from the−1 subsite was also analyzed. corresponding X-ray structure with PDB code 1UG6 (resolution Very few QM/MM studies have compared hydrolysis and 0.99 Å) (Lokanath et al., in press). For the incorporation of the transglycosylation reactions for retaining glycosidases and, to substrates into the enzyme active site, the 1UG6 coordinates the best of our knowledge, none has focused on family were overlaid with those of the Osβ-gly structure crystalized GH1. BB1K:AMBER QM/MM calculations on a β-galactosidase with non-hydrolyzed cellotetraose (CTT) substrate(PDB code: belonging to family GH2 were used to study the catalytic 3F5J, resolution 1.95 Å) (Chuenchor et al., 2011). In this way, mechanism of glycosylation, hydrolysis and transglycosylation we took advantage of the position of CTT and modified it reactions, the latter leading to different regioproducts (Bráa with PyMol program1 , by building the fucose moiety from the et al., 2010). All TSs were characterized as dissociative and CTT glucose ring closer to the catalytic residues (Glu338 and with the proton still attached to the acidic residue, although, Glu164 of the 1UG6 wild-type structure). The second ring of the unfortunately, more detailed geometric information on the CTT substrate was used to build up the (p-nitrophenol (pNP)) transglycosylation TS was not given. The 2-OH group H-bond leaving group of the pNP-Fuc substrate. Two crystallographic with the nucleophile was also related to pyranosyl ring distortion water molecules from the wild-type structure that clashed to the and TS stabilization. The energy barriers for transglycosylation added substrate were removed; the remaining of crystallographic were found to be higher than for hydrolysis by ∼2–4 kcal/mol, waters were kept in the model. Missing hydrogen atoms were and the origin of the observed regioselectivity was found to added to the model and the ionizable residues were protonated be thermodynamic more than kinetic. Recently, the hydrolysis using Propka3.0 (Olsson et al., 2011) at pH = 7.0, except for and transglycosylation reactions catalyzed by a family GH3 β- Glu338 and Glu392 that were modeled as deprotonated, Glu164 glucosidase have been studied by umbrella sampling calculations as protonated and His119 was protonated at Nδ. All molecular at the SCC-DFTB/CHARMM level (Geronimo et al., 2018). Both dynamics (MD) simulations were done using the ff14SB (Maier reactions were found to have similar free energy barriers (∼18 et al., 2015) Amber force field for the protein; GLYCAM06j kcal/mol) but quite different TSs; the predicted TS for hydrolysis (Kirschner et al., 2008) and GAFF (Wang et al., 2004) atom types also differed from the ones previously described. For hydrolysis, and parameters were employed for the sugar moieties and the the TS had the glycosyl-water bond practically formed (1.50 pNP group, respectively. In addition, several missing parameters ± 0.04 Å), resulting in a reduced ionic character of the sugar corresponding to a bond, and several intramolecular angles and which, according to the Cremer-Pople polar coordinates (Cremer dihedrals of the pNP-Fuc substrate were completed searching and Pople, 1975) that describe ring conformations, presented for the corresponding atom types translation in GAFF and a 4 C1 conformation. For transglycosylation, an earlier TS was GLYCAM06j force field databases (those parameters are given predicted (with the new bond only partially formed, 1.9 ± 0.1 in Supplementary Scheme S1 and Supplementary Table S1). Å), with the transferred sugar more positively charged and in a Three sodium ions were added to neutralize the system, which 4 H conformation. In both cases, proton transfer to the catalytic 3 was finally solvated in a cubic box of preequilibrated TIP3P base occurred late, especially in the hydrolysis step, where after approaching the base catalyst at 1.6 Å. 1 The PyMOL Molecular Graphics System, Version 2.0 Schrödinger, LLC. Frontiers in Chemistry | www.frontiersin.org 4 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation water molecules, with no <15 Å from the edge of the water All simulations were carried out using AMBER16 software box to the nearest protein atom. The resulting system contained (GPU (CUDA) version of the PMEMD (Götz et al., 2012; 70,227 atoms. Salomon-Ferrer et al., 2013) package. QM/MM Calculations Molecular Dynamics (MD) Simulation QM/MM calculations were performed with the modular program Details package ChemShell (Sherwood et al., 2003; Metz et al., 2014) The solvated Ttβ-gly/pNP-Fuc model system was relaxed using TURBOMOLE-V6.3 (Ahlrichs et al., 1989) to obtain by performing energy minimization at the MM level with the QM energies and gradients at the DFT level. The PBE0 the steepest descent and conjugate gradient methods. The functional (Adamo and Barone, 1999) was used, as it gave errors minimization protocol consisted in a three-step process, firstly smaller than 0.5 kcal/mol in a recent benchmarking modeling restraining all the system but the ligand (pNP-Fuc), secondly glycosidic bond hydrolysis by glycosidases (Pereira et al., 2017) restraining only the protein atoms (except hydrogen atoms) and a and inclusion of D3 correction did not affect significantly the third step without restraints. Once the system was relaxed, an MD calculated energies. MM energies and gradients were evaluated simulation was performed under periodic boundary conditions by DL_POLY (Smith and Forester, 1996), which was accessed (PBC). The MD protocol consisted in a heating simulation (from through the ChemShell package, using the AMBER force field. 0 to 300 K) during 200 ps under NVT conditions using Langevin The electrostatic embedding scheme (Bakowies and Thiel, 1996) dynamics. Following, an equilibration period of 400 ps was was used within the QM/MM approach to let the MM point performed. The first 200 ps were run under NPT conditions to charges to polarize the electronic density of the QM region. No reach a system density of around 1 g cm−3 and using a isotropic cutoffs were introduced for the nonbonding MM and QM/MM weak-coupling algorithm and the Berendsen barostat (Berendsen interactions. In all the QM/MM calculations, the cubic box of et al., 1984) at 1 atm. Then, the system volume was fixed and the water molecules was simplified to a 30 Å sphere surrounding next 200 ps were run under NVT conditions. An equilibration of the full protein (Figure 2). All residues and water molecules 10 ns and the final production run of 100 ns were performed at within 15 Å from the anomeric center were included in the 300 K using the NVT ensemble without any restraints. Along the optimization process as active region (around 2,100 atoms) while MD, the covalent bonds containing hydrogen were constrained the remaining atoms were kept fixed. using the SHAKE algorithm (Ryckaert et al., 1977), and the Two different QM regions were studied (Figure 3): the particle-mesh Ewald method (Essmann et al., 1995) was used to small QM region containing 84 atoms in the glycosylation treat long-range electrostatic interactions. A 1fs time step was step (the pNP-Fuc substrate, four nearby waters, and the side- used in all the MD trajectories. The analysis of the MD simulation chains of Glu164, Glu338, and Tyr284), and a large QM was carried out using the standard tools of the AMBER 16 region with 138 atoms (including the small QM region and, package (Case et al., 2005) [cpptraj program (Roe and Cheatham, in addition, the side-chains of Asn282, Asn163, Arg75, Gln18, 2013)] and VMD (Humphrey et al., 1996). Based on this analysis, and Glu392). For the hydrolysis and transglycosylation QM/MM we established the criteria for the selection of a few frames of descriptions the O-pNP group was substituted in the QM the MD trajectory that served as starting point for the QM/MM regions by a water molecule and the BnON(Me)-Glc molecule, calculations of the glycosylation (G) reaction. respectively. Three water molecules were added (besides the To obtain the starting structure for the QM/MM calculations nucleophilic water, WAT431) in the QM zones for the hydrolysis of the hydrolysis reaction (H), the product geometry of the calculations whereas only one QM water molecule was left in glycosylation process was modified by deleting all the pNP the transglycosylation study. The number of QM atoms in the coordinates except those of the just formed hydroxyl group small/large regions for the hydrolysis and transglycosylation that was converted into a water molecule using the PyMol steps are 67/121 and 100/154, respectively. The methodology program. The system was solvated with a cubic box of TIP3P based on link atoms (three in the QM(small)/MM models and water molecules, and their coordinates were relaxed at the MM eight in the QM(large)/MM partitions) was used to define the level by constraining all protein atoms and the just built water QM/MM boundary with the charge-shift approach (Claeyssens molecule (acceptor). This was followed by a short MD simulation et al., 2006). The total charge of both QM regions is−1 in the (210 ps) to allow the rearrangement of water molecules in the three catalytic reaction steps. active site. From the last frame of that MD trajectory a QM/MM QM/MM optimizations were carried out employing the minimization was carried out to locate a reactant minimum to limited-memory Broyden-Fletcher-Goldfarb-Shanno (L- initiate the hydrolysis reaction. BFGS) (Liu and Nocedal, 1989) algorithm combined with For the transglycosylation reaction (T), the coordinates of the the Hybrid Delocalized Internal Coordinate Scheme (Billeter above mentioned CTT substrate were used as template to build et al., 2000) as implemented in Chemshell (Metz et al., 2014). those of the BnON(Me)-Glc acceptor substrate, enforcing the 4- Reaction paths were scanned by performing harmonically OH hydroxyl group from the glucose moiety to overlay with the restrained optimizations along a properly defined reaction hydroxyl coordinates of the pNP leaving group from the previous coordinate for each chemical process in steps of 0.2 Å. Reaction glycosylation step. Three water molecules were deleted due to coordinates were defined as linear combinations of the three clashes with the acceptor molecule. The system was relaxed as main distances involved in each step, a definition that we have done for the hydrolysis reactant. successfully used when modeling other carbohydrate-active Frontiers in Chemistry | www.frontiersin.org 5 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation FIGURE 2 | Representation of the enzymatic model used for the QM/MM calculations. The location of the QM region is highlighted. The inset shows as an example the arrangement of the QM atoms for the glycosylation reaction at the active site. corresponding to the maximum of the potential energy profile was taken as indicative of the corresponding transition state. The PBE0 hybrid functional and two basis sets (SVP and TZVP) were used for the description of the small QM region, whereas the large QM region was always treated at the PBE0/TZVP level. In addition, QM(PBE0/TZVP)/MM single point energy calculations on the QM(PBE0/SVP)/MM geometries along the hydrolysis and transglycosylation potential energy profiles were also performed. Natural population analysis (NPA) charges (Reed et al., 1985) were determined from QM/MM calculations with the QM region described at the PBE0/TZVP level. The contribution of the different residues to the QM/MM potential energy barriers of the hydrolysis and transglycosylation steps was examined by setting their point charges to zero in additional single point energy calculations at the QM(PBE0/TZVP)/MM level along the QM(PBE0/TZVP)/MM reaction paths. Umbrella sampling at the QM(SCC-DFTB)/AMBER level FIGURE 3 | Representation of the QM regions used in the QM/MM study of was performed to compute the free energy profile for the the glycosylation reaction. The small QM region includes: the substrate atoms hydrolysis and transglycosylation steps, using the dynamics and water molecules (in black), the catalytic residues (in blue) and the Tyr184 module within ChemShell. The reaction coordinates, defined residue (in red). For the hydrolysis and transglycosylation the O-pNP group as before, were scanned at 0.1 Å intervals using a harmonic was substituted by a water molecule and the BnON(Me)-Glc acceptor substrate, respectively. Three additional waters were included for the force constant of 300 kcal/mol·Å2 . The same force constant hydrolysis and just one for the transglycosylation step. The large QM region was used for restraining the distance between OE1GLU338 comprises for the three reaction steps the corresponding small QM region plus and O2FUC to 2.4 Å in the hydrolysis step to correct for the protein residues depicted in green. undesired proton transfers occurring during the simulations. Still, some of the hydrolysis data had to be discarded in the analyses presented below. Forty-five windows were used enzymes (Gómez et al., 2015). Thus, RC = [d(C1FUC - for hydrolysis and fifty for transglycosylation. Each simulation O4DpNP )-d(C1FUC -OE2GLU338 )-d(HGLU164 -O4DpNP )] for consisted of 20 ps of equilibration and 60 ps of production/data glycosylation, RC = [d(C1FUC -OE2GLU338 )-d(C1FUC -OWAT431 )- collection, under the NVT ensemble and using the Nosé- d(H1WAT431 -OE2GLU164 )] for hydrolysis and RC = [d(C1FUC - Hoover thermostat (Nosé, 1984; Hoover, 1985). All atoms from OE2GLU338 )-d(C1FUC -O4GLC )-d(H4OGLC -OE2GLU164 )] for residues beyond 20 Å of C1FUC were frozen. The histograms transglycosylation. For each reaction step, the structure show sufficient overlap of the windows for the chosen stepsize Frontiers in Chemistry | www.frontiersin.org 6 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation and force constant. The umbrella integration analysis method TABLE 1 | Analysis over the MD simulation of the most populated (Kästner and Thiel, 2005, 2006) was used to compute the free hydrogen-bonds between pNP-Fuc, protein residues, and water molecules at the active site. energy profiles. Acceptor Donor H-bond Average occupancy % Distance (Å) RESULTS AND DISCUSSION Glycosylation Step ONO2pNP Solvent 99.6 2.79 OE1GLU164 ND2ASN282 92.5 2.81 Molecular Dynamics Simulations O4DpNP OE2GLU164 90.6 2.76 To obtain a first structural insight into the dynamics of OE2GLU338 OHTYR284 89.0 2.76 the Ttβ-gly/pNP-Fuc Michaelis complex, the root-mean-square deviations (RMSDs) of the protein backbone and of the heavy OE1GLN18 O3FUC 79.3 2.72 atoms of the pNP-Fuc substrate were calculated along the 100 OE1GLU392 O4FUC 62.5 2.80 ns of production trajectory (Supplementary Figure S1). The O2FUC ND2ASN163 46.8 2.90 protein structure is well equilibrated. The substrate RMSD (with OE1GLU338 O2FUC 46.6 2.84 an average value of 0.78 ± 0.37 Å) presents an oscillatory O3FUC NE2HIE119 34.4 2.87 pattern because the pNP phenyl group is flipping between two O5FUC Solvent 24.0 2.85 conformations all along the trajectory. As it will be detailed O4FUC Solvent 20.4 2.84 below, the Fuc ring is more rigidly bound at the−1 subsite than O3FUC NE1TRP393 16.6 2.91 the pNP group is at the larger +1 subsite, thus having this latter OE2GLU338 O2FUC 16.3 2.85 moiety of the substrate more room to move along the simulation. OE1GLU338 NH1ARG75 13.6 2.91 To analyze in more detail the molecular interactions more H-bond occupancies are defined as the fraction of frames the bond is present. Only H- likely to play a role in the catalytic mechanism, a search of bond occupancies above 10% are included. The average donor to acceptor heavy atoms distances of the bonds when present are also given. the most populated hydrogen bonds between the substrate, protein residues and water molecules at the active site was carried out along the MD simulation. In Table 1 a summary of these hydrogen bonds is presented (a more complete analysis each Fuc···water interaction), and the axial 4-OH of Fuc is H- is given in Supplementary Tables S2–S5). The fraction of bonded to OE1GLU392 (62.5% occupancy). This latter interaction frames the hydrogen bond is present (the H-bond occupancy), structurally differs from that reported by Badieyan et al. (2012) along with the corresponding donor···acceptor average distance, when cellobiose is hydrolyzed by Osβ-gly, as in Glc this hydroxyl are given. An occupancy higher than 90% is obtained for group is in the equatorial position, which makes it interact with the interaction between the carboxylate oxygen OE1 of the the other oxygen of the Glu392 carboxylate. catalytic acid/base residue Glu164 and HD22 of Asn282. This To summarize, the glycosyl moiety of the substrate bound interaction helps to maintain (with an occupancy of 90.6%) the at the−1 subsite establishes a complex H-bond network with hydrogen bond between the proton of Glu164 and the glycosidic the protein and some water molecules that contributes to its oxygen (O4DpNP ). This is a very important interaction for the stabilization. It is worth highlighting that all the first- and second- glycosylation process because in this way Glu164 remains well- shell protein residues around the−1 subsite identified in this positioned (average heavy atoms distance of 2.76 Å) for proton H-bond analysis are highly conserved residues among the GH1 donation to the leaving group. The catalytic nucleophilic residue family (Bráa et al., 2010; Camargo et al., 2019). In contrast, the Glu338 mainly interacts through its carboxylate oxygen OE2 H-bond analysis did not identify any conserved residue around with Tyr284 (88.9% occupancy) so that the negative charge the +1 site so confirming the lack of specificity for the aglycone on the nucleophile becomes stabilized. In addition, this same moiety within the GH1 family as previously observed (Bráa carboxylate group of the nucleophilic residue Glu338 (through et al., 2010). The pNP molecule only establishes H-bonds with its OE1 oxygen) establishes nearly half of the time another two water molecules all along the MD trajectory through its hydrogen bond with the glycosyl 2-OH group of the substrate. NO2 group. In addition to those H-bond interactions, the Fuc Some oscillation of this 2-OH group between OE1GLU338 ring remains surrounded all along the trajectory by hydrophobic (46.6% occupancy) and OE2Glu338 (16.3% occupancy) is actually interactions established with conserved residues such as Trp385, observed. As commented in the Introduction, this interaction Trp393, Trp120, and also Phe401, whereas the pNP leaving group is key for the stabilization of the glycosylation transition state, basically interacts with Trp312. especially for retaining β-glycosidases because the interaction To initiate the QM/MM study of the glycosylation catalytic with the nucleophile will affect oxocarbenium cation formation. mechanism, two MD snapshots of the Ttβ-gly/pNP-Fuc Moreover, OE1Glu338 is also H-bonded with Arg75 side chain and Michaelis complex were selected. Those structures were filtered a water molecule. On the other hand, the 2-OH group of the Fuc by imposing the presence of the direct H-bonds (O4DpNP - ring also interacts with Asn163 side-chain (with an occupancy HGLU164 and OE2Glu338 -H2OFUC ) between the two catalytic of 46.8%). The 3-OH group of Fuc establishes a H-bond with residues, Glu338 and Glu164, with the pNP-Fuc substrate. In HE2HIS119 (34.4% occupancy) and with OE1GLN18 (79.3%). The addition, the presence of the H-bond between Tyr284 and 4-OH and 5-OH of Fuc are hydrogen bonded as acceptors to Glu338 was also verified in the selected frames. The relevance the same water molecule (with an occupancy of around 20% for of this interaction (and of describing it at the QM level) for the Frontiers in Chemistry | www.frontiersin.org 7 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation energetics of the glycosylation step and hydrolysis steps has been TABLE 2 | QM/MM potential energy barriers (1V6= ) and reaction energies (1VR ) recently highlighted by Badieyan et al. (2012) in their QM/MM calculated for the glycosylation step for the two complexes studied (GI and GII ). study of Osβ-gly mentioned above. This is the reason why we GI GII have included Tyr284 in the QM(small)/MM partition. Finally, in the two frames selected it was also verified that the additional QM level 1V6= 1VR 1V6= 1VR predominant H-bond interactions highlighted in Table 1 were PBE0/SVPa 28.4 7.1 30.2 0.4 also present. PBE0/TZVP//PBE0/SVPa 28.8 5.1 32.1 −0.1 PBE0/TZVPa 27.9 6.3 30.2 0.6 QM/MM Calculations PBE0/TZVPb 22.0 4.8 20.4 6.6 In this section the results of the QM/MM calculations on the glycosylation reaction will be presented. This catalytic step a QM(small)/MM. consists in the cleavage of the glycosidic bond of the pNP-Fuc b QM(large)/MM. molecule to form a covalent glycosyl-enzyme intermediate with Energies are given in kcal/mol. Glu338. Despite pNP is a good leaving group, acid catalysis is required for this mechanistic step, a role that is carried out by barriers reflect that they are very sensitive to different QM/MM the protonated Glu164 (Figure 1). partitioning schemes and that non-catalytic residues (other than The two selected frames from the MD trajectory (denoted Tyr284) might play a role in the stabilization of the glycosylation as GI and GII in this paper) were optimized using the two transition state, as previously stated. As for the reaction energies, QM/MM partitions described above and the corresponding the glycosylation process turns out to be more endoergic for minima were located and characterized. From those minimum frame GI than for frame GII when the QM(small)/MM partition energy structures, the potential energy profiles along the is used. However, both glycosylation reactions become clearly glycosylation step were calculated as a function of the reaction endoergic when the QM region is enlarged. Interestingly, the coordinate RC = [d(C1FUC -O4DpNP )-d(C1FUC -OE2GLU338 )- endothermic nature of the glycosylation reaction was recently d(HGLU164 -O4DpNP )]. In Table 2 the QM/MM potential energy observed in previous studies where their authors suggest that an barriers (obtained from the corresponding energy profiles) endothermic glycosylation might be a prerequisite for efficient and the reaction energies (considering the optimized glycosyl- transglycosylation (Raich et al., 2016). enzyme product) for the glycosylation step at the QM = In Figure 4 the evolution of the four main interatomic PBE0/SVP and PBE0/TZVP levels using the QM(small)/MM distances involved in the glycosylation reaction (C1FUC - and QM(large)/MM partitions are given for both frames. The O4DpNP , C1FUC - OE2GLU338 , HGLU164 - O4DpNP , and HGLU164 corresponding energy values obtained by means of PBE0/TZVP - OE2GLU164 ) are plotted for the GI frame along the reaction single-point energy calculations on the PBE0/SVP(small)/MM coordinate calculated with the QM(large)/MM partition and geometries are also given. It can be observed that the potential at the PBE0/TZVP level for the QM region. In Figure 5 energy barrier using the QM(small)/MM partition takes values the molecular representation of the corresponding reactant, in between 28 and 30 kcal/mol when the smaller basis set is transition state and product are depicted. It can be observed used. The energy correction using the larger basis set tends to that the C1FUC - OE2GLU338 distance between the nucleophile increase those energy barriers (by 0.4–2 kcal/mol) especially for residue and the anomeric carbon gradually diminishes from 3.33 the GII frame. However, when the geometry optimization is Å at reactants (RC= −3.90 Å) to a value of 2.61 Å at the carried out with the TZVP basis set, no significant differences transition state (RC= −2.10 Å) (what corresponds to a bond not are obtained with respect to the SVP energy values. These energy yet formed), and finally takes a value of 1.50 Å at the covalent barriers are in the range of that calculated with a large QM glycosyl-enzyme intermediate (RC= 0.70 Å) in which this bond region for the glycosylation reaction involving cellobiose and is already completely formed. The length of the glycosidic bond Osβ-gly (Badieyan et al., 2012), and also of those calculated for C1FUC - O4DpNP is 1.41 Å at the reactants, increases to 2.33 Å at other retaining GHs (Petersen et al., 2009, 2010; Biarnés et al., the TS (what corresponds to the glycosylic bond already broken) 2011). However, a large effect is observed in the present system and finally takes a value of 3.17 Å at the product structure in when the large QM region is used with the TZVP basis set. which the leaving group has been released. For the second frame Lower potential energy barriers (by 5.9 and 9.8 kcal/mol for (GII ) the distances are very similar (Supplementary Table S6). frames GI and GII , respectively) were obtained for both frames. These structural results are coherent with the dissociative nature Notice that, as pNP is a better leaving group than the Glc of of the glycosylation transition state. The oxocarbenium character cellobiose modeled in the glycosylation reaction catalyzed by of this TS is also made apparent by a higher degree of double Osβ-gly (Badieyan et al., 2012), a lower energy barrier should be bond character of the C1FUC -O5FUC bond (the C1FUC -O5FUC expected. In particular, the potential energy barriers calculated at distance diminishes 0.14 Å from the reactant to the TS structure QM(large)/MM level and the TZVP basis set are 22.0 kcal/mol for both frames). In addition, the acidic residue Glu164 gradually (for GI ) and 20.4 kcal/mol (GII ). These values are in better approaches the glycosidic bond along the two reaction paths and in reasonable qualitative agreement with the value of 17.1 (the O4DpNP - HGLU164 hydrogen bond diminishes from 1.99 kcal/mol for the phenomenological free energy of activation Å (GI ) and 2.02 Å (GII ) at the reactants to 1.82 and 1.88 derived from the experimental kcat at 40◦ C (Teze, 2012). The Å at the GI and GII TS, respectively) although the proton significant differences between the calculated potential energy is not yet transferred. Proton transfer takes place when the Frontiers in Chemistry | www.frontiersin.org 8 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation FIGURE 4 | Evolution of C1FUC -O4DpNP (in blue), C1FUC -OE2GLU338 (in orange), HGLU164 -O4DpNP (in purple), and HGLU164 -OE2GLU164 (in green) distances (in Å) along the glycosylation reaction. The results correspond to frame GI calculated with the QM(large)/MM partition and at the QM(PBE0/TZVP) level for the QM region. OE2GLU338 -C1FUC distance is ∼2 Å, and is followed by complete formation of the C1FUC - OE2GLU338 bond, the Tyr284 hydroxyl OE2GLU338 -C1FUC bond formation. This is consistent with a group separates from Glu338 and clearly approaches (with a change in pKa of Glu164 that makes proton transfer easier as reduction in the bond distance of around 1.3 Å at the TSs and the negative charge of the leaving group gradually increases and of around 1.9 Å at the product structures) the O5FUC atom. The that of the nucleophile gradually decreases with the progress of expected oxocarbenium ion-like character of the TS is favored by the glycosylation reaction. When a smaller QM region is used the coplanarity of the C2, C1, O5, and C5 atoms of the fucose within the QM/MM approach, and independently of the basis ring which adopts a slightly distorted chair conformation close to set employed for the QM description, the transition state can a half-chair (4 H3 ) and with some 4 E character. also be described as dissociative (Supplementary Table S7 and At this point it is interesting to comment that it has Supplementary Figure S2). Curiously, the O4DpNP - HGLU164 H- been shown that GHs usually bind the substrate in a subtle bond distance at the reactant varies by 0.42 Å between the two but critical distorted conformation (Davies et al., 2012). Such frames considered (GI and GII ), and a significant difference in conformational change brings the glycosidic bond closer to the this distance is also observed at the corresponding TSs. When equatorial orientation and to a more planar ring structure (so applying the QM(large)/MM partition, these differences between closer to the oxocarbenium TS geometry) that facilitates in-line frames disappear. With the smaller QM region, the nucleophilic nucleophilic attack. A possible exception to this has been recently attack by Glu338 is completed at lower RC values than with the reported for family GH3 enzymes, for which the glucose bound QM(large)/MM partition and quite before the proton is fully at the−1 subsite seems to adopt a 4 C1 conformation (maybe transferred. Thus, the description at the QM level of the H- slightly distorted to 4 H5 ) in the Michaelis complex according to bonds that Glu338 and Glu164 establish with their neighboring QM/MM calculations and crystallographic data (Hrmova et al., residues seems to affect significantly their nucleophilic and 2001; Geronimo et al., 2018). Different conformational itineraries acid/base characteristics. (from reactants to TS and to products) have been associated with The inclusion of other protein residues besides the catalytic different GHs (Ardèvol and Rovira, 2015). Thus, for retaining β- ones (Glu338, and Glu164) and Tyr184 in the zone treated glucosidases acting on glucose substrates the 1 S3 → 4 H3 → quantum-mechanically (as in our QM(large)/MM model) also 4 C itinerary has been described for the glycosylation step. The 1 affects the stabilization by the enzyme environment of the QM/MM study by Badieyan et al. (2012) of the glycosylation oxocarbenium ion and the negatively charged glycosyl oxygen step catalyzed by rice Osβ-gly using cellobiose as substrate anion. Moreover, as mentioned in the Introduction, one agrees with this itinerary. However, in the present study, in of the most stabilizing protein-substrate interactions at the which we have pNP-Fuc instead of Glc(β1-4)β-Glc as substrate, glycosylation TS is the H-bond between the 2-OH group of we were unable to characterize a Michaelis complex in a 1 S3 the sugar moiety and the nucleophilic residue Glu338. This conformation. The QM(large)/MM optimized reactant has a stabilizing contribution correlates with a reduction of the slightly distorted 4 C1 conformation (with θ and φ angles values OE1GLU338 - H2OFUC H-bond distance, from 2.09 Å and 2.37 Å of 216 and 27◦ , respectively). Attempts to obtain a different at reactants to 1.82 and 1.85 Å at the TS of GI and GII , respectively conformation by maintaining the ring restrained previous to the (Supplementary Table S6). The Fuc at the −1 subsite is also full QM(large)/MM minimization only produced an alternative stabilized along the glycosylation pathway by the HD21ASN163 minimum close to a 4 E (slightly 4 H3 ) conformation. This - O2FUC interaction as the corresponding H-bond distance structure presents a small increase in the C1FUC – O4DpNP diminishes around 0.10 Å from reactants to the transition state distance (1.49 Å) and slightly shorter C1FUC – OE2GLU338 structures. Tyr284 plays a relevant role along the glycosylation (3.29 Å), C1FUC – O5FUC (1.36 Å), and HGLU164 – O4DpNP process. As the OE2GLU338 atom loses negative charge upon (1.77 Å) distances than the previous one. The potential energy Frontiers in Chemistry | www.frontiersin.org 9 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation favorable conformation that the enzyme will bind and, thus, it has no catalytic significance. Badieyan et al. (2012) reported the importance of having Glu440 of Osβ-gly (that interacts with the 4-OH and 6-OH groups of Glc) in the QM region in order to characterize the 1 S3 conformation. In the present study, the QM(large)/MM partition includes the equivalent Glu392 residue, amongst many others. At the origin of these differences may be the fact that we are modeling Fuc instead of Glc (they differ in the position of the 4-OH group and Fuc has no 6- OH group but a methyl in this position); or that our substrate is not an oligosaccharide but a synthetic phenyl derivative that may introduce less steric requirements than a monosaccharide at the +1 subsite. Nevertheless, as glycosidic bond breakage occurs early in the reaction, the ring is distorted early along the calculated reaction coordinate with subsequent attack by the nucleophilic Glu338 occurring in the expected orientation. Hydrolysis and Transglycosylation Steps (Deglycosylation) Once the pNP molecule leaves the enzyme active site, a water molecule can compete as nucleophile in a hydrolysis reaction with the acceptor ligand (BnON(Me)-Glc), that also attacks as a nucleophile on the anomeric carbon of the fucose molecule in a transglycosylation reaction leading to the N-methyl-O-benzyl-N- (β-D-fucopyranosyl(1-4)β-D-glucopyranosyl)-hydroxylamine product. The attacking water molecule as well as the BnON(Me)- Glc ligand are activated by a proton transfer to the Glu164 residue that acts as a base along those two nucleophilic processes. Starting from the minimized MD structure (see Methodology section) in which a water molecule (WAT431) was located at the active site, close enough to the anomeric carbon of fucose ring, we calculated the QM/MM energy profile for the hydrolysis process using a reaction coordinate defined as RC = [d(C1FUC -OE2GLU338 )-d(C1FUC -OWAT431 )-d(H1WAT431 - OE2GLU164 )]. As mentioned above, a H-bond network of three water molecules, in addition to the nucleophilic WAT431, were included in the QM region for the hydrolysis simulation. On the other hand, for the transglycosylation step, initial energy profiles starting from minimized MD structures gave highly overestimated energy barriers (>34 kcal/mol), suggesting some structural deficiencies in the model. With the aim of improving this model, we followed a strategy successfully used by our group in previous studies of enzyme catalyzed reactions involving sugar-transfer between donor and acceptor substrates (Albesa- Jové et al., 2015; Mendoza et al., 2016, 2017). The refinement protocol consists on running a QM(SCC-DFTB)/MM MD of, FIGURE 5 | Representative structures for reactant (A), transition state (B), and in this case, the transglycosylation product, followed by scan product (C) for the glycosylation reaction. The depicted geometries calculations of frames from this MD to generate new models correspond to frame GI within the QM(large)/MM partition and with the QM for the transglycosylation reactant that is then minimized at the region described at the PBE0/TZVP level. Distances are in Å. QM(PBE0)/MM level. By doing this, interactions between donor and acceptor (or also involving water molecules) that are difficult to capture in the MM MD seem to be better represented. In the profile calculated from this minimum gives a barrier height of present system, the structural changes observed upon refinement only 10 kcal/mol (Supplementary Figure S3), which is much involve the reorientation of the 4-OH hydroxyl of the acceptor smaller than the experimental value. This suggests that, although Glc toward the Fuc ring. This displaces a water molecule that was such conformation may be possible, it may not be the most located between the two and that most likely hindered acceptor Frontiers in Chemistry | www.frontiersin.org 10 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation approach during reaction (Supplementary Figure S4). In the presented and the corresponding molecular structures new reactant structure, the Glc hydroxyl will be able to better of the different reactants, transition states and products interact with Fuc as transglycosylation takes place. Substrate- are depicted in Figure 6 (for QM(small)/MM results see substrate interactions have been shown to have a very important Supplementary Table S8). At reactants, with QM(large)/MM, contribution in the sugar transfer catalyzed by retaining the initial location of the water molecule in the hydrolysis glycosyltransferases (Gómez et al., 2012, 2013; Mendoza et al., process corresponds to a slightly more distant position from 2017) and could, thus, be also relevant for transglycosylation. the anomeric carbon (3.45 Å) than in the two reactants of the Two frames (called here TI and TII ) were generated by transglycosylation process (3.27 and 3.35 Å). The hydrolysis this approach. From the two optimized QM/MM reactant values reported for Glc hydrolysis by rice β-glucosidase (Osβ-gly) structures, the transglycosylation potential energy profiles for the are slightly larger (3.62–3.75 Å) (Badieyan et al., 2012; Wang corresponding β(1-4) glycosidic linkages were calculated using as et al., 2013), as occurs at the same level of theory when the reaction coordinate the following expression RC = [d(C1FUC - QM(small)/MM partition is used (Supplementary Table S8). OE2GLU338 )-d(C1FUC -O4GLC )-d(H4OGLC -OE2GLU164 )]. In this Also at TS structures the C1FUC -OAcc bond, that is still forming, model only one water molecule remains in the QM region. We is slightly longer for hydrolysis (with distances of 2.27, 2.11, did not analyze other glycosidic linkages between the Fuc and and 2.12 Å at H, TI and TII transition states, respectively), BnON(Me)-Glc because it has been proven experimentally that while the nucleophile Glu338 side chain is completely displaced the native enzyme does not produce other product regioisomers from the sugar in both reactions but slightly closer from C1FUC when these two substrates are used (Teze et al., 2013). in the case of hydrolysis (with C1FUC -OE2GLU338 distances In Table 3 the potential energy barriers and the of 3.12, 3.42, and 3.34 Å at H, TI , and TII transition states, reaction energies corresponding to the hydrolysis and the respectively). For Glc hydrolysis by Osβ-gly, slightly shorter transglycosylation steps are given for the two QM/MM partitions distances have been reported, especially for the C1GLC -ONu one used in this study and at the different QM levels of electronic (1.89–2.14 and 2.45–2.51 Å, for C1GLC -OAcc and C1GLC -ONu , structure theory analyzed. The first general observation is respectively) (Badieyan et al., 2012; Wang et al., 2013). Along the that the potential energy barrier for the hydrolysis step is nucleophilic attacks, the proton from the water molecule or the very similar to that of the transglycosylation one, as expected BnON(Me)-Glc acceptor comes closer to the general base residue because the experimental result for this reaction catalyzed by Glu164 but it has not been transferred yet at the transition states WT Ttβ-gly gives a transglycosylation yield of 36%, which (HAcc -OE2GLU164 distances between 1.42 and 1.48 Å). This is implies differences in the energy barriers no larger than 1 in agreement with the hydrolysis description given on previous kcal/mol. With the small QM region at the PBE0/SVP level, the works (Badieyan et al., 2012; Wang et al., 2013; Geronimo et al., hydrolysis potential energy barrier is 2.1 and 0.9 kcal/mol higher 2018). Interestingly, this proton transfer is more advanced in the than for the transglycosylation frames TI and TII , respectively. hydrolysis and transglycosylation TSs that it was the Glu164 to According to this result, the transglycosylation reaction would pNP one in the glycosylation TS (d(HGLU164 -O4DpNP ∼1.8Å). be faster than the hydrolysis process. This trend is inverted when Notice that this is the first time that the transglycosylation step single-point energy calculations are carried out at a higher level is studied for wild-type family GH1, thus direct comparison (QM(PBE0/TZVP)), being the hydrolysis energy barrier now for this step with previous works is not possible. For family slightly lower (by 0.7 (TI ) and 0.9 (TII ) kcal/mol). The correct GH2, potential energy scans followed by TS characterization trend between the hydrolysis and transglycosylation barriers identified a dissociative TS for transglycosylation and with is confirmed when a bigger basis set is used in the QM/MM the proton not yet transferred; however, detailed structural optimizations and also when the QM region is enlarged. As information (i.e., distances) was not provided (Bráa et al., 2010). for the glycosylation step, the values of the energy barriers Besides, for the GH2 enzyme a magnesium ion interacting calculated with the smaller QM region are around 6 kcal/mol with the acid/base residue participates in catalysis, which is higher than those obtained with the larger QM region. Those likely to introduce some structural differences when compared energy differences reflect again that the two reaction processes to the present system (e.g., the hydrolysis TS presented the are very sensitive to different QM/MM partitioning schemes. anomeric C1 atom equidistant from the nucleophile and the Despite this, the difference between transglycosylation and attacking water (2.25 Å), which differs from all the results hydrolysis energy barriers is similar (1.2 kcal/mol for TI and 2.6 commented above). kcal/mol for TII ). Thus, the energy trends might be reproduced Despite the sensitivity observed for these TS geometries to at a lower computational cost (QM(PBE0/TZVP)/MM but the DFT functional, the basis set, the QM/MM partition, the QM(small)/MM). As for the reaction energies, the hydrolysis characterization method, the substrates or the specific enzyme, process is more exoergic than the transglycosylation reaction. it is clear that both the hydrolysis and the transglycosylation This result is in agreement with previous results by Bráa transition states are very loose structurally, and they have a et al. (2010) in their study of the catalytic mechanism of a dissociative nature. From the distances reported in Table 4, β-galactosidase and with the known thermodynamically favored which are thus to be taken as orientative, some trends can be hydrolysis over synthesis reaction (Planas and Faijes, 2002). extracted at the QM(PBE0/TZVP)/MM with QM(large)/MM In Table 4 the most relevant interatomic distances (calculated level. For example, the sum of the C1FUC -OAcc and C1FUC - with the QM(large)/MM at the QM(PBE0/TZVP)/MM level) OE2GLU338 distances for the deglycosylation TSs (5.39, 5.53, along the hydrolysis and transglycosylation pathways are and 5.46 Å for H, TI , and TII transition states, respectively) Frontiers in Chemistry | www.frontiersin.org 11 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation TABLE 3 | Potential energy barriers (1V6= ) for the hydrolysis (H) and transglycosylation processes (TI and TII , corresponding to the two different frames studied) along with the corresponding reaction energies (1VR ) with the two QM/MM partitions and at the different QM levels employed. PBE0/SVPa PBE0/TZVP//PBE0/SVPa PBE0/TZVPa PBE0/TZVPb 1V6= 1VR 1V6= 1VR 1V6= 1VR 1V6= 1VR H 25.4 −16.7 22.7 −17.8 24.6 −15.8 18.3 −17.4 TI 23.3 −9.7 23.4 −8.5 25.6 −8.3 19.5 −12.5 TII 24.5 −8.9 23.6 −4.6 27.9 −5.8 20.9 −4.0 a QM(small)/MM. b QM(large)/MM. All energies are given in kcal/mol. TABLE 4 | Distances (in Å) between selected atoms involved in each reaction step for the reactant, transition state (TS) and product of the hydrolysis (H) and transglycosylation (TI and TII , corresponding to the two different frames studied) steps. Reactant TS Product H TI TII H TI TII H TI TII d(C1FUC -OE2GLU338 ) 1.49 1.51 1.49 3.12 3.42 3.34 3.32 3.37 3.31 d(C1FUC -OAcc ) 3.45 3.27 3.35 2.27 2.11 2.12 1.40 1.42 1.43 d(C1FUC -O5FUC ) 1.36 1.35 1.36 1.24 1.25 1.25 1.39 1.38 1.38 d(HAcc -OE2GLU164 ) 1.81 1.80 1.74 1.48 1.48 1.42 0.96 0.96 0.96 d(HTYR284 -O5FUC ) 2.08 2.19 2.13 3.00 3.22 3.19 3.94 3.98 3.81 d(HTYR284 -OE2GLU338 ) 2.46 2.46 2.48 1.90 1.86 1.87 1.82 1.83 1.82 d/HARG75 -OE1GLU338 ) 2.19 2.25 2.25 2.09 2.27 2.10 2.05 2.11 2.08 d(H2OFUC -OE1GLU338 ) 1.78 1.78 1.75 1.89 1.84 1.75 2.05 2.15 1.92 d(H4OFUC -OGLU392 ) 2.29 2.73 4.13 2.16 2.56 4.10 1.93 2.00 4.00 d(HASN163 -O2FUC ) 1.97 1.94 1.93 1.96 1.94 1.93 1.96 2.06 1.98 d(H3OGLC -O4FUC ) – 2.86 2.83 – 2.74 2.45 – 3.11 2.35 d(H3OGLC -O5FUC ) – 3.96 4.06 – 3.43 3.55 – 1.88 2.35 d(H3OGLC -O4GLC ) – 2.34 2.33 – 2.41 2.41 – 2.53 2.57 d(O3GLC -H2WAT433 ) – 2.00 2.23 – 2.03 2.25 – 2.01 2.25 d(OWAT431 -H2WAT432 ) 1.82 – – 1.86 – – 2.77 – – d(H2WAT431 -OWAT433 ) 2.02 – – 2.19 – – 2.33 – – d(H2WAT432 -O5FUC ) 3.58 – – 3.33 – – 1.95 – – d(H1WAT432 -O4 FUC ) 2.28 – – 2.25 – – 2.59 – – The results correspond to QM(large)/MM and at the PBE0/TZVP level. The Acc subscript refers to the acceptor water and glucose moieties in hydrolysis and transglycosylation processes, respectively. are all significantly larger than the sum of the analogous (RC values of −0.63, −0.17, and −0.20 Å for H, TI , and TII TS C1FUC -O4DpNP and C1FUC -OE2GLU338 distances at the structures, respectively). glycosylation TS (which is 4.97 Å on average). On the other In an attempt to determine the robustness of the geometrical hand, the differences between C1FUC -OAcc and C1FUC - differences found between the hydrolysis and transglycosylation OE2GLU338 distances in the deglycosylation TSs (0.85, 1.31, TSs, the corresponding free energy profiles at the QM(SCC- and 1.22 Å for H, TI , and TII TS structures, respectively) are DFTB)/MM level were obtained by umbrella sampling molecular much larger than the corresponding glycosylation value (0.27 dynamics simulations (Supplementary Figures S5, S6). Both Å). Notice also that the hydrolysis value is significantly smaller reactions show very similar free energy barriers (21.10 and 21.93 than the transglycosylation ones, and of the same order as kcal/mol for H and T, respectively), being the transglycosylation the values obtained for Glc hydrolysis by Osβ-gly (0.31–0.62 one slightly higher. At the corresponding TSs, the average C1FUC - Å) (Badieyan et al., 2012; Wang et al., 2013). These results OAcc , C1FUC -OE2Glu338 and HAcc -OE2Glu164 distances are 1.95 ± are in agreement with experimental observations indicating 0.15, 2.60 ± 0.18, and 1.52 ± 0.26 Å for hydrolysis, and 2.11 ± a more SN 2 character for the glycosylation step than for the 0.05, 2.31 ± 0.08, and 1.32 ± 0.03 Å for transglycosylation. Thus, deglycosylation one (Kempton and Withers, 1992). Moreover, the dissociative character of the TSs is maintained, as well as the transglycosylation is predicted to have a transition state more late proton transfer, although at this level of theory the hydrolysis advanced on the reaction coordinate than the hydrolysis one TS would be more advanced on the reaction coordinate (RCH Frontiers in Chemistry | www.frontiersin.org 12 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation FIGURE 6 | Representative structures for reactant (R), transition state (TS), and product (P) for hydrolysis (A) and transglycosylation reactions (frame TI (B) and frame TII (C). The depicted geometries were calculated within the QM(large)/MM partition and with the QM region described at the PBE0/TZVP level. Distances are in Å. = −0.87 Å) than the transglycosylation one (RCT = −1.04 Å). ones (sum of C1FUC -OAcc and C1FUC -OE2GLU338 distances) is However, it is important to note the significant fluctuations on already present at the potential energy maxima, indicating that the distances measured for the hydrolysis TS, which may be due it is likely a consequence of the choice of QM level and not to the flatness of the free energy profile at this zone and are to the umbrella sampling simulations. QM(SCC-DFTB)/MM indicative of a wide ensemble of available TS configurations with umbrella sampling simulations have been reported for hydrolysis different reaction coordinate values. For the transglycosylation and transglycosylation of Glc substrates catalyzed by a family reaction, these fluctuations are much smaller. Interestingly, GH3 β-glucosidase. Transglycosylation TS presented C1GLC - the transglycosylation QM(SCC-DFTB)/MM potential energy OAcc , C1GLC -ONuc and HAcc -OBase distances of 3.0 ± 0.2, 1.9 maximum is more advanced than the hydrolysis one (as with ± 0.1, and 1.6 ± 0.1 Å, which are relatively similar to those QM=PBE0/TZVP), which suggests that the apparent switch in obtained here at the QM=PBE0/TZVP level. However, the TS localization along the reaction coordinate is a consequence hydrolysis TS had the C1GLC -OAcc bond practically formed of the introduction of dynamics. On the other hand, the more (1.5 ± 0.04 Å) and the proton at a distance of 2.06 ± 0.08 compact TS geometries when compared to the QM=PBE0/TZVP Å from the catalytic base (although right before the TS this Frontiers in Chemistry | www.frontiersin.org 13 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation distance had shortened to 1.4–1.5 Å) (Geronimo et al., 2018). et al., 2015). It is interesting to notice that the hydrogen atom Therefore, water was predicted not to be activated by the basic of this 4-OHFUC group is in turn interacting with Glu392; residue, which is not what we observe in our calculations. The the H4OFUC -OE1GLU392 interaction becomes shorter along the authors also related these geometric parameters to the lack of hydrolysis and the TI paths (it is not present along the TII oxocarbenium ion-like character of this TS and a preference energy profile because a water molecule is located as a H- of the GH3 −1 subsite for a TS in a 4 C1 conformation. bonding bridge between O4FUC and OE1GLU392 ). Moreover, (Geronimo et al., 2018). In the present system, the C1FUC -O5FUC the presence of Fuc instead of Glc as substrate (differing on bond becomes shorter from reactants to the corresponding the 4-hydroxyl configuration) has provoked a rotation of the transition states accordingly to its partial double bond character Glu392 sidechain. The analysis of electrostatic contributions along the nucleophilic attack. Hence, the positive charge on to the potential energy barrier of residues at the −1 subsite C1FUC increases in between 0.24 and 0.27 au from reactants (Figure 7) show that Glu392 is stabilizing both the hydrolysis to the transition state structures (Supplementary Table S9 and transglycosylation TSs, being the effect larger for hydrolysis. for QM(PBE0/TZVP)/MM and QM(large)/MM) while O5FUC In fact, this residue is the one showing the largest difference in the fucose ring becomes less negative (around 0.06 au) between the two competitive reactions, suggesting that altering (with the QM(small)/MM partition the change on C1FUC is this interaction could be a way to increase the T/H ratio. This slightly bigger (0.29–0.30 au), Supplementary Table S10). Those observation could be related to the reported transglycosylation charge variations confirm the oxocarbenium ion-like character increase in the Asn390Ile and Phe401Ser mutants, both residues observed in deglycosylation reactions (hydrolysis as well as in contact with Glu392; Asn390 is H-bonding to Glu392 so transglycosylation) in the GH1 β-glucosidases family. The ring that its mutation is likely to alter the Glu392-Fuc interaction. puckering is very similar in both reactions, with Cremer-Pople φ Three other residues interact with Fuc hydroxyl groups: values of 41◦ (for H), 44 and 47◦ (for TI and TII ) and φ angles In this case, an according to the electrostatic interactions of 224◦ (H), 224, and 221◦ (for TI and TII ), which correspond to analysis, they provide a small destabilization effect in this a conformation between 4 H3 and 4 E. For the transglycosylation deglycosylation steps: Asn163, which interacts with the 2-OHFUC TS in family GH3 β-glucosidase mentioned above, a 4 H3 group, with similar values for hydrolysis and transglycosylation; conformation was predicted. (Geronimo et al., 2018). His119, that interacts with O3FUC and O2FUC , and provides It can be observed from the QM(PBE0/TZVP)/MM data a slightly larger destabilization effect for transglycosylation; in Supplementary Table S9 that there are only very small and Gln18, that interacts with O3FUC with similar effect for differences on the atomic charges of the atoms involved in the both TSs. hydrolysis and transglycosylation reactions, in accordance with In the hydrolysis and transglycosylation steps here studied, the small differences obtained in the energy barriers for the the nucleophile Glu338 presents strong interactions with the two kinds of deglycosylation processes. The lower energy barrier 2-OHFUC group, Tyr284, Arg75, and Asn282. The H2OFUC - for hydrolysis might be due then to the more negative charge OE1GLU338 distance increases as the Glu338 residue moves away on the oxygen atom of the nucleophilic water in comparison with the breakage of the C1FUC -OE2GLU338 , although the H- to the negative charge value on the O4 atom of the glucose bond interaction is maintained along the three reaction pathways acceptor molecule. and the negative charge of Glu338 also increases (around 0.3 au As can be seen in Figure 6 and Table 4, some interactions from reactants to the transition states, Supplementary Table S9). between the substrates or with the surrounding network of Concomitantly with those molecular changes, Tyr384 side-chain water molecules also vary along the two reactions studied. For moves away from O5FUC while it approaches OE2GLU338 atom transglycosylation, the hydrogen of the 3-OH group of the (which, as mentioned, is increasing its negative charge). It has Glc acceptor substrate (adjacent to the attacking 4-OHGLC ), in been suggested that the short hydrogen bond with the hydroxyl going from reactant to TS, approaches O5FUC (by 0.53 and group of the equivalent tyrosine residue of a β-galactosidase 0.38 Å for TI and TII , respectively) and O4FUC (by 0.12 and could facilitate the elimination of the nucleophile protein residue 0.51 Å for TI and TII , respectively). In the hydrolysis reaction, (Bráa et al., 2010). The analyses depicted in Figure 7 show a water molecule (WAT432) establishes a bridge between the this increasing stabilization by Tyr284 as reaction proceeds (by oxygen atom of the attacking water (WAT431) and that of the 4.76 and 6.12 kcal/mol at the hydrolysis and transglycosylation 4-OHFUC group. The OWAT431 -H2WAT432 and the H1WAT432 - TSs and 6.31 and 8.26 kcal/mol at their respective products). O4FUC distances show little variation from reactant to TS, The electrostatic stabilization of the TSs as compared to the but at the product the OWAT431 -H2WAT432 distance increases reactants provided by Arg75, that tends to slightly approach the by 0.91 Å due to a rotation of WAT432 (as OWAT431 has OE1GLU338 as reaction proceeds, is the largest one (by 8.24 and formed the new bond with C1FUC ), which now points its 8.06 kcal/mol for hydrolysis and transglycosylation, respectively). H2WAT432 toward O5FUC (d(H2WAT432 -O5FUC ) is 3.58 Å at These results qualitatively agree with the much lower activity reactants, 3.33 Å at TS and 1.95 Å at products). Thus, WAT432 and 3-OHGLC seem to share a similar role in stabilizing the measured when these two residues are mutated. Finally, Asn282 hydrolysis and transglycosylation TS, respectively, via interaction is predicted to have a moderate destabilizing effect, especially for with the oxocarbenium ion. This kind of interactions have the products. This residue interacts with both catalytic Glu338 been identified in retaining glycosyltransferases to help the and Glu164, thus the results indicate that its interaction with reaction, with the difference that due to the different disposition Glu164, which will go from a total −1 au charge to neutral, is the of the substrates they serve leaving group departure (Gómez predominant one. Frontiers in Chemistry | www.frontiersin.org 14 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation FIGURE 7 | Electrostatic Interaction Analysis (in kcal/mol) of residues in the−1 subsite for (A) hydrolysis and both transglycosylation (B) TI and (C) TII frames studied in this work. Calculations were performed on the structures obtained with QM(small)/MM and at the QM = (PBE0/TZVP) level. Values, with respect to reactant, are given at different points along the reaction coordinate (RC), the transition states (TS) and products. Red: HIS 119, dark green: GLN 18, purple: ASN 163, dark blue: ASN 282, orange: TYR 284, light blue: GLU 392, pink: ARG 75, light green: CYS 164. The residues that directly interact with the nucleophile (Tran et al., 2010; Bissaro et al., 2015a,b). In the present system (Tyr284, Asn282, Arg75) or with the acid/base catalyst (Asn163; the acceptor substrate can interact with Trp312. On the other which also interacts with O2FUC and Arg75) are part of hand, mutation of the catalytic acid/base Glu to Asp in NkBgl the ones identified to increase the T/H ratio when properly also catalyzed unforeseen transglycosylation (Jeng et al., 2012). mutated (Teze et al., 2014) Bissaro et al. recently proposed The crystal structure of the mutant revealed a water molecule that mutations impairing the optimized electron displacement close to the shorter Asp chain and the anomeric carbon; thus, system that the enzyme has for catalyzing hydrolysis (which it was proposed that Asp would play its role through this includes the nucleophile and the oxocarbenium TS), results nearby water molecule. A structurally conserved water molecule in the destabilization of the deglycosylation TSs and leads observed in a wild-type enzyme belonging to family GH101, was to accumulation of the glycosyl-enzyme intermediate; such also proposed to enable a Grotthuss proton shuttle between the accumulation has been linked to an increased T/H ratio protein residue carboxylate and the glycosidic oxygen (Gregg (Bissaro et al., 2015a). Hydrophobic interactions (e.g., C-H – π et al., 2015). Still, the particular reasons why all these mutations interactions) augmenting the affinity of the enzyme for glycoside- in the −1 subsite favor transglycosylation over hydrolysis have like acceptors is another factor that can favor transglycosylation. not been fully unveiled. Frontiers in Chemistry | www.frontiersin.org 15 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation The small differences in energy barriers these mutagenesis energy barrier is predicted to be 1.2–2.6 kcal/mol higher experiments imply and the high dimensionality of enzymatic than the hydrolysis one. Hydrolysis and transglycosylation systems, make it a great challenge to reproduce the effect of transition states are both very loose, with strong oxocarbenium- mutation computationally or to clearly identify its origin. For like character and the proton not yet transferred. Their ring example, in the case of Tyr284 a transglycosylation yield of 76 conformation falls between 4 H3 and 4 E and they present very % (but a lower rate) was measured for the Tyr284Phe mutant as similar atomic charges, except for the more negative charge compared to a 36 % for the WT enzyme (Teze et al., 2014). This of the oxygen atom of the attacking water when compared to was interpreted as a destabilizing effect produced by the mutation that of the attacking 4-OH group of Glc. Structural differences that affected more hydrolysis than transglycosylation. Our appear between the glycosylation and the deglycosylation steps, qualitative analysis of electrostatic interactions, though, estimates with the former having a more pronounced SN 2 character. The that in WT Ttβ-gly Tyr284 provides a slightly higher stabilization proton transfer is also more advanced in the deglycosylation in the transglycosylation case. Attempts to reproduce the TSs that was the proton transfer in the glycosylation one. experimental observations by directly introducing the Tyr284Phe Differences between the hydrolysis and the transglycosylation mutation on frames H, TI , and TII and recalculating the TSs are also observed, but they are very dependent on the level potential energy profiles also failed. Thus, further computational of theory used. At the QM(PBE0/TZVP)/MM level, the potential studies of mutants taking into account structural rearrangements energy maxima for transglycosylation is more advanced on upon mutation (even if small) would be needed to try to the reaction coordinate than the hydrolysis one: has shorter reproduce these trends. However, and more importantly, these O4GLC -C1FUC (forming bond) distance and longer OE2GLU338 - analyses are not able to capture the subtle changes in energy C1FUC (breaking) distance, although the HAcc proton is closer barriers involved in these cases, which are below the current to the Glu164 base in hydrolysis. QM(SCC-DFTB)/MM free “chemical accuracy.” energy maxima show the inverted situation, although large geometric fluctuations are seen for the hydrolysis TS, indicating CONCLUSIONS a large ensemble of TS configurations covering a range of reaction coordinate values. Interestingly, some interactions Here we have presented a QM(DFT)/MM study of the between the substrates (3-OHGLC group) or with neighboring glycosylation, hydrolysis and transglycosylation steps catalyzed water molecules (WAT432) are identified that stabilize the by wild type Thermus thermophilus β-glycosidase (family GH1), a oxocarbenium TSs through interaction with O5FUC and O4FUC . retaining glycosyl hydrolase for which a transglycosylation yield An analysis of electrostatic interactions of the substrates with of 36 % has been determined experimentally for pNP-Fuc and the residues at the−1 subsite has also been performed. However, BnON(Me)-Glc substrates. correlation with experimental mutagenesis studies found to alter In the glycosylation step (first step of the double-displacement the T/H ratio is very challenging due to the small differences mechanism of retaining GHs), the pNP-Fuc substrate establishes on energy barriers involved and potential structural (even if a strong network of hydrogen bonds at the−1 subsite and small) rearrangements upon mutation. The analysis suggests, is found to be bound in an only slightly distorted 4 C1 though, that perturbing the Glu392-Fuc interaction could conformation, contrary to what has been reported when the increase the T/H ratio, even by direct mutation of this residue substrate is a Glc disaccharide. At the transition state for or indirectly as reported experimentally in the N390I and this step, that has a strong oxocarbenium character (close to F401S cases. 4 H conformation but with some 4 E characteristics), proton 3 transfer from Glu164 to the leaving group has not yet started AUTHOR CONTRIBUTIONS but it takes place at a OE2GLU338 -C1FUC distance of ∼2 Å with the QM(large)/MM partition and the PBE0/TZVP QM SR-T performed the calculations. SR-T, ÀG-L, and LM wrote the level. The potential energy barrier at this level of theory manuscript. SR-T, JL, ÀG-L, and LM contributed to the analysis (20.4–22.0 kcal/mol) is in qualitative agreement with the and discussion of the results. experimental barrier (17.1 kcal/mol). The use of a smaller QM region increases this value by more than 5 kcal/mol. ACKNOWLEDGMENTS The size of the QM region also changes the location of the TS along the reaction coordinate. Thus, it seems to affect The authors acknowledge the financial support from the Spanish the nucleophilic and acid/base characteristics of Glu338 and Ministerio de Ciencia e Innovación through contract number Glu164, respectively. CTQ2017-83745-P. LM thanks the Universitat Autònoma de For the hydrolysis and transglycosylation reactions Barcelona for the Talent program. (competing reactions in the deglycosylation step), a significant reduction of the barrier heights is also observed with the SUPPLEMENTARY MATERIAL QM(large)/MM partition, although the T vs. H trends are well reproduced with the QM(small)/MM one provided that The Supplementary Material for this article can be found the PBE0/TZVP level is used to obtain the energies (or also online at: https://www.frontiersin.org/articles/10.3389/fchem. for geometry optimization). The transglycosylation potential 2019.00200/full#supplementary-material Frontiers in Chemistry | www.frontiersin.org 16 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation REFERENCES Cremer, D., and Pople, J. A. (1975). General definition of ring puckering coordinates. J. Am. Chem. Soc. 97, 1354–1358. doi: 10.1021/ja00839a011 Adamo, C., and Barone, V. (1999). Toward reliable density functional methods Danby, P. M., and Withers, S. G. (2016). Advances in enzymatic glycoside without adjustable parameters: the PBE0 model. J. Chem. Phys. 110, 6158–6170. synthesis. ACS Chem. Biol. 11, 1784–1794. doi: 10.1021/acschembio.6b00340 doi: 10.1063/1.478522 David, B., Irague, R., Jouanneau, D., Daligault, F., Czjzek, M., Sanejouand, Y.-H., Ahlrichs, R., Bär, M., Häser, M., Horn, H., and Kölmel, C. (1989). Electronic et al. (2017). Internal water dynamics control the transglycosylation/hydrolysis structure calculations on workstation computers: the program system balance in the agarase (AgaD) of Zobellia galactanivorans. ACS Catal. 7, turbomole. Chem. Phys. Lett. 162, 165–169. doi: 10.1016/0009-2614(89)85118-8 3357–3367. doi: 10.1021/acscatal.7b00348 Albesa-Jové, D., Mendoza, F., Rodrigo-Unzueta, A., Gomollón-Bel, F., Cifuente, J. Davies, G. J., Planas, A., and Rovira, C. (2012). Conformational analyses O., Urresti, S., et al. (2015). A native ternary complex trapped in a crystal reveals of the reaction coordinate of glycosidases. Acc. Chem. Res. 45, 308–316. the catalytic mechanism of a retaining glycosyltransferase. Angew. Chemie Int. doi: 10.1021/ar2001765 Ed. 54, 9898–9902. doi: 10.1002/anie.201504617 Dion, M., Fourage, L., Hallet, J.-N., and Colas, B. (1999). Cloning and Ardèvol, A., and Rovira, C. (2015). Reaction mechanisms in carbohydrate-active expression of a β-glycosidase gene from Thermus thermophilus. Sequence and enzymes: glycoside hydrolases and glycosyltransferases. Insights from ab Initio biochemical characterization of the encoded enzyme. Glycoconj. J. 16, 27–37. quantum mechanics/molecular mechanics dynamic simulations. J. Am. Chem. doi: 10.1023/A:1006997602727 Soc. 137, 7528–7547. doi: 10.1021/jacs.5b01156 Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G. Badieyan, S., Bevan, D. R., and Zhang, C. (2012). Probing the active site chemistry (1995). A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577–8593. of β-glucosidases along the hydrolysis reaction pathway. Biochemistry 51, doi: 10.1063/1.470117 8907–8918. doi: 10.1021/bi300675x Faijes, M., and Planas, A. (2007). In vitro synthesis of artificial polysaccharides Bakowies, D., and Thiel, W. (1996). Hybrid models for combined quantum by glycosidases and glycosynthases. Carbohydr. Res. 342, 1581–1594. mechanical and molecular mechanical approaches. J. Phys. Chem. 100, doi: 10.1016/j.carres.2007.06.015 10580–10594. doi: 10.1021/jp9536514 Feng, H.-Y., Drone, J., Hoffmann, L., Tran, V., Tellier, C., Rabiller, C., Berendsen, H. J. C., Postma, J. P. M., Van Gunsteren, W. F., Dinola, A., and Haak, et al. (2005). Converting a β-Glycosidase into a β-Transglycosidase by J. R. (1984). Molecular dynamics with coupling to an external bath. J. Chem. directed evolution. J. Biol. Chem. 280, 37088–37097. doi: 10.1074/jbc.M5028 Phys. 81, 3684–3690. doi: 10.1063/1.448118 73200 Biarnés, X., Ardèvol, A., Iglesias-Fernández, J., Planas, A., and Rovira, C. Field, R. A. (2011). Challenging reaction equilibria. Nat. Chem. Biol. 7, 658–659. (2011). Catalytic itinerary in 1,3-1,4-β-glucanase unraveled by QM/MM doi: 10.1038/nchembio.668 metadynamics. Charge is not yet fully developed at the oxocarbenium ion-like Geronimo, I., Payne, C. M., and Sandgren, M. (2018). Hydrolysis and transition state. J. Am. Chem. Soc. 133, 20301–20309. doi: 10.1021/ja207113e transglycosylation transition states of glycoside hydrolase family Billeter, S. R., Turner, A. J., and Thiel, W. (2000). Linear scaling geometry 3 β-glucosidases differ in charge and puckering conformation. optimisation and transition state search in hybrid delocalised internal J. Phys. Chem. B 122, 9452–9459. doi: 10.1021/acs.jpcb.8b coordinates. Phys. Chem. Chem. Phys. 2, 2177–2186. doi: 10.1039/a909486e 07118 Bissaro, B., Durand, J., Biarnés, X., Planas, A., Monsan, P., O’Donohue, M. J., Gómez, H., Lluch, J. M., and Masgrau, L. (2013). Substrate- et al. (2015a). Molecular design of non-leloir furanose-transferring enzymes assisted and nucleophilically assisted catalysis in bovine α1,3- from an α- l -arabinofuranosidase: a rationale for the engineering of evolved galactosyltransferase. mechanistic implications for retaining transglycosylases. ACS Catal. 5, 4598–4611. doi: 10.1021/acscatal.5b00949 glycosyltransferases. J. Am. Chem. Soc. 135, 7053–7063. doi: 10.1021/ja4 Bissaro, B., Monsan, P., Fauré, R., and O’Donohue, M. J. (2015b). Glycosynthesis 024447 in a waterworld: new insight into the molecular basis of transglycosylation Gómez, H., Mendoza, F., Lluch, J. M., and Masgrau, L. (2015). QM/MM in retaining glycoside hydrolases. Biochem. J. 467, 17–35. doi: 10.1042/BJ201 studies reveal how substrate–substrate and enzyme–substrate interactions 41412 modulate retaining glycosyltransferases catalysis and mechanism. Adv. Bráa, N. F., Fernandes, P. A., and Ramos, M. J. (2010). QM/MM studies on the β- Protein Chem. Struct. Biol. 100, 225–254. doi: 10.1016/bs.apcsb.2015. galactosidase catalytic mechanism: hydrolysis and transglycosylation reactions. 06.004 J. Chem. Theory Comput. 6, 421–433. doi: 10.1021/ct900530f Gómez, H., Polyak, I., Thiel, W., Lluch, J. M., and Masgrau, L. (2012). Camargo, H. A., Caputa, P., Das, D., Heller, M. P., and Jefferson, R. (2019). Retaining glycosyltransferase mechanism studied by QM/MM methods: Complexity as a novel probe of quantum quenches: universal scalings and lipopolysaccharyl-α-1,4-galactosyltransferase C transfers α-galactose via an purifications. Phys. Rev. Lett. 122:081601. doi: 10.1103/PhysRevLett.122.08 oxocarbenium ion-like transition state. J. Am. Chem. Soc. 134, 4743–4752. 1601 doi: 10.1021/ja210490f Cantarel, B. L., Coutinho, P. M., Rancurel, C., Bernard, T., Lombard, V., and Götz, A. W., Williamson, M. J., Xu, D., Poole, D., Le Grand, S., and Walker, Henrissat, B. (2009). The Carbohydrate-Active EnZymes database (CAZy): R. C. (2012). Routine microsecond molecular dynamics simulations with an expert resource for glycogenomics. Nucleic Acids Res. 37, D233–D238. AMBER on GPUs. 1. generalized born. J. Chem. Theory Comput. 8, 1542–1555. doi: 10.1093/nar/gkn663 doi: 10.1021/ct200909j Case, D. A., Cheatham, T. E. III., Darden, T., Gohlke, H., Luo, R., Merz, K. M. Jr, Gregg, K. J., Suits, M. D. L., Deng, L., Vocadlo, D. J., and Boraston, et al. (2005). The Amber biomolecular simulation programs. J. Comput. Chem. A. B. (2015). Structural analysis of a family 101 glycoside 26, 1668–1688. doi: 10.1002/jcc.20290 hydrolase in complex with carbohydrates reveals insights into its Chuenchor, W., Pengthaisong, S., Robinson, R. C., Yuvaniyama, J., Svasti, J., mechanism. J. Biol. Chem. 290, 25657–25669. doi: 10.1074/jbc.M115.6 and Ketudat Cairns, J. R. (2011). The structural basis of oligosaccharide 80470 binding by rice BGlu1 beta-glucosidase. J.Struct.Biol. 173, 169–179. Hart, G. W., and Copeland, R. J. (2010). Glycomics hits the big time. Cell 143, doi: 10.1016/j.jsb.2010.09.021 672–676. doi: 10.1016/j.cell.2010.11.008 Claeyssens, F., Harvey, J. N., Manby, F. R., Mata, R. A., Mulholland, Henrissat, B. (1991). A classification of glycosyl hydrolases based on amino acid A. J., Ranaghan, K. E., et al. (2006). High-accuracy computation of sequence similarities. Biochem. J. 280(Pt. 2), 309–316. doi: 10.1042/bj2800309 reaction barriers in enzymes. Angew. Chemie Int. Ed. 45, 6856–6859. Hoover, W. G. (1985). Canonical dynamics: equilibrium phase-space distributions. doi: 10.1002/anie.200602711 Phys. Rev. A 31, 1695–1697. doi: 10.1103/PhysRevA.31.1695 Cobucci-Ponzano, B., and Moracci, M. (2012). Glycosynthases as tools for the Hrmova, M., Varghese, J. N., De Gori, R., Smith, B. J., Driguez, H., and production of glycan analogs of natural products. Nat. Prod. Rep. 29:697. Fincher, G. B. (2001). Catalytic mechanisms and reaction intermediates along doi: 10.1039/c2np20032e the hydrolytic pathway of a plant β-D-glucan glucohydrolase. Structure 9, Cobucci-Ponzano, B., Strazzulli, A., Rossi, M., and Moracci, M. (2011). 1005–1016. doi: 10.1016/S0969-2126(01)00673-6 Glycosynthases in biocatalysis. Adv. Synth. Catal. 353, 2284–2300. Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: visual molecular doi: 10.1002/adsc.201100461 dynamics. J. Mol. Graph. 14, 33–38. doi: 10.1016/0263-7855(96)00018-5 Frontiers in Chemistry | www.frontiersin.org 17 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation Jeng, W.-Y., Wang, N.-C., Lin, C.-T., Chang, W.-J., Liu, C.-I., and Wang, A. H.- Petersen, L., Arde, A., Rovira, C., and Reilly, P. J. (2010). Molecular mechanism J. (2012). High-resolution structures of Neotermes koshunensis β-glucosidase of the glycosylation step catalyzed by golgi. J. Am. Chem. Soc. 132, 8291–8300. mutants provide insights into the catalytic mechanism and the synthesis doi: 10.1021/ja909249u of glucoconjugates. Acta Crystallogr. Sect. D Biol. Crystallogr. 68, 829–838. Petersen, L., Ardèvol, A., Rovira, C., and Reilly, P. J. (2009). Mechanism of cellulose doi: 10.1107/S0907444912013224 hydrolysis by inverting GH8 endoglucanases: A QM/MM Metadynamics Study. Kästner, J., and Thiel, W. (2005). Bridging the gap between thermodynamic J. Phys. Chem. B 113, 7331–7339. doi: 10.1021/jp811470d integration and umbrella sampling provides a novel analysis method: Planas, A., and Faijes, M. (2002). Glycosidases and Glycosynthases in enzymatic “Umbrella integration.” J. Chem. Phys. 123:144104. doi: 10.1063/1.20 synthesis of oligosaccharides. An overview. Afinidad 59, 295–313. 52648 Plou, F. J., de Segura, A. G., and Ballesteros, A. (2007). “Application of Glycosidases Kästner, J., and Thiel, W. (2006). Analysis of the statistical error in umbrella and Transglycosidases in the Synthesis of Oligosaccharides,” in Industrial sampling simulations by umbrella integration. J. Chem. Phys. 124:234106. Enzymes (Dordrecht: Springer), 141–157. doi: 10.1063/1.2206775 Raich, L., Borodkin, V., Fang, W., Castro-López, J., Van Aalten, D. M. F., Hurtado- Kempton, J. B., and Withers, S. G. (1992). Mechanism of Guerrero, R., et al. (2016). A trapped covalent intermediate of a glycoside Agrobacterium.beta.-glucosidase: kinetic studies. Biochemistry 31, 9961–9969. hydrolase on the pathway to transglycosylation. insights from experiments and doi: 10.1021/bi00156a015 quantum mechanics/molecular mechanics simulations. J. Am. Chem. Soc. 138, Kiessling, L. L., and Splain, R. A. (2010). Chemical approaches 3325–3332. doi: 10.1021/jacs.5b10092 to glycobiology. Annu. Rev. Biochem. 79, 619–653. Reed, A. E., Weinstock, R. B., and Weinhold, F. (1985). Natural population doi: 10.1146/annurev.biochem.77.070606.100917 analysis. J. Chem. Phys. 83, 735–746. doi: 10.1063/1.449486 Kirschner, K. N., Yongye, A. B., Tschampel, S. M., González-Outeiriño, J., Daniels, Roe, D. R., and Cheatham, T. E. (2013). PTRAJ and CPPTRAJ: software for C. R., Foley, B. L., et al. (2008). GLYCAM06: a generalizable biomolecular processing and analysis of molecular dynamics trajectory data. J. Chem. Theory force field. Carbohydrates. J. Comput. Chem. 29, 622–655. doi: 10.1002/jcc. Comput. 9, 3084–3095. doi: 10.1021/ct400341p 20820 Ryckaert, J. P., Ciccotti, G., and Berendsen, H. J. C. (1977). Numerical Liu, D. C., and Nocedal, J. (1989). On the limited memory BFGS method for integration of the cartesian equations of motion of a system with large scale optimization. Math. Program. 45, 503–528. doi: 10.1007/BF015 constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327–341. 89116 doi: 10.1016/0021-9991(77)90098-5 Lokanath, N. K., Shiromizu, I., Miyano, M., Yokoyama, S., Kuramitsu, S., and Sajith, S., Priji, P., Sreedevi, S., and Benjamin, S. (2016). An overview on Kunishima, N. (in press). Structure of Beta-Glucosidase at Atomic Resolution fungal cellulases with an industrial perspective. J. Nutr. Food Sci. 6:461. from Thermus Thermophilus Hb8. doi: 10.4172/2155-9600.1000461 Mackenzie, L. F., Wang, Q., Warren, R. A. J., and Withers, S. G. Salomon-Ferrer, R., Götz, A. W., Poole, D., Le Grand, S., and Walker, R. C. (2013). (1998). Glycosynthases: mutant glycosidases for oligosaccharide Routine microsecond molecular dynamics simulations with AMBER on GPUs. synthesis. J. Am. Chem. Soc. 120, 5583–5584. doi: 10.1021/ja98 2. Explicit solvent particle mesh ewald. J. Chem. Theory Comput. 9, 3878–3888. 0833d doi: 10.1021/ct400314y Maier, J. A., Martinez, C., Kasavajhala, K., Wickstrom, L., Hauser, K. E., and Saumonneau, A., Champion, E., Peltier-Pain, P., Molnar-Gabor, D., Hendrickx, Simmerling, C. (2015). ff14SB: improving the accuracy of protein side chain and J., Tran, V., et al. (2016). Design of an α-l-transfucosidase for the backbone parameters from ff99SB. J. Chem. Theory Comput. 11, 3696–3713. synthesis of fucosylated HMOs. Glycobiology 26, 261–269. doi: 10.1093/glycob/ doi: 10.1021/acs.jctc.5b00255 cwv099 Malet, C., and Planas, A. (1998). From β-glucanase to β-glucansynthase: Seeberger, P. H., and Cummings, R. D. (2017). “Glycans in Biotechnology and glycosyl transfer to α-glycosyl fluorides catalyzed by a mutant the Pharmaceutical Industry,” in Essentials of Glycobiology, eds E. A. Varki, endoglucanase lacking its catalytic nucleophile. FEBS Lett. 440, 208–212. Cummings, R. D., Esko, J. D., et al. (Cold Spring Harbor, NY: Cold Spring doi: 10.1016/S0014-5793(98)01448-3 Harbor Laboratory Press). Mendoza, F., Gómez, H., Lluch, J. M., and Masgrau, L. (2016). α1,4- N - Sherwood, P., de Vries, A. H., Guest, M. F., Schreckenbach, G., Catlow, Acetylhexosaminyltransferase EXTL2: the missing link for understanding C. R. A., French, S. A., et al. (2003). QUASI: A general purpose glycosidic bond biosynthesis with retention of configuration. ACS Catal. 6, implementation of the QM/MM approach and its application to problems 2577–2589. doi: 10.1021/acscatal.5b02945 in catalysis. J. Mol. Struct. Theochem 632, 1–28. doi: 10.1016/S0166-1280(03) Mendoza, F., Lluch, J. M., and Masgrau, L. (2017). Computational insights into 00285-9 active site shaping for substrate and regiospecificity in the EXTL2 retaining Smith, W., and Forester, T. R. (1996). DL_POLY_2.0: a general-purpose glycosyltransferase. Organic Biomol. Chem. doi: 10.1039/C7OB01937H parallel molecular dynamics simulation package. J. Mol. Graph. 14, 136–141. Metz, S., Kästner, J., Sokol, A. A., Keal, T. W., and Sherwood, P. (2014). doi: 10.1016/S0263-7855(96)00043-4 ChemShell—a modular software package for QM/MM simulations. Teze, D. (2012). Recherche des Bases Moléculaires de L’équilibre Transglycosylation Wiley Interdiscip. Rev. Comput. Mol. Sci. 4, 101–110. doi: 10.1002/w / Hydrolyse d’une Glycoside Hydrolase de la Famille 1. Université de Nantes, cms.1163 Faculté des Sciences et des Techniques. 1–198. Moracci, M., Trincone, A., Perugino, G., Ciaramella, M., and Rossi, M. (1998). Teze, D., Daligault, F., Ferrières, V., Sanejouand, Y.-H., and Tellier, C. Restoration of the activity of active-site mutants of the hyperthermophilic (2015). Semi-rational approach for converting a GH36 α-glycosidase β-glycosidase from sulfolobus solfataricus: dependence of the mechanism into an α-transglycosidase. Glycobiology 25, 420–427. doi: 10.1093/glycob/ on the action of external nucleophiles. Biochemistry 37, 17262–17270. cwu124 doi: 10.1021/bi981855f Teze, D., Dion, M., Daligault, F., Tran, V., André-Miral, C., and Tellier, Nosé, S. (1984). A unified formulation of the constant temperature molecular C. (2013). Alkoxyamino glycoside acceptors for the regioselective dynamics methods. J. Chem. Phys. 81, 511–519. doi: 10.1063/1.4 synthesis of oligosaccharides using glycosynthases and transglycosidases. 47334 Bioorg. Med. Chem. Lett. 23, 448-−451. doi: 10.1016/j.bmcl.2012. Olsson, M. H. M., Søndergaard, C. R., Rostkowski, M., and Jensen, J. H. (2011). 11.065 PROPKA3: consistent treatment of internal and surface residues in empirical Teze, D., Hendrickx, J., Czjzek, M., Ropartz, D., Sanejouand, Y.-H., Tran, V., pKa predictions. J. Chem. Theory Comput. 7, 525–537. doi: 10.1021/ct1 et al. (2014). Semi-rational approach for converting a GH1 β-glycosidase into 00578z a β-transglycosidase. Protein Eng. Des. Sel. 27, 13–19. doi: 10.1093/protein/ Pereira, A. T., Ribeiro, A. J. M. M., Fernandes, P. A., and Ramos, gzt057 M. J. (2017). Benchmarking of density functionals for the kinetics Tran, V., Hoffmann, L., Rabiller, C., Tellier, C., and Dion, M. (2010). and thermodynamics of the hydrolysis of glycosidic bonds catalyzed Rational design of a GH1 -glycosidase to prevent self-condensation by glycosidases. Int. J. Quantum Chem. 117:e25409. doi: 10.1002/qua. during the transglycosylation reaction. Protein Eng. Des. Sel. 23, 43–49. 25409 doi: 10.1093/protein/gzp068 Frontiers in Chemistry | www.frontiersin.org 18 April 2019 | Volume 7 | Article 200 Romero-Téllez et al. QM/MM β-Glycosidase: Hydrolysis vs. Transglycosylation Wang, J., Hou, Q., Dong, L., Liu, Y., and Liu, C. (2011). QM/MM studies on the Conflict of Interest Statement: The authors declare that the research was glycosylation mechanism of rice BGlu1 β-glucosidase. J. Mol. Graph. Model. 30, conducted in the absence of any commercial or financial relationships that could 148–152. doi: 10.1016/j.jmgm.2011.06.012 be construed as a potential conflict of interest. Wang, J., Hou, Q., Sheng, X., Gao, J., Liu, Y., and Liu, C. (2013). Theoretical study on the deglycosylation mechanism of rice BGlu1 β-glucosidase. Int. J. Quantum Copyright © 2019 Romero-Téllez, Lluch, González-Lafont and Masgrau. This is an Chem. 113, 1071–1075. doi: 10.1002/qua.24131 open-access article distributed under the terms of the Creative Commons Attribution Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). License (CC BY). The use, distribution or reproduction in other forums is permitted, Development and testing of a general amber force field. J. Comput. Chem. 25, provided the original author(s) and the copyright owner(s) are credited and that the 1157–1174. doi: 10.1002/jcc.20035 original publication in this journal is cited, in accordance with accepted academic Zechel, D. L., and Withers, S. G. (2000). Glycosidase mechanisms : anatomy of a practice. No use, distribution or reproduction is permitted which does not comply finely tuned catalyst. Acc. Chem. Res. 33, 11–18. doi: 10.1021/ar970172+ with these terms. Frontiers in Chemistry | www.frontiersin.org 19 April 2019 | Volume 7 | Article 200