ARTICLE https://doi.org/10.1038/s41467-019-09691-z OPEN Discovery of processive catalysis by an exo-hydrolase with a pocket-shaped active site Victor A. Streltsov1,15, Sukanya Luang2,15, Alys Peisley 1, Joseph N. Varghese1, James R. Ketudat Cairns 3, Sebastien Fort 4, Marcel Hijnen5, Igor Tvaroška6, Ana Ardá7, Jesús Jiménez-Barbero7, Mercedes Alfonso-Prieto 8, Carme Rovira 8,9, Fernanda Mendoza 10,11, Laura Tiessler-Sala11, José-Emilio Sánchez-Aparicio 11, Jaime Rodríguez-Guerra 11,12, José M. Lluch 11,13, Jean-Didier Maréchal11, Laura Masgrau 11,13 & Maria Hrmova 2,14 1234567890():,; Substrates associate and products dissociate from enzyme catalytic sites rapidly, which hampers investigations of their trajectories. The high-resolution structure of the native Hordeum exo-hydrolase HvExoI isolated from seedlings reveals that non-covalently trapped glucose forms a stable enzyme-product complex. Here, we report that the alkyl β-D-glucoside and methyl 6-thio-β-gentiobioside substrate analogues perfused in crystalline HvExoI bind across the catalytic site after they displace glucose, while methyl 2-thio-β-sophoroside attaches nearby. Structural analyses and multi-scale molecular modelling of nanoscale reactant movements in HvExoI reveal that upon productive binding of incoming substrates, the glucose product modifies its binding patterns and evokes the formation of a transient lateral cavity, which serves as a conduit for glucose departure to allow for the next catalytic round. This path enables substrate-product assisted processive catalysis through multiple hydrolytic events without HvExoI losing contact with oligo- or polymeric substrates. We anticipate that such enzyme plasticity could be prevalent among exo-hydrolases. 1 Commonwealth Scientific and Industrial Research Organisation, Materials Science and Engineering, Parkville Victoria 3052, Australia. 2 School of Agriculture, Food and Wine, University of Adelaide, Waite Campus, Glen Osmond South Australia 5064, Australia. 3 School of Chemistry and Center for Biomolecular Structure, Function and Application, Suranaree University of Technology, Nakhon Ratchasima 30000, Thailand. 4 University Grenoble Alpes, Centre de Recherches sur les Macromolécules Végétales, Grenoble cedex 9 38041, France. 5 GE Healthcare Life Sciences, Paramatta NSW 2150, Australia. 6 Institute of Chemistry, Slovak Academy of Sciences, Bratislava 84538, Slovak Republic. 7 Centre for Cooperative Research in Biosciences, Derio-Bizkaia ́ 48160, Spain. 8 Departament de Quimica Inorgànica i Orgànica, Universitat de Barcelona, Barcelona 08028, Spain. 9 Institució Catalana de Recerca i Estudis Avançats, Barcelona 08010, Spain. 10 Departamento de Ciencias Químicas, Universidad Andrés Bello, Sede Concepción, Talcahuano 4260000, Chile. 11 Departament de Química, Universitat Autònoma de Barcelona, Bellaterra 08193, Spain. 12 Institute of Chemical Research of Catalonia, The Barcelona Institute of Science and Technology, Tarragona 43007, Spain. 13 Institut de Biotecnologia i de Biomedicina, Universitat Autònoma de Barcelona, Bellaterra 08193, Spain. 14 School of Life Sciences, Huaiyin Normal University, Huaian 223300, China. 15These authors contributed equally: Victor A. Streltsov and Sukanya Luang. Correspondence and requests for materials should be addressed to M.H. (email: maria.hrmova@adelaide.edu.au) NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z E nzymes are biological catalysts that are fundamental to life. non-reducing termini of oligo- and polymeric β-D-glucosides to Enzymes afford enormous accelerations to chemical reac- support seed germination and plant development9. tions compared to uncatalysed reactions. It has become In this work, we examine product and substrate pathways increasingly recognised that the largest contribution to the along the catalytic cycle of a plant exo-hydrolase HvExoI, and enzyme catalytic power arises from the electrostatic environment how this enzyme utilises its plasticity and that of glycoside of polar active sites1,2. Enzymes use protein architectures to substrates. We use the high-resolution X-ray crystallography precisely position a set of amino acid residues to catalyse inter- supported with enzyme kinetics, mass spectrometry, nuclear conversions of substrates into products. This structure-based view magnetic resonance (NMR) spectroscopy, multi-scale molecular is supported by thermodynamic and kinetic models that treat modelling employing docking19, Molecular dynamics (MD) enzymes in complex with substrates and transition states in simulations, Genetic Algorithms with Unrestricted Descriptors succession3. Only rarely these models consider substrate and for Intuitive Molecular Modelling (GaudiMM)20 and Protein product diffusion as a part of their catalytic mechanisms, as these Energy Landscape Exploration (PELE)21 calculations, to reveal processes proceed rapidly, although on occasions products are the Glc product displacement route, and to extract how each seen entrapped in active sites4. Yet, the lack of descriptions of hydrolytic event including Glc release is precisely coordinated substrate associative and product dissociative pathways creates a with the incoming substrate association and hydrolysis. This leads knowledge gap5. us to generate a variant enzyme and probe it for Glc entrapment A vital conclusion drawn from the structural studies with the and the S-glycoside analogue binding. Along this journey with native β-D-glucan glucohydrolase, isoenzyme ExoI (HvExoI) HvExoI, we observe a remarkable phenomenon, coin the term isolated from barley seedlings, was that the glucose (Glc) product ‘substrate-product assisted processive catalysis’, and describe how released from β-D-glucoside substrates remains entrapped in product and substrate trajectories are coordinated via this the enzyme active site until an incoming substrate binds6–9, mechanism during catalysis. We discover that through chemical presumably lowering the energy barrier to facilitate Glc dis- signalling, the incoming substrate evokes the formation of an placement. At this stage, the mechanism of Glc displacement and autonomous and transient lateral cavity that serves as a conduit how it is linked to the catalytic cycle of HvExoI remained for the Glc product displacement from the active site. Crucially, unanswered. Notably, no other native GH3 structures with these co-operative reactant pathways, where polymeric substrates entrapped Glc are available, because these recombinant enzymes never dissociate from the enzyme, sustain consecutive relocations not exposed to presumably high enough Glc concentration, of substrates and products via dedicated routes, ensuring that this do not entrap Glc during protein maturation and secretion. exo-hydrolase adopts processive catalysis. However, several GH3 enzyme complexes with in crystallo- perfused Glc10,11 or Glc-derivatives12,13 are available. Addition- Results ally, the native GH78 α-L-rhamnosidase with a deep pocket- GC/MS proves that Glc is bound to crystalline native HvExoI. shaped active site also holds the entrapped Glc molecule, which We purified22, crystallised23 and chemically analysed crystals to was not perfused in crystals14. reconcile our structural observations with Glc entrapment in We have previously observed that mechanism-based inhibitors native HvExoI purified from barley seedlings6–9. Several crystals conduritol B epoxide and 2′,4′-dinitrophenyl 2-deoxy-2-fluoro-β- removed from the mother liquor were extracted and fractionated D-glucoside dislodge Glc from the −1 subsite and form stable by normal-phase high performance liquid chromatography cyclitol esters via an α-anomeric linkage with Asp2857. Similarly, (HPLC) with evaporative light scattering detection (Fig. 1a; top transition-state ion-like gluco-phenylimidazole mimics displace panel, solid line). The mother liquor from these crystals was also Glc and bind in the −1 subsite with their Glc component dis- analysed, both without and augmented with authentic Glc torted in the 4E envelope conformation15,16. Non-hydrolysable (Fig. 1a; top panel, dashed and dotted lines, respectively). The S-glycoside analogues also remove Glc from the −1 subsite and materials from these preparations eluting between 6.0 and 7.3 bind across the active site7,17. Because HvExoI hydrolyses various min during HPLC separation were pooled, reduced and per- positional isomers of β-D-glucosides, we rationalised this obser- acetylated. The total ion chromatograms of the material extracted vation by the non-reducing Glc moiety being locked into a fixed from crystals and from the mother liquor augmented with Glc position via hydrogen bonds (H-bonds) with the residues in the showed a peak with the retention time of 23.26 min, corre- −1 subsite, whereas the reducing-end moiety is free to adopt sponding to glucitol hexa-acetate (Fig. 1b; 1st and 3rd middle multiple orientations in the +1 subsite. This allows the remainder panels). Mass fragmentation spectra of these materials (Fig. 1c; of the saccharide substrates beyond the +1 subsite to project bottom panels) showed identical profiles of total ions, demon- unencumbered7,9,17. strating that the chemical substance eluting at 23.26 min was Glc. The fundamental function of hydrolysis of oligo- and poly- We could not detect Glc in the mother liquor, from which crystals saccharides is essential for the understanding of the global carbon had been removed (Fig. 1b; middle panel), while authentic Glc cycle, which forms the basis of a multibillion-dollar biotechno- added to the same mother liquor was readily identifiable. logical industry. In this context, the HvExoI enzyme6 that belongs to the glycoside hydrolase (GH) family 318, serves as one of the archetypal models for hydrolysis of oligo- and polysaccharides. Crystal structures reveal the Glc product in native HvExoI. Structural information of GH3 enzymes is available through 36 The key observation that led us to explore the precise spatial unique entries in the Protein Data Bank (PDB) from all phyla, disposition of the Glc product in the active site and the Glc although there is only one plant structure available (that of displacement route in HvExoI, was based on our original obser- HvExoI), notably in complex with a variety of oligosaccharide vation that the Glc product remains entrapped in the enzyme’s substrate analogues and mimics7–9. HvExoI folds in an (α/β)8 active site that is isolated from the native source of young barrel (domain 1) and an (α/β)6 sandwich (domain 2) that at their seedlings6,8. We anticipated that the last remaining Glc molecule interface cradle a 13 Å-deep pocket- or crater-shaped active site. originating from oligo- and polysaccharide substrates, stays This site accommodates the catalytic nucleophile Asp285 and associated with the enzyme in plant tissues. To this end, previous catalytic acid/base Glu491 in the −1 subsite, while the aromatic crystal structures of native HvExoI at 2.20 Å and 2.70 Å (in-house clamp of Trp286 and Trp434 indole moieties constitutes the X-ray source data) resolution showed a single Glc molecule neighbouring +1 subsite6–9. HvExoI releases a single Glc from bound in the −1 subsite6,8. In the current work we refined native 2 NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z ARTICLE a 24 22 ELSD (mV) 20 Pooled 18 Authentic glucose 16 1 2 3 4 5 6 7 min b 20,000 17.01 15,000 Exogenous myo-inositol 10,000 22.54 23.26 5000 * 4000 Abundance 3000 Myo-inositol 2000 17.01 22.54 1000 4000 3000 Glucitol hexa-acetate 2000 17.01 23.26 * 1000 4 8 12 16 20 24 28 min c 800 115 600 103 128 145 157 400 187 171 217 200 260 Abundance 290 100 120 140 160 180 200 220 240 260 280 300 320 340 115 2000 1600 128 145 103 187 157 1200 171 218 800 289 260 400 220 230 273 100 120 140 160 180 200 220 240 260 280 300 320 340 m/z 4C Fig. 1 Glc trapped in native HvExoI is bound non-covalently in the 1 conformation. a HPLC chromatogram of material extracted from crystals (solid line), mother liquor (dashed) and mother liquor augmented with authentic Glc (dotted). Arrow points to the peak of authentic Glc in augmented mother liquor. Fractions (left right arrow) were pooled for GC/MS analysis. b Total ion chromatograms obtained by GC/MS of HPLC-eluted materials (6–7.3 min) containing alditol acetates. Materials from crystals (top), mother liquor (middle) and mother liquor augmented with Glc (bottom). Myo-inositol (20 ng) served as an internal standard (top and middle panels). Numbers near peaks indicate retention times. Substance eluted in 17.01 min peak corresponds to a contaminating plasticiser. Material extracted from crystals and eluting at 23.26 min (top panel) is identical to that of glucitol hexa-acetate at 23.26 min (bottom panel) on a Prevail Carbohydrate ES column (peaks indicated by asterisks). c Fragmentation mass spectra of glucitol hexa-acetate: top, material extracted from crystals; bottom glucitol hexa-acetate eluting at 23.26 min NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z HvExoI at 1.65 Å resolution (synchrotron data), revealing that a Native with trapped Glc Glc is bound at 0.5 occupancy at each −1 and +1 subsites in the Tyr253 Tyr253 4C chair conformation, suggesting that it may be mobile between 1 Trp286 Trp286 both subsites (Fig. 2a; Supplementary Fig. 1a). Glu491 Glc′ Glu491 Asp285 Asp285 Arg158 Arg158 His207 His207 Glc is absent in recombinant HvExoI but could be perfused in. Although Glc has always been detected in the native enzyme Trp434 Trp434 isolated from barley seedlings6,8, it was uncertain if it could be Lys206 Glc Gly57 Lys206 Gly57 observed in recombinant HvExoI produced in Pichia pastoris. –1 subsite +1 subsite Asp95 Asp95 The 1.45 Å structure of the recombinant enzyme, which is kine- tically and structurally nearly identical to the native one24,25, b Native perfused with 3dGlc failed to show the presence of Glc in the active site, where we identified glycerol and up to seven water molecules mimicking Trp286 Trp286 the positions of hydroxyl groups of Glc bound in native HvExoI Tyr253 Tyr253 (Fig. 3a; Supplementary Fig. 2a). However, after the crystals of Glu491 Asp285 3dGlc Glu491 Asp285 recombinant HvExoI were perfused with Glc at a near saturating Trp434 Trp434 concentration, the 1.55 Å structure revealed two Glc molecules Gly57 Gly57 one each in the −1 and +1 subsites, adopting alternate 4C1 chair Arg158 Arg158 (occupancy 0.8) and 1S3 skew boat (occupancy 0.2) conforma- His207 Glc +1 subsite His207 tions with classical Cremer–Pople ring-puckering parameters26. Lys206 Asp95 Lys206 Asp95 –1 subsite In the −1 subsite, we detected the network of 12–13 mono- and co-operative bi- and tri-dentate H-bonds27,28; of interest were c Native perfused with octyl-O-Glc those forming short H-bonds of 2.5–2.7 Å between Oδ1 and Oδ2 of Asp95, and the C6-OH and C4-OH groups of the Glc moiety. We attempted to position into the electron density other skew Octyl-O-Glc Trp286 Trp286 boat conformers of Glc (3S5, 1S5) with various occupancies, Asp285 Asp285 however, convergence in the refinement was reached with those His207 His207 Glu491 Glu491 of 4C1 and 1S3. We concluded that during maturation in P. pas- toris, Glc cannot be entrapped and thus observed, because the Lys206 Lys206 intracellular concentration in this host is not high enough during Trp434 Trp434 –1 subsite +1 subsite enzyme maturation or secretion. To demonstrate that Glc was Asp95 Asp95 not displaced from the active site by the cryoprotectant, data were collected from recombinant crystals at ambient temperatures d Native perfused with PEG (n = 5–10) without applying vitrification solutions. Here, no Glc was Tyr253 Tyr253 Trp286 Trp286 observed in the active site, but was readily observable after per- Asp285 Asp285 fusion of crystals with Glc. His207 PEG 1/1′ His207 Glc Lys206 Lys206 SPR analysis of Glc and substrate binding to HvExoI. To Lys206 Glu491 Glu491 determine the strength of Glc binding and to compare it with that Asp95 Gly57 Asp95 Gly57 of the substrate thio-analogue methyl 6-thio-β-gentiobioside Arg158 Trp434 PEG 2 Arg158 Trp434 (G6SG-OMe), we used the recombinant HvExoI enzyme. The –1 subsite Gln435 Gln435 steady-state affinity29 KD value of 0.008 × 10–3 M for G6SG-OMe +1 subsite revealed that it bound tighter than the Glc hydrolytic product Fig. 2 Native HvExoI with trapped Glc, and with perfused 3dGlc, octyl-O- (KD of 0.16 × 10–3 M) (Table 1; Supplementary Fig. 3). Glc and PEG. a Stereo view of native HvExoI with Glc (Glc and Glc′, carbons: yellow sticks) at 0.5 occupancy, which oscillates between the −1 NMR spectroscopy of Glc binding to recombinant HvExoI. To and +1 subsites. b Stereo view of native HvExoI with Glc (carbons: yellow identify conformational states of bound Glc with the measured sticks) and 3dGlc (carbons: orange sticks) at 1.0 occupancies, bound in the KD value of 0.16 × 10–3 M (Table 1; Supplementary Fig. 3), we −1 and +1 subsites, respectively. c Stereo view of native HvExoI with octyl- used 1H saturation transferred difference (STD)30 and transferred O-Glc (carbons: orange sticks) at 1.0 occupancy, bound across the −1 and nuclear Overhauser effect spectroscopy (trNOESY)31 in solution +1 subsites. d Stereo view of native HvExoI with two PEG (n = 5–10) NMR. STD spectra recorded with Glc in a 60-fold excess relative molecules (carbons: orange sticks; PEG 1 in two alternate conformations to HvExoI detected Glc bound non-covalently. As only STDs for at occupancies 0.5 each and PEG 2 at occupancy 1.0), bound in the +1 and protons of β-anomeric Glc (and not for α-form) were observed, putative +2 subsites. Grey, red, and blue represent carbon, oxygen and while one would expect an α:β anomeric ratio of ~2:3 in solution nitrogen atoms, respectively. Water molecules are shown as red or equilibrium, we concluded that the enzyme specifically binds magenta (bound alternate water molecules when the ligand is missing and β-anomeric Glc (Fig. 4a); this observation supports the retaining are not numbered) spheres in complexes with Glc and PEG. In the PEG hydrolytic mechanism of HvExoI22. To identify the conformation complex magenta-coloured water molecules mimic positions of OH groups of Glc, trNOESY experiments under a relatively low HvExoI of Glc (added for illustration in magenta lines). Separations of less than versus Glc ratio revealed weak trNOEs between H1-H3 and 3.50 Å from active site residues (carbons: grey sticks) are shown as H1-H5 protons of β-anomeric Glc (Fig. 4b), consistent with dashed lines. Derived |2mFobs − DFcalc| electron density maps are contoured β-anomeric Glc bound in the low-energy 4C1 chair conformation. at 1σ (blue mesh). In |2mFobs − DFcalc|, Fobs and Fcalc are observed and To corroborate these findings, thiocellobiose which cannot be calculated X-ray structure factor amplitudes, where m is the figure of merit hydrolysed by HvExoI7, was added in excess relative to HvExoI. and D is the estimated coordinate error Here, clear negative NOEs between specific thiocellobiose proton pairs were observed (Fig. 4c), while strong trNOEs of the 4 NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z ARTICLE a Recominant in ligand-free form Table 1 Binding (KD) and inhibition (Ki) parameters for Tyr253 Tyr253 HvExoI by ligands and substrate analogues Gol 2/2′/2″ Trp286 Trp286 Asp285 Asp285 Ligand/ KDa/Kib (M) ΔG c Chemical His207 His207 inhibitor (kJ × mol−1) formulad Glca 0.16 × 10–3 −21.7 C6H12O6 Gol 1 Lys206 Lys206 G6SG-OMea 0.008 × 10–3 −29.1 C13H24O10S Arg158 –1 subsite Gly57 Arg158 Gly57 3dGlcb 9.8 × 10–3 −11.7 C6H12O5 9.4 × 10–3 Trp434 Trp434 Asp95 +1 subsite Asp95 4dGlcb −11.8 C6H12O5 Octyl-O-Glcb 0.13 × 10–3 −22.6 C14H28O6 b Recombinant perfused with Glc Octyl-S-Glcb 1.1 × 10–3 −17.2 C14H28O5S G2SG-OMeb 2.55 × 10–3 −15.1 C13H24O10S Tyr253 Tyr253 Glc 2/2′ aUsing Surface Plasmon Resonance with recombinant HvExoI Glu491 Glu491 bUsing inhibition kinetics with the 4-nitrophenyl β-D-glucopyranoside substrate8 with Trp286 Trp286 Asp285 Asp285 native HvExoI cCalculated according to ΔG = −RT ln [K ] or ΔG = −RT ln [K ]29 Arg158 Arg158 i D dChemical structures are shown in Supplementary Table 1 His207 His207 Trp434 Trp434 Glc 1/1′ Gly57 Gly57 Lys206 Lys206 +1 subsite –1 subsite Conformational Free Energy Landscape of Glc bound to Asp95 Asp95 HvExoI. To further assess conformational states of bound Glc in the high-resolution structure of recombinant HvExoI in complex c Recombinant perfused with G6SG-OMe with Glc, quantum mechanics/molecular mechanics (QM/MM) metadynamics simulations were performed to asses Conforma- His207 Arg158 Asp285 Trp286 Arg158 His207 Asp285 Trp286 tional Free Energy Landscape (FEL) maps of β-D-Glc bound in the active site of HvExoI (Fig. 5; Supplementary Fig. 2b). These Lys206 Lys206 analyses indicated that although Glc was engaged in extensive Trp434 Trp434 H-bond networks, the β-D-glucopyranose ring adopted 4C1 chair G6SG-OMe Gly57 Gly57 and 1S3/B3,O (−1 subsite) or 4C1 chair and B3,O (+1 subsite) Asp95 +1 subsite Asp95 distorted skew boat/boat conformations (Fig. 5; left and middle –1 subsite panels). More precisely in the −1 subsite, H-bonds with Asp285 and Arg158 restrained the 2-OH group in the equatorial position, d Recombinant perfused with G2SG-OMe and stabilised the 4C1 conformation by 8 kcal/mol more than that Glu287 Glu287 of 1S3 or B3,O, with a conformational energy barrier change of G2SG-OMe 11 kcal/mol. Conversely, hydroxyl groups of Glc in the +1 subsite interacted with Asp95, Arg158, Lys206, His207, Asp285, and Trp286 Trp286 Glu491, stabilising the 4C1 Glc conformation by 2 kcal/mol more Glu491 Glu491 PEG 1 than that of B3,O, with a conformational energy barrier change of 6 kcal/mol. Notably, these findings concurred with conforma- Glu491 Glu491 tional FEL of isolated Glc33, in which 4C1 and 1S3/B3,O con- Gol 1 Trp434 +2 subsite Trp434 formers were most stable (Fig. 5; right panel). Nonetheless, FEL –1 subsite Gly57 Gly57 of Glc bound in the −1 and +1 subsites exhibited exclusively +1 subsite these two minima due to a restricted conformational space within the enzyme active site (Figs. 2a, 3b, 5) as demonstrated in Fig. 3 Recombinant HvExoI and with perfused Glc, G6SG-OMe, or G2SG- other glycosyl hydrolases34. OMe. a Stereo view of recombinant HvExoI in a ligand-free form. Two glycerol molecules (carbons: green sticks) (Gol 1 at occupancy 0.5; Gol 2 in three alternate conformations) are bound in the −1 and +1 subsites. Inhibition kinetics of substrate analogues in native HvExoI. To b Stereo view of recombinant HvExoI with two Glc molecules in the 4C1 shed light on the molecular basis of the Glc displacement route, conformation (carbons: yellow sticks) at 0.8 occupancy and in the 1S3 we selected a series of deoxy-Glc (dGlc) and alkyl-glucoside conformation (carbons: cyan sticks) at 0.2 occupancy, bound in the −1 derivatives, and substrate mimics that could potentially displace and +1 subsites. c Stereo view of recombinant HvExoI with G6SG-OMe Glc6. To assess the strength of inhibition, we determined (carbons: orange sticks) at 1.0 occupancy, bound across the −1 and +1. dissociation constants of enzyme-inhibitor complexes (Ki) for d Stereo view of recombinant HvExoI with G2SG-OMe (carbons: orange 3-deoxy-glucose (3dGlc), 4-deoxy-glucose (4dGlc), n-octyl β-D- sticks) at 0.7 occupancy, bound across the +1 and putative +2 subsites. glucopyranoside (octyl-O-Glc), n-octyl 1-thio-β-D-glucopyrano- Water molecules are shown as red spheres. Separations of less than 3.50 Å side (octyl-S-Glc), a hydrophilic polymer polyethylene glycol from the active site residues (carbons: grey sticks) are shown as dashed (PEG; putative fractional polysaccharide mimic35) and the thio- lines. Derived |2mFobs − DFcalc| electron density maps are contoured at analogue methyl 2-thio-β-sophoroside (G2SG-OMe). All inhibi- 1σ for G2SG-OMe and Gol (blue mesh) in (a–c). In (d), electron density tors showed competitive inhibition with Ki values between 9.8 × maps are contoured at 1σ for PEG (orange mesh) and Gol (green mesh) 10–3 M and 0.13 × 10–3 M (Table 1). The most effective was octyl- in the −1 subsite O-Glc with Ki of 0.13 × 10–3 M, while surprisingly G2SG-OMe was weak with Ki of 2.55 × 10–3 M, while PEG was non-inhibi- tory, as expected. anomeric proton of non-reducing Glc (H-1′) with intra-residual H-3′ and H-5′ (Supplementary Fig. 4; Supplementary Note 1) Glc is not removed from native HvExoI by dGlc derivatives. To implied that the non-reducing β-D-glucopyranose ring of thio- test the hypothetical Glc displacement route by using dGlc deri- cellobiose was bound to HvExoI in the 4C1 conformation32. vatives, crystals were perfused with saturating concentrations of NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z a STD H6′β (×50) H5β + H3β H4β + H4α H6α + H5α H3α H6′α H2α H1β H6β H2β Reference 4.5 4.0 3.5 3.0 2.5 ppm b 4.40 4.35 4.30 F1 [ppm] H1/H5 H1/H3 3.8 3.6 3.4 3.2 3.0 2.8 F2 [ppm] c 4.3 H1/H5 H1/H3 H1β H1′ 4.4 4.5 H1′/H3′ H1′/H5′ H1′/H2′ H1′/H4 3.8 3.6 3.4 3.2 3.0 2.8 F2 [ppm] Fig. 4 Recombinant HvExoI recognises β-D-Glc in the 1 conformation. a 4C 1H STD NMR spectrum (top) and reference spectrum (below), acquired with recombinant HvExoI and Glc (40 µM and 2.4 mM, respectively) at 600 MHz and 283 K. Only protons of β-D-Glc (squares) showed STDs. Off and on resonance frequencies were set at 100 ppm and 0.65 ppm, respectively. The Gaussian shaped pulse (30 ms) was used for selective irradiation with total saturation time of 2 s. b trNOESY spectrum of recombinant HvExoI and Glc (57 µM and 285 µM), acquired at 800 MHz, 283 K and with 300 ms mixing time. Only NOEs observed were weak trNOEs of Glc between H1/H3 and H1/H5. c trNOESY spectrum of recombinant HvExoI with Glc and thiocellobiose (57 µM, 285 µM, and 171 µM, respectively), acquired at 800 MHz, 283 K and with 300 ms mixing time. All NOEs were negative. trNOEs defining the 4C1 conformation of non-reducing Glc of bound thiocellobiose are in boxes. Blue lines in (b) and (c) refer to residual noise signals 3dGlc (Fig. 2b; Supplementary Fig. 1b) and 4dGlc. The structures via C6-OH with Nε1 of Trp434 and Oε2 of Glu491, and via C4- of both complexes showed that dGlc derivatives could not dis- OH with the N atom of Gly57 (Fig. 2b; Supplementary Fig. 1b). place Glc from the −1 subsite, but that Glc consolidated to the Notably, as 3dGlc rotated to a different position relative to that of −1 subsite. The well-defined electron density map indicated that Glc in native HvExoI (such that the C4-OH group of 3dGlc 3dGlc was orientated parallel to Trp286 and Trp434 in the overlapped the position of C3-OH of Glc in the native structure), +1 subsite, through π-system stacking contacts, forming H-bonds it could establish H-bonds with Gly57 and Glc in the −1 subsite. 6 NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z ARTICLE a Glc bound in –1 subsite b Glc bound in +1 subsite c Isolated Glc 1 1 1 qy qy qy 3,O B O 3 3,O S2 S1 3,O B O B 3 O 3 S2 S1 S2 S1 B2,5 B1,4 B2,5 B2,5 B1,4 B1,4 5 4 1 S5 4 S1 4 1 C1 5 S1 C1 0 1 C1 5 S1 0 S5 0 S5 1,4 1,4 B 2,5 B 1,4 2,5 B 2,5 B B B 1 2 0 1 2 SO S3 SO 1 2 S3 B3,O S3 SO B3,O B3,O qx qx qx –1 –1 –1 –1 0 1 –1 0 1 –1 0 1 ΔG ΔG ΔG (kcal/mol) 0 4 8 12 16 20 24 28 32 (kcal/mol) 0 5 10 15 20 25 (kcal/mol) 0 4 8 12 16 20 24 28 32 Fig. 5 Conformational FEL maps of β-D-Glc bound in the active site of HvExoI. Conformational FEL of β-D-Glc in the −1 (a) and +1 (b) subsites are based on ab initio QM/MM metadynamics. c Conformational FEL map of isolated Glc34 (reproduced with permissions of the Journal of American Chemical Society; all rights reserved; https://pubs.acs.org/doi/10.1021/jacs.5b01156). Energy level lines are separated by 1 kcal/mol Alkyl-glucosides and PEG remove Glc from the active site. rigidity, and hence it projected with the reducing-end Glc moiety Contrary to what we observed with dGlc derivatives, n-octyl β-D- into the solvent at the putative +2 subsite, while the non- glucosides with tighter Ki values (Table 1) replaced Glc and reducing-end was held in the adjacent +1 subsite. For the bound across the active site (Fig. 2c; Supplementary Fig. 1c). −1 subsite, we assigned glycerol and PEG molecules that formed While aliphatic chains of alkyl β-D-glucosides were threaded H-bonds with the G2SG-OMe non-reducing end that also through the Trp286 and Trp434 aromatic clamp and did not established the H-bond with Gly57; the reducing moiety of interact with a surrounding environment, Glc moieties estab- G2SG-OMe interacted via C4-OH with Oε2 of Glu287 and via lished a network of 10–11 mono- and co-operative bi- and tri- C3-OH with Nε1 of Trp286 (Fig. 3d; Supplementary Fig. 2d). dentate H-bonds, like those observed in native (with trapped Glc) In the G6SG-OMe structure we noted significant elevations of B- or recombinant Glc-perfused HvExoI (Figs. 2, 3). Similarly, after factor values in two loops (Thr214-Glu228 and Glu491-Asn498) applying PEG at saturating concentrations to native HvExoI near the catalytic site entry (Supplementary Note 2); these crystals, the 1.80 Å PEG-perfused structure revealed that Glc was B-factor elevations were not seen in the G2SG-OMe structure. replaced by five water molecules occupying the positions of the C1-OH to C6-OH groups of Glc in the −1 subsite (Fig. 2d; HvExoI-G2SG-OMe complex is an intermediate for Glc exit Supplementary Fig. 1d). The PEG molecule bound between path. We considered the HvExoI in complex with G2SG-OMe Trp286 and Trp434 in two alternate conformations made a (PDB 6MD6) (Fig. 3d; Supplementary Fig. 2d) to embody the water-mediated H-bond with Gly57, which was reminiscent of attributes of an intermediate enzyme-substrate complex, and the the bond formed between C4-OH of 3dGlc and Gly57. The sec- disposition of G2SG-OMe that of an incoming substrate. Hence, ond PEG molecule bound near Trp434 did not contact the we used this complex for investigations of the hypothetical Glc enzyme, but instead interacted with the first PEG via water- displacement route via multi-scale molecular modelling employ- mediated H-bonds (Fig. 2d; Supplementary Fig. 1d). ing docking19 and MD simulations, followed by GaudiMM20 and PELE21 pathway calculations. G6SG-OMe and G2SG-OMe bind to HvExoI in two distinct Respective reciprocal docking of the Glc product and the β-D- modes. To understand how (1,6)- and (1,2)-linked substrate thio- glucopyranosyl-(1,2)-D-glucose (G2OG) or β-D-glucopyranosyl- analogues displace the Glc product from the active site, we per- (1,3)-D-glucose (G3OG) substrates into HvExoI:G2OG-OMe or fused recombinant crystals of HvExoI and observed that the HvExoI:Glc complexes, combined with MD simulations indicated sugars adopted different poses (Fig. 3c, d; Supplementary Fig. 2c, that the existence of ternary HvExoI:Glc:G2OG or HvExoI:Glc: d). It is worth noting that G6SG-OMe binds to the enzyme with G3OG complexes (designated 1–6) was plausible, where Glc was KD of 0.008 × 10–3 M, whereas G2SG-OMe binds about 320-fold bound in the −1 subsite, and G2OG or G3OG were attached at weaker (Ki of 2.55 × 10–3 M) (Table 1). We observed that in the +1 and putative +2 subsites (Fig. 6a; Supplementary Figs. 5, 6, 8, 1.57 Å structure of HvExoI with G6SG-OMe, the ligand formed Supplementary Notes 3–4, Supplementary Data 1). This well-defined electron densities for both Glc moieties across the suggested that if the incoming substrate binds at the +1 to active site (Fig. 3c; Supplementary Fig. 2c), and that the 1.68 Å +2 subsites, while Glc is still trapped in the −1 subsite, an structure with G2SG-OMe (refined at 0.7 occupancy at both alternative Glc exit path, other than through the +1 subsite needs subsites) had a less-defined density for the ligand, which bound to be considered. We also searched for substrate binding modes to the enzyme in a markedly different position (Fig. 3d; Supple- different to those at the +1 to +2 subsites, using docking of the mentary Fig. 2d). For G6SG-OMe, we observed short H-bonds in G2OG or G3OG substrates in structures based on MD the −1 subsite, which resembled those observed in the native simulations of the native HvExoI:Glc complex (Supplementary (Fig. 2a; trapped Glc) or recombinant (Fig. 3a; perfused with Glc) Figs. 5–8; Supplementary Note 4; Supplementary Data 1–2). enzymes. In contrast, G2SG-OMe could not slide past the However, no stable binding sites were found other than those +1 subsite in the pocket, likely due to a thio-glycosidic bond at +1 and +2 subsites. In addition, docking and MD simulations NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z a 3D model with Glc and G3OG b GaudiMM calculations Tyr253 VdW overlap (Å3) 240 Arg291 200 160 Asp285 120 Trp286 G3OG 80 Average = 69.82 40 Average = 31.74 0 0 2 4 6 8 10 12 14 Glc Glu220 Separations (Å) VdW overlap (Å3) –1 subsite +1 subsite +2 subsite 20 16 12 Average = 8.02 8 Glu491 Trp434 4 Average = 0.5 Glu161 0 0 2 4 6 8 10 12 14 Arg158 Lys493 Separations (Å) c Cork-like motion of domains 1 and 2 d PELE binding energies during Glc displacement 10 h e 0 f g –10 PELE energy (kcal/mol) –20 –30 –40 –50 Complex 3 Complex 2 –60 Complex 1 –70 –80 0 5 10 15 20 Separations (Å) e f Step 1 of the Glc displacement route Step 2 of the Glc displacement route Separations ~2–3 Å Separations ~5–7 Å Trp286 Trp286 Asp285 Asp285 Glc G3OG G3OG G3OG G3OG Glc Glc Glc Arg158 Trp434 Trp434 Glu491 Glu491 Arg158 g h Step 3 of the Glc displacement route Step 4 of the Glc displacement route Separations ~8–11 Å Separations ~13–17 Å Trp286 Asp285 Asp285 Trp286 Glc Glc Glc G3OG G3OG G3OG G3OG Trp434 Arg158 Glc Bulk solvent Glu491 Glu491 Arg158 unequivocally revealed that the HvExoI:Glc:G2OG and HvExoI: cavity adjacent to the −1 and +1 subsites. We observed that Glc:G3OG complexes (Fig. 6a) converged and were stable, when only one Glc was bound at the −1 subsite, Tyr253 remained and could be used for investigations of Glc displacement. MD H-bonded to the carbonyl oxygen of the Trp286 backbone for simulations also revealed that stable binding of substrates the full extent of simulation. The sidechains of the catalytic required Trp434 adopting the parallel orientation to that of nucleophile Asp285 and acid/base Glu491, together with Arg158 Trp286, where C-H/π interactions mediated binding. Notably, (Fig. 6a; dashed lines) physically separated this cavity from the binding of the G2OG or G3OG substrates triggered the −1 subsite. In complexes 1, 3 and 4, the Tyr253 side chain moved conformational change of Tyr253, modifying the buried lateral away from this lateral cavity enlarging it, while in complex 2, 8 NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z ARTICLE Fig. 6 Displacement pathway of Glc, based on crystal structures and molecular modelling. a Converged HvExoI:Glc:G3OG complex obtained by docking. Selected residues (carbons: atomic and purple sticks), Glc (carbons: yellow sticks) and G3OG (carbons: orange sticks) are shown. Glc in the −1 subsite is separated from the lateral cavity by the Arg158-Glu491-Asp285 (carbons: purple sticks) toll-like barrier (dashed lines). Surface morphology coloured by electrostatic potentials: white, neutral; blue, +5 kT e−1; red, −5 kT e−1. b Glc cannot exit, when crystallographic coordinates are considered by GaudiMM (grey symbols). Glc displacement only occurs when protein atoms are relocated to the lowest-energy normal mode (black). Clashes of Glc are shown as the function of separations from the Glc geometric centre in the −1 subsite. Displacement routes present fewer clashes, when MD structures are considered (green) and protein atoms move (dark red). c Lowest-energy normal mode shows a cork-like motion of domains 1 and 2. d Glc displacement routes calculated by PELE using ternary HvExoI:Glc:substrate complexes 1–3. PELE binding energies, plotted as a function of Glc separation from the −1 subsite, show that Glc traverses via the lateral cavity at the separation of 8 Å from protein surface. Full, semi-dashed, dashed and dotted lines specify separation ranges of 2–3 Å, 5–7 Å, 8–11 Å and 13–17 Å (referring to line thickness in (e–h)) between respective centres of masses and C4-OH groups of Glc in initial and final structures. e–h Four steps along the Glc displacement route based on complexes 1–3 (Glc carbons: yellow spheres). Left and right images are rotated by 200o (e–g) and −50o (h) along y-axes. e Glc is in the −1 subsite. f Movements of Arg158, Asp285 and Glu491 allow Glc traversing into the lateral cavity. g Glc in the lateral cavity is partly exposed to the bulk solvent. h Glc is fully exposed to the bulk solvent. Black dotted ellipsoid indicates the transient solvent-exposed aperture of the lateral cavity that facilitates Glc displacement this side chain rotated to the bottom of the cavity making it Ser252, Tyr253-backbone, Ser254, Asp285, Arg291, Glu491 and shallower (Supplementary Figs. 6, 8; Supplementary Notes 3–4). Thr492 (Fig. 6; Supplementary Fig. 9). This cavity was partially To find out if binding of substrates by HvExoI lacking the separated from the solvent by the Asn219, Glu220 and bound Glc product would lead to tighter or weaker binding Arg291 sidechains, with the two latter residues forming a salt than that with Glc included, we carried out further calculations bridge with Oε1 and Oε2 of Glu220 to NH1 of Arg291 at and compared binding affinities. Docking calculations of G2OG, separations of 3.10 Å/2.65 Å, 2.77 Å/3.56 Å and 4.65 Å/2.84 Å for G3OG and β-D-glucopyranosyl-(1,6)-D-glucose (G6OG) disac- complexes 1, 2, and 3, respectively. The abundance of data points charide substrates at the −1 and +1 subsites, predicted binding in the PELE plot at separations between 8–11 Å (Fig. 6d) with higher affinities (Goldscore scoring function values of 66 for indicated that there was a local energy minimum for Glc binding G2OG, 76 for G3OG and 74 for G6OG), when the Glc product in this cavity (Fig. 6d; subpanel g). For complex 2, PELE was absent in the −1 subsite. However, when the Glc product was calculations suggested a more transient passage of Glc through included in the −1 subsite, docking of G2OG, G3OG and G6OG the lateral cavity, most likely due to a shallower profile resulting in the +1 and putative +2 subsites predicted binding with lower from the specific Tyr253 conformation, compared to other affinities (Goldscore scoring function values of 60 for G2OG, complexes. Nevertheless, in all cases, if Glc was to advance to the 57 for G3OG and 61 for G6OG). These data indicated that bound lateral cavity from the −1 subsite, it must traverse the space Glc lowered binding energies for incoming substrates as they between Asp285, Glu491 and Arg158, which likely represent a had no access to the higher affinity −1 subsite. toll-like barrier (Fig. 6a, dashed lines; Fig. 6e); this Glc passage corresponded to fewer data points at separations between 5–7 Å Discovery of substrate-product assisted processivity in HvExoI. (Fig. 6d; Supplementary Data 4). With Glc bound in the Initially, we used the GPathFinder extension of the GaudiMM −1 subsite, the shortest separations between the Arg158- platform20 to reveal potential Glc exit routes from the active site Asp285, Arg158-Glu491 and Glu491-Asp285 sidechains were into the bulk solvent, whereby we considered steric clashes. The those at 5.10 Å/5.42 Å/5.17, 3.09 Å/3.16 Å/4.42 Å and 5.39 Å/6.32 initial set of simulations was performed with the HvExoI:Glc: Å/6.65 Å in respective complexes 1, 2 and 3. Notably, conforma- G2OG ternary complex 1, based on the coordinates of HvExoI tional changes of Arg158, Asp285 and Glu491 sidechains in complex with G2SG-OMe and Glc docked in the −1 subsite. facilitated Glc movement (Fig. 6e, f), and altered the separations No exit path was detected in this case, as Glc never overcame specified above to 7.07 Å/5.33 Å/6.21 Å, 3.00 Å/2.97 Å/6.82 Å and steric clashes to vacate the −1 subsite (Fig. 6b; top panel; Sup- 6.96 Å/7.23 Å/7.71 Å, when Glc bypassed the toll-like barrier. plementary Data 3). However, further calculations with the pro- Once in the lateral cavity, Glc was free to exit from the cavity into tein backbone displaced by 2.2 Å along the lowest-energy normal the bulk solvent via a transient and autonomous aperture. This mode strikingly identified the Glc displacement path. As a bonus, transient opening was formed through rotations and backbone these calculations also revealed a cork-like motion between fluctuations of Glu220, Arg291, Thr492 and Lys493, and domains 1 and 2, which was correlated with the open and closed surrounding residues in the vicinity of the bound β-D-glucoside states of the active site (Fig. 6c). The same ligand path was also molecule (Fig. 6h). As Glc migrated across the lateral cavity obtained with the MD-derived HvExoI:Glc:G3OG ternary com- (Fig. 6e, f), it established H-bonds36 with the protein residues and plex 3 (Fig. 6b; bottom panel), highlighting the importance of incoming substrates to maintain energetic favorability (Supple- incorporating local and collective protein motions in these mentary Tables 2–4). This suggested that the hydrophilic calculations. environment of the toll-like barrier and the lateral cavity may To explore the Glc displacement route in detail, the ligand have evolved for this exit route and may be evolutionarily migration calculations of the PELE approach21 were used, conserved in GH3 enzymes (Supplementary Fig. 9). Finally, and whereby structures generated at each step along the path (by most importantly, the Glc displacement route raised the sequential ligand and protein geometric perturbations followed possibility that this trajectory may facilitate processive catalysis by energy minimisation using the OPLS-AA force field) were in HvExoI and other GH3 exo-hydrolases with pocket-shaped evaluated and accepted or rejected according to a Metropolis active sites. criterion at a given temperature. Using this approach, the HvExoI:Glc:G2OG and HvExoI:Glc:G3OG complexes 1–3 con- Proof of concept using the R158A/E161A variant. Ligand verged, and revealed that Glc egressed from the −1 subsite to the migration calculations using PELE (Fig. 6) and the conservation adjacent lateral cavity that was enlarged after Tyr253 rotated. This patterns of the toll-like barrier (Supplementary Fig. 9) led us cavity was defined by an ensemble of 14 residues: Trp156, designing the double non-conservative R158A/E161A variant to Arg158, His207, Phe208-backbone, Asp211, Asn219, Glu220, critically asses the roles of these residues in binding Glc and the NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z a Variant R158A/E161A in ligand-free form observation in native HvExoI isolated from seedlings6, which renders it an archetype model to examine product/substrate Asp285 Trp286 Asp285 Trp286 dissociations/associations in a pocket-shaped active site. Most His207 Glu491 (A, B) His207 Glu491 (A, B) approaches towards the descriptions of reactant movements are based on MD simulations of substrate and product binding R158A Gol 2 R158A and unbinding, including random collisions or diffusion40. Lys206 Gol 1 Trp434 Lys206 Trp434 Other studies reveal substrate/product migration pathways that Gly57 Gly57 accompany conformational changes of protein backbones or E161A Asp95 E161A Asp95 secondary structures41,42. –1 subsite +1 subsite When we first described the HvExoI structure6–9, we were b Variant R158A/E161A perfused with Glc unable to explain why the trapped Glc product remained bound in the active site and has not diffused away, and what was the Asp285 Trp286 Asp285 Trp286 implication of this unassuming observation. The obvious expla- Gol 2 nation for not seeing naturally bound products or co-factors His207 Glu491 (A, B) His207 Glu491 (A, B) in structures of other GH enzymes is that these proteins are R158A R158A generated in recombinant hosts, where intracellular concentra- Lys206 Gol 1 Trp434 Lys206 Trp434 tions of potential enzyme reactants are not high enough during Gly57 Gly57 protein maturation. To this end, solving crystal structures E161A Asp95 E161A Asp95 of enzymes purified from native sources offers an additional –1 subsite +1 subsite information that could be beneficial for the understanding of catalytic cycles. c Variant R158A/E161A perfused with G6SG-OMe When we first detected the Glc product entrapment in native G6sG-OMe HvExoI, we hypothesised that the role of Glc may be linked to Trp286 Trp286 Asp285 Asp285 the pre-organised state of the active site to maintain efficient His207 Glu491 His207 Glu491 catalysis required for the growth of a plant embryo9. One way R158A R158A of preserving this pre-organised state would be to keep entropic Lys206 Lys206 costs of catalysed reactions low, where an entrapped product Trp434 Trp434 Gly57 Gly57 from a previous hydrolytic cycle could lower overall entropic E161A Asp95 E161A Asp95 demands for binding of incoming substrates1,43. The Glc product –1 subsite +1 subsite retention in the active site may also disfavour binding of incorrect Fig. 7 R158A/E161A of HvExoI and with perfused Glc or G6SG-OMe. substrates through product dissociation rates that would govern a Stereo view of a ligand-free form of R158A/E161A. Two glycerol the selectivity of substrate binding44. Others may argue that molecules (carbons: green sticks) at 1.0 occupancies are bound in the −1 although most enzymes have evolved to use conformational and +1 subsites. b Stereo view of R158A/E161A with perfused Glc, which adjustments to favour tight binding of correct substrates and failed to bind, instead two glycerol molecules (carbons: green sticks) are product release, some enzymes may have not acquired this asset. bound at 1.0 occupancies in the −1 and +1 subsites. c Stereo view of This assumption was confirmed by perfusing dGlc derivatives R158A/E161A with G6SG-OMe (carbons: orange sticks) at 1.0 occupancy, that akin of reaction products could not remove Glc, but substrate bound across the −1 and +1 subsites. Separations of less than 3.30 Å from analogues and mimics could. electronegative atoms of active site residues (carbons: grey) and water As mentioned above, despite the plethora of structural data for oxygens, are shown as dashed lines. Derived |2mFobs − DFcalc| electron HvExoI6–9,15–17, we could not explain, how Glc egresses from the density maps are contoured at 1σ (blue mesh) 13 Å-deep pocket-shaped active site, although we detected cross- talk between residues entrapping Glc in the −1 subsite, and those in the +1 subsite delineated by Trp286 and Trp434. Function of G6SG-OMe substrate thio-analogue (Fig. 7; Supplementary Trp residues in retaining Glc or substrates is unsurprising, as Fig. 10). The 1.65 Å and 2.21 Å structures without and with these residues have extended heterocyclic indole ring systems perfused Glc showed two glycerol molecules, one each in the −1 with amphipathic characteristics, high de-localised electron den- and +1 subsites (Fig. 7a, b; Supplementary Fig. 10a, b), and the sities and permanent dipole moments due the intra-ring nitro- absence of Glc despite perfusing crystals at saturating con- gens, that form H-bonds29. In all 22 structures of HvExoI we so centrations. On the other hand, the 2.30 Å structure of R158A/ far resolved (Supplementary Table 5), Trp286 and Trp434 at the E161A perfused with G6SG-OMe strikingly revealed this sugar +1 subsite showed no structural heterogeneity, and at separations bound across the active site, but in a different pose compared to of about 4 Å from ligands were compliant with the provision of that of the wild-type (WT) enzyme (cf. Figures 3c, 7c; Supple- two sets of C-H/π interactions45 depending on the stereo- mentary Figs. 2c and 10c). In the variant structure the position chemistry of bound ligands7,17. Crystalline HvExoI was also of the Glc moiety in the −1 subsite matched that of WT, but perfused with the (3)-β-D-Glc-S-(1,3)-β-D-Glc-(1,)14–22 polymeric the disposition of the reducing-end Glc moiety was flipped in its substrate46, but we did not observe potential binding sites beyond position between the indole moieties of Trp286 and Trp434, such the −1 and +1 subsites. that intra-ring oxygen of the β-D-glucopyranose moiety pointed The time-space averaged vision of crystal structures presented to Trp434. Obviously, this G6SG-OMe pose was due to Ala here was supported by the multi-scale molecular modelling, replacements for Arg158 and Glu161. based on docking19, MD simulations, QM/MM, GaudiMM20 and PELE21, that provided the first clear view of succession of events during Glc displacement. This approach also hinted that trapped Discussion Glc may play a surprising role in catalysis and led to the Substrates and products in enzyme active sites bind and unbind at discovery of processivity by this exo-hydrolase. In the light of this fast rates that challenge investigations of their trajectories, leading discovery we coin the term ‘substrate-product assisted processive to the lack of a deep understanding of these hallmarks of enzyme catalysis’, due to the key role that a substrate and a product play catalysis37–39. Here, we capitalise on the Glc product entrapment in the evocation of the catalytic pathway. 10 NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z ARTICLE a Trp286 Trp286 Trp286 Trp286 Trp434 Trp434 Trp434 Trp434 –1 subsite +1 subsite –1 subsite +1 subsite –1 subsite +1 subsite –1 subsite +1 subsite Step 1: Step 2: Step 3: Step 4: Dimeric substrate binds in Dimeric substrate is hydrolysed, After new dimeric Dimeric substrate progresses productive mode aglycon diffuses away, substrate binds, to catalytic site. Glc exits in –1 and +1 subsites. Glc product is bound non-covalently Glc product is via lateral cavity formed and moves between consollidated transiently near +1 subsite. –1 and +1 subsites. to –1 subsite. New hydrolytic event can begin. b c Trp286 Trp286 Trp286 Trp286 Product Product Trp434 Trp434 Trp434 Substrate –1 subsite +1 subsite –1 subsite +1 subsite –1 subsite +1 subsite Substrate Step 1: Step 2: Step 3: Trp434 Polysaccharide Non-reducing Glc residue Glc product modifies its substrate binds is cleaved off and binding patterns and in productive stays bound in –1 subsite. exits via lateral cavity mode in The polymeric substrate, formed transiently –1, +1 and +2 short of Glc residue, near +1 subsite. subsites. is attached in New hydrolytic event +1 and +2 subsites. can begin. Fig. 8 The mechanism of substrate-product assisted processive catalysis by HvExoI. a Mechanism of Glc displacement with a disaccharide. After the disaccharide (empty blue and filled grey squares) bound in −1 and +1 subsites (step 1) is hydrolysed and an aglycon diffuses away (step 2), Glc (cyan square) remains non-covalently trapped and oscillates (double arrow) between the −1 and +1 subsites. Glc (yellow square) is consolidated to the −1 subsite after an incoming substrate binds (step 3), which later advances to the catalytic site. Glc (red square) modifies its binding patterns and exits (large arrow) via an autonomous and transient lateral cavity (cylinder in dashed lines) formed near the catalytic site (step 4); a next hydrolytic cycle can begin. b Mechanism of Glc displacement with a polysaccharide to facilitate substrate-product assisted processive catalysis. After non-reducing (blue square) and penultimate (filled square) residues bind in a productive mode at the −1 and +1 subsites (step 1), the non-reducing Glc residue (yellow square) is cleaved off, with the remainder of the substrate attached (step 2). Glc (red square) modifies its binding patterns and is released (large arrow) via a lateral cavity (cylinder in dashed lines, step 3). Here, the next hydrolytic cycle continues with the same polysaccharide, where the polysaccharide, short of one Glc (filled squares) advances into the catalytic site after uninterrupted binding. Non-reducing (blue) and hydrolysed (cyan-yellow-red in (a), yellow-red in (b)) Glc is shown. c Structural basis of substrate-product assisted processive catalysis. Left: The substrate (orange spheres) slides into the active site (protein surface: domain 1-light grey, domain 2-dark grey). After hydrolysis, the Glc product (yellow spheres) exits from the −1 subsite through the lateral cavity into the bulk solvent. Right: detail of the active site. Circle line with arrows indicates directionality during processive catalysis. The image was generated using the coordinates of HvExoI:G6SG-OMe and HvExoI:G2SG-OMe complexes (Fig. 3c–d; orange spheres), and the Glc displacement snapshots based on PELE calculations (Fig. 6e–h; yellow spheres) We demonstrate that the substrate binding and product dis- (yellow square) is consolidated to the −1 subsite. This is linked to placement routes proceed through stages that are carefully Tyr253 rotation, which enlarges the lateral cavity adjacent to the orchestrated in succession (Fig. 8). We suggest that a disaccharide −1 and +1 subsites (step 3). Next, the Glc product (red square) (Fig. 8a) binds near the catalytic site at the +1 and putative traverses from the −1 subsite through rotations of the Arg158 +2 subsites most likely through random collisions, from where and Asp285 sidechains and associated backbone atoms, and it progresses to the −1 and +1 subsites. Following binding in bypasses the toll-like barrier shaped by Asp285, Glu491 and a productive mode (step 1), the disaccharide (empty blue and Arg158, via this autonomous lateral cavity (step 4, cavity shown filled grey squares) is poised for hydrolysis. After the non- as a barrel in dashed lines). This cavity is formed transiently near reducing Glc is cleaved off and aglycon diffuses away, the Glc the +1 subsite, through which the Glc product exits (arrow) via product is entrapped non-covalently, whereby it may oscillate an aperture into the bulk solvent (Fig. 8a; step 4). Strikingly, when between the −1 and +1 subsites (cyan square) (step 2). After the the toll-like barrier is disrupted in the R158A/E161A variant, Glc next disaccharide binds at +1 and +2 subsites, the Glc product no longer binds in the active site (Fig. 7), which underlines the NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z vital role of the barrier in Glc entrapment. The sequence of events The discovery of substrate-product assisted processive catalysis involving incoming substrate binding and Glc displacement are in HvExoI prompts investigations of the evolutionary origin of captured in Supplementary Movie 1. this mechanism. The availability of the presented experimental We propose a similar substrate/product progression route and computational tools alongside a rich source of information (Fig. 8b) takes place with a polysaccharide (empty blue and fil- for the GH3 family enzymes that originate from all phylla9,18 led grey squares), often with differently linked β-D-glucosyl (currently around 23,000 entries), will now allow the study of residues8,19, through processive catalysis. Here, (step 1) after Glc evolution of this catalytic mechanism. Preliminary analyses (blue square) is cleaved off from the polysaccharide, (step 2) the of 500 sequences related to HvExoI revealed that the residues Glc product (yellow square) remains bound at the −1 subsite, underlying this mechanism are conserved in land plants, but while the remainder of the polysaccharide (grey squares) stay absent in red and green algae, suggesting that this mechanism is attached at the +1 and putative +2 subsites. As the Glc product rather ancient and has evolved in land plants about 470 million (red square) exits (arrow) through the lateral cavity (step 3, years ago. shown as a barrel in dashed lines), the polysaccharide can We conclude that through the descriptions of associative and advance to the −1 and +1 subsites, such that an uninterrupted dissociative reactant trajectories that we explored using our inter- hydrolysis of the polymer continues until it is completely disciplinary approach, we are now in a better position to hydrolysed. In both instances, after Glc egresses (Fig. 8a, b), the understand how reactants advance in active sites, and how to next hydrolytic cycle that is facilitated by the dedicated substrate improve enzyme catalytic rates, stability, product inhibition and binding and product displacement routes, can continue. Conse- drug discovery. quently, this coordinated system with oligo- and polymeric sub- strates, resembles a loop or a conveyer belt that directionally shuffles substrates and ensures efficient catalysis (Fig. 8c). Methods Materials. 3-Deoxy-glucose (3dGlc) and 4-deoxy-glucose (4dGlc) were obtained Although in all computations, we modelled Glc in 4C1 con- from Glycoteam GmbH (Hamburg, Germany). N-octyl β-D-glucopyranoside formation, we do not rule out that the Glc distortion to a higher (octyl-O-Glc) and n-octyl 1-thio-β-D-glucopyranoside (octyl-S-Glc) were from energy of the 1S3 conformer, as it occurs in crystal structures, may Anatrace Inc. (Maume, OH, USA). 4-Nitrophenyl β-D-glucopyranoside (4NPGlc), glycerol, myo-inositol and other chemicals were obtained from Sigma-Aldrich contribute to displacement of the Glc product or its progression. (USA) or as described previously6,7,17. CH3CN was from BDH Laboratory Supplies To our knowledge, this type of processive catalysis47–49 has not (UK), and polyethylene glycol 400 (PEG400; n = 5–10) was from Fluka Biochemica been described in any exo-acting hydrolase with a pocket-shaped (Germany). active site like that of HvExoI. In this work, we define processive catalysis in broader terms, rather than those based on hydrolytic Synthesis of methyl 6-thio-β-gentiobioside. The synthesis proceeded in two profiles of products50. Considering that HvExoI operates on steps as follows. polysaccharides such as (1,3;1,4)-β-D-glucans contained in plant First, methyl 6-thio-β-gentiobioside heptaacetate (chemical structure in cell walls9,51, this finding is unsurprising, as the enzyme’s efficient Supplementary Table 1): 2,3,4,6-Tetra-O-acetyl-1-thio-β-D-glucopyranose (845 mg, hydrolysis is the prerequisite for seed germination and the 2.32 mmol)56, 1,4-dithioerythritol (357 mg, 2.32 mmol) and cysteamine (1179 mg, 2.32 mmol) were successively added to a solution of methyl 2,3,4-tri-O-acetyl-6- development of a plant embryo9,23,51. Processive catalysis (non- deoxy-6-iodo-β-D-glucopyranoside (1 g, 2.32 mmol) in hexamethylphosphoramide dissociative sequential degradation) has been accepted for endo- (HMPA) (7 mL)57. The mixture was kept for 1 h at room temperature then acting hydrolases with cleft- or groove-shaped catalytic sites, such precipited into ice-water (100 mL). The solid was collected on Celite, washed with as (1,3)- and (1,3–1,4)-β-D-glucan endohydrolases52, and (1,4)-β- water and then dissolved in CH2Cl2. The organic phase was washed with water, dried over Na2SO4 and concentrated. Column chromatography in silica gel D-glucan endohydrolases53, and for exo-acting chitinases50 and (EtOAc-light petroleum 1:1) afforded methyl 6-thio-β-gentiobioside heptaacetate cellobiohydrolases54 with ridge- or tunnel-shaped sites. Con- (1.3 g, 84% yield). [α]25 D -20.0 (c 0.73, CHCl3); ESI-HRMS(+) calcd for versely, β-D-glucosidases and exo-acting hydrolases with funnel-, C27H38O17NaS 689.17219; found 689.17081; 1H NMR (CDCl3): δ 5.16 (t, 1 H, J = crater- or pocket-shaped6,55 active sites are deemed to be non- 9.5 Hz), 5.15 (t, 1 H, J = 9.5 Hz), 5.0 (t, 1 H, J = 9.8 Hz), 4.95 (t, 1 H, J = 10.0 Hz), processive48,49. 4.93 (dd, 1 H, J = 8.0 Hz, J = 10 Hz), 4.88 (t, 1 H, J = 9.5 Hz), 4.65 (d, 1 H, J = 10.0 Hz), 4.39 (d, 1 H, J = 8.0 Hz), 4.22 (dd, 1 H, J = 5.0 Hz, J = 12.5 Hz), 4.12 Our concept of substrate and product processive co-operative (dd, 1 H, J = 2.0 Hz, J = 12.5 Hz), 3.65 (m, 2 H), 3.51 (s, 3 H), 2.79 (m, 2 H), routes with pocket-shaped exo-hydrolases6–9 is associated with 2.08–1.97 (m, 21 H); 13C NMR (CDCl3): δ 169.6–169.3, 101.6, 83.9, 76.0, 74.4, 73.8, chemical signalling, where entrapped products, bound sub- 72.7, 71.9, 71.3, 70.2, 68.2, 62.1, 57.1, 31.2, 20.77–20.72. strates and enzymes co-operatively create pathways dedicated Second, methyl 6-S-β-D-glucopyranosyl-6-thio-β-D-glucopyranoside (methyl 6- thio-β-gentiobioside; G6SG-OMe) (chemical structure in Supplementary Table 1): to product dissociation. The hallmark and advantage of pro- sodium methylate 1 M (300 µL) was added to a solution of methyl 6-thio-β- cessive catalysis is a high efficiency, where oligo- or poly- gentiobioside heptaacetate (300 mg, 0.45 mmol) in MeOH (30 ml). After 4 h at saccharides are directionally threaded without the catalyst losing room temperature, the mixture was neutralised with Amberlite IR 120 H+ resin, contact with substrates to keep binding entropic cost low, filtered and concentrated. The residue was dissolved with water, and freeze dried affording methyl 6-thio-β-gentiobioside (166 mg, 99% yield). [α]25 D -26.2 (c 1.13, while performing multiple hydrolytic events on the same sub- H2O); ESI-HRMS( + ) calcd for C13H24O10NaS 395.09824; found 395.09750; 1H strate molecule. While this work was focussed on a plant GH3 NMR (D2O): δ 4.69 (d, 1 H, J = 10.0 Hz), 4.43 (d, 1 H, J = 8.0 Hz), 3.94 (dd, 1 H, hydrolase, we suggest that substrate-product assisted processive J = 2 Hz, J = 12.5 Hz), 3.75 (dd, 1 H, J = 5.5 Hz, J = 12.0 Hz), 3.68 (m, 1 H), 3.62 catalysis may be more prevalent among exo-hydrolases with (s, 3 H), 3.54–3.24 (m, 8 H), 2.97 (dd, 1 H, J = 8.0 Hz, J = 14.5 Hz); 13C NMR pocket-shaped active sites, irrespective of their substrate speci- (D2O): δ 103.3, 86.0, 79.8, 77.2, 76.0, 75.5, 73.2, 72.5, 72.4, 69.5, 60.9, 57.4, 31.3. ficity. To this end, our findings could promote further investi- gations of enzymes involved in biomass degradation, where one Synthesis of methyl 2-thio-β-sophoroside. The synthesis proceeded in three of the challenges is to circumvent the bottleneck of product steps as follows. inhibition, a problem that could be addressed through the First, methyl 3,4,6-tri-O-acetyl-2-O-trifluoromethanesulfonyl-β-D- mannopyranoside (chemical structure in Supplementary Table 1): Distilled better understanding of enzyme catalytic cycles. More broadly, pyridine (0.5 mL) and trifluoromethanesulfonic anhydride (113 µL, 3.6 equiv) substrate-product assisted processive catalysis in enzymes with were successively added to an ice-cold solution of methyl 3,4,6-tri-O-acetyl-β- pocket-shaped catalytic sites is likely to have significance D-mannopyranoside58 (60 mg, 0.19 mmol) in anhydrous CH2Cl2 (3 mL). The in other enzymes due to the plasticity of protein structures. It mixture was stirred at 0 °C for 30 min and at room temperature for 1 h. The solution was diluted with water and extracted with CH2Cl2. The organic phase was is exactly this plasticity that accounts in HvExoI for product washed with aqueous KHSO4 (20%, v/v), saturated aqueous NaHCO3, dried over dissociation, although in other enzymes the atomic details of anhydrous Na2SO4 and concentrated. The crude triflate (84 mg) was obtained in a reactants movements may be different. quantitative yield and directly used without purification. 12 NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z ARTICLE Second, methyl 3,4,6-tri-O-acetyl-2-S-(2,3,4,6-tetra-O-acetyl-β-D- cation-exchange and immobilised metal affinity chromatography[24,25. Purified glucopyranosyl)-2-thio-β-D-glucopyranoside (methyl 2-thio-β-sophoroside recombinant wild-type (WT) and variant HvExoI were crystallised via a hanging- heptaacetate) (chemical structure in Supplementary Table 1): diethylamine (1 mL) drop vapour-diffusion method, using macro- and cross-seeding with WT native under Argon was added to a stirred solution of 2,3,4,6-tetra-O-acetyl-1-thio-β- crystal seeds. The fully-grown crystals of recombinant HvExoI reached up to D-glucopyranose59 (91 mg, 0.22 mmol), to which crude methyl 3,4,6-tri-O-acetyl-2- 500 × 250 × 375 μm sizes after 5 to 14 days and were used for diffraction25. O-trifluoromethanesulfonyl-β-D-mannopyranoside (84 mg, 0.19 mmol) in anhydrous DMF (10 mL) was added. The mixture was stirred overnight at room temperature and then concentrated. A solution of the residue in CH2Cl2 was Crystal structure determination. Enzyme crystals with lengths of up to 400 μm washed in water, dried over anhydrous Na2SO4 and concentrated. Flash in the longest dimensions were transferred in 100 mM HEPES-NaOH buffer, chromatography on silica gel (ethyl acetate/petroleum ether, 1:1) afforded methyl pH 7.0, containing 1.2% (w/v) PEG (n = 5–10) and 1.7 M ammonium 2-thio-β-sophoroside heptaacetate (26 mg, 21%). Analytical data agreed with those sulphate (solution A), where 3-deoxy-glucose (3dGlc), 4-deoxy-glucose (4dGlc), reported previously60. ESI-MS: m/z = 689 [M+Na]+; 1H NMR (CDCl3) δ 5.12 n-octyl β-D-glucopyranoside (octyl-O-Glc), n-octyl 1-thio-β-D-glucopyranoside (t, 1 H, J = 9.3 Hz), 5.04 (t, 1 H, J = 9.4 Hz), 4.97 (m, 3 H), 4.73 (d, 1 H, J = (octyl-S-Glc), methyl 2-thio-β-sophoroside (G2SG-OMe), methyl 6-thio-β- 10.2 Hz), 4.36 (d, 1 H, J = 8.5 Hz), 4.30 (dd, 1 H, J = 4.6, J = 12.3 Hz), 4.26 (dd, 1 H, gentiobioside (G6SG-OMe) or Glc were supplied during ligand perfusions. The J = 4.8, J = 12.5 Hz), 4.13 (m, 2 H), 3.66 (2ddd, 2 H), 3.56 (s, 3 H), 3.06 (dd, 1 H, final concentrations of ligands were between 5–20 mM and perfusion proceeded J = 8.7 Hz, J = 10.7 Hz), 2.08–1.99 (7 s, 21 H); 13C NMR (CDCl3) δ 170.57–169.22 during 5–720 min. Alternatively, the crystals were deposited in a solution A with (CO), 103.79, 83.79, 76.03, 74.02, 73.11, 71.77, 71.48, 69.51, 68.25, 62.25, 62.15, the concentration of PEG increased to the 35% (v/v). After various perfusion 57.40, 49.90, 20.86–20.77 (CH3). times at 4 ± 2 °C crystals were cryo-protected with 15% (v/v) glycerol in the Third, methyl 2-S-β-D-glucopyranosyl-2-thio-β-D-glucopyranoside (methyl 2- solution A and mounted on a synchrotron goniometer in a stream of N2 gas thio-β-sophoroside; G2SG-OMe) (chemical structure in Supplementary Table 1): at 100 K (Oxford Instruments, UK). X-ray diffraction data (native, and with 1 M sodium methylate (200 µL) was added to a solution of methyl 2-thio-β- 3dGlc and 4dGlc) were collected on the undulator beamline BioCARS 14-ID-B sophoroside heptaacetate (26 mg, 0.04 mmol) in MeOH (5 mL) and stirred at at the University of Chicago Advanced Photon Source (USA) equipped with a ambient temperature for 2 h. The mixture was neutralised with the Amberlite IR bent cylindrical Si-mirror (Rh coating) Diamond(111) double-bounce mono- 120 H+ resin, filtered and concentrated. After lyophilisation, pure methyl 2-thio-β- chromator and focussing to the MARCCD-165 detector. X-ray diffraction data sophoroside was isolated with 97% yield (14.5 mg). [α]D −15.02 (c 0.3, H2O), ESI- sets of crystals with PEG and with octyl-O-Glc and octyl-S-Glc were collected MS: m/z = 395 [M+Na]+; 1H NMR (CDCl3) δ 4.76 (d, 1 H, J = 10.0 Hz), 4.57 (d, 1 using the multipole wiggler beamline BL5 at the Photon Factory (Japan), which H, J = 9.0 Hz), 3.96 (m, 2 H), 3.78 (m, 2 H), 3.66 (dd, 1 H, J = 8.8 Hz, J = 10.8 Hz), was fitted with a collimating mirror, double-crystal Si(111) monochromator and 3.61 (s, 3 H), 3.56–3.33 (m, 6 H), 2.85 (dd, 1 H, J = 9.2 Hz, J = 10.6 Hz); 13C NMR a focusing to the ADSC Quantum 315 CCD detector. X-ray diffraction data from (CDCl3) δ 102.36, 84.43, 79.76, 77.02, 75.61, 74.97, 72.82, 70.63, 69.32, 60.77 (2 C), recombinant WT and variants were collected at the MX1 and MX2 beamlines 57.23, 51.26. of the Australian Synchrotron (Australia) at 100 K (Oxford Instruments) or at ambient temperature (291 K) with a collimating mirror, double-crystal Si(111) monochromator 03BM1 dipole/bending magnet and ADSC Quantum 210r Steady-state affinity binding with recombinant HvExoI. Interaction between Detector. In the latter case, crystals were mounted using transparent polymer recombinant HvExoI and analytes Glc (molecular mass 180.2 g), G6SG-OMe capillaries from the RT Screening kit (MiTeGen). All data were collected at (372.4 g) and PEG (n = 5–10); (380–420 g) were performed using Surface Plasmon 0.5–1° oscillations throughout 180–720°. Data were processed using the Resonance (SPR) at 25 °C using the Biacore X100 V2.0.1 instrument (GE DENZO/SCALEPACK HKL 2000 suite of programmes61. Autoindexing deter- Healthcare, USA) equipped with the plus package. Purified recombinant HvExoI mined that space group P43212 of crystals were consistent with the tetragonal at 99% homogeneity (detected by SDS-PAGE) was covalently attached to a space group P43212 in all instances. The structures of the enzyme-ligand CM5 chip (GE Healthcare) in 10 mM sodium acetate buffer, pH 5.0. The car- complexes were refined using CCP4 REFMAC562 and PHENIX63. We used the boxymethyl groups present on the dextran layer of the CM5 chip were activated solved structures (PDB accessions: 1EX1, 1IEQ, 1IEV, 1IEW, 1IEX and 1J8V) with 1-ethyl-3-(3-dimethylaminopropyl) carbodiimide hydrochloride (EDC) and without ligands, ions, glycerol and water molecules, as starting models for N-hydroxysuccinimide (NHS) (GE Healthcare). Any unoccupied reactive sites refinements of HvExoI in complex with ligands7,15–17. The iterative model were blocked by a 7 min injection of 1 M ethanolamine hydrochloride-NaOH, pH building using XtalView/MIFIT64, Coot65 and refinements with REFMAC562 8.5. A ligand density of ~15.000 RU of recombinant HvExoI was achieved. Due to allowed tracing of all the residues. Following convergence in a standard refine- the fast rate constants of interactions, it was impossible to obtain individual rate ment, a further improvement of about 2% in the Rwork/Rfree factors ratio was constants ka and kd. Hence equilibrium analyses were performed for all protein- achieved by refining domain 1, (residues 1–357) and domain 2, (residues analyte interactions. For each interaction serial dilutions of analytes covering the 374–559), as two independent anisotropic entities with translation-libration- 10 · KD − 0.1 · KD concentration ranges were prepared in 10 mM HEPES-NaOH screw (TLS) motion66,67. Electron densities for the ligands were well-defined in buffer, pH 7.4 with 150 mM NaCl. Samples were injected onto the surface of a the active site regions at 3σ level in m|Fo|-|Fc| maps and water molecules were chip at a rate of 30 µl/min for 90 s, after which the formed enzyme/analyte located automatically with CCP4 ARP at levels higher than 3σ. Water molecules complexes could dissociate for 90 s. As all analytes dissociated completely from the were retained if they satisfied H-bond criteria and if their 2 m|Fo|-D|Fc| electron surface of the recombinant enzyme immobilised on a chip, no regeneration steps density maps were confirmed after refinements, where m is the figure of merit were needed. All responses were double referenced using a blank-immobilised and D is an estimated coordinate error. During model building and refinement, reference surface (NHS + EDC-activated/deactivated) and a series of blank injec- 5% of the data were flagged for cross-validation to monitor the progress of tions conducted to account for the system and injection artefacts. For each analyte refinements, using Rfree statistics68,69. The PROCHECK programme69 was used at least one concentration was injected in duplicate. Triplicate experiments to check the geometrical quality of the models. Ramachandran plots70,71 of all were performed, and the double-referenced data were globally fitted using a 1:1 structures showed that 99.8% of residues were found in the most favourable, equilibrium interaction model and the Biacore T200 v2.0 analysis software additionally allowed and generously allowed regions of the plot, with well- (GE Healthcare). Steady-state affinity binding constants KD (Supplementary Fig. 3) defined electron density for Ile432, the only residue located in disallowed were calculated 30 s after the onset of injections using a 5 s window in an average regions. Data collection, cell parameters and refinement statistics of all structures response in the fit. are summarised in Supplementary Tables 7–10. Graphics images were prepared with PyMol (Schrödinger, USA); electrostatic potentials were generated with a solvent probe radius at 1.40 Å, and protein and solvent dielectric constants Inhibition constants Ki with native HvExoI. Ki values of deoxy-glucose and alkyl 2 and 78, respectively. β-D-glucoside derivatives, and G2SG-OMe with native HvExoI were derived from hydrolytic rates. Enzyme inactivation was monitored at 30 °C by incubating 8 μmoles of the enzyme in 100 mM sodium acetate buffer, pH 5.25, containing GC/MS with evaporative light scattering detection. Six crystals of native 0.8–6.6 mM 4NPGlc as a substrate, 160 μg/ml BSA, and 3dGlc, 4dGlc, octyl-O-Glc, HvExoI with well-defined edges in a longest axis dimension of ~80–120 μm were octyl-S-Glc and G2SG-OMe. Each inhibitor was tested at six concentrations at collected from a hanging drop in a cat whisker, deposited in a test tube, washed in 0.4–3 times the Ki values in duplicate. We used the v = Vmax S / [KM × (1 + I/Ki) + solution A, pelleted by low speed centrifugation (500 × g, 1 min) and dissolved S] equation to calculate Ki constant; Ki is defined as a dissociation constant of in 100 μl of 90% (v/v) ethanol. The mixture was treated for 5 min at 80–85 °C, held the enzyme-inhibitor complex. The enzyme activity was monitored at 410 nm on ice for 1 h and centrifuged (4000 × g, 10 min). The supernatant was recovered, and Dixon plots were used to determine Ki under the pseudo-first-order rate the pellet washed twice in 100 μl of re-distilled water (conductance 24 μS), and the limiting conditions by a non-linear regression analysis, as described previously22. supernatants were combined and freeze dried. The extract obtained from crystals, a freeze-dried mother liquor, and a mother liquor with authentic Glc (100 ng) being added, were dissolved in 25 μl 65% (v/v) CH3CN and separated by high perfor- Expression and crystallisation of recombinant HvExoI. The optimised HvExoI mance liquid chromatography (HPLC) coupled with the Evaporative Light Scat- cDNA (GenBank Accession GU441535; was subcloned into the pPICZαBNH8 tering Detector (ELSD model 800, Alltech Associates Inc.). The samples were vector, from which the protein was expressed in Pichia pastoris strain SMD1168H eluted isocratically with 65% (v/v) CH3CN at a flow rate of 0.6 ml/min. The (Invitrogen, Carlsbad, CA, USA). The R158A/E161A variant (primers listed in fractions eluting between 6–7.3 min were collected and freeze dried. A model 1090 Supplementary Table 6) was prepared by site-directed mutagenesis and expressed HPLC liquid system controlled by ChemStation software (Agilent Technologies) in Pichia24. Recombinant HvExoI enzymes were purified using SP-Sepharose was used for chromatography, and separation was carried out on a Prevail NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z Carbohydrate ES column, 5 μm, 250 × 4.6 mm (Alltech Associates Inc.). The conformations and protein-sugar interactions were consistently maintained during column temperature during separations was 21 °C. The column eluent was split classical MD. in 6.5 (collect): 1 (detector) ratio. The detector drift tube was at 40 °C, the nitrogen inlet pressure was 1.5 bars and the signal GAIN setting was set at 4. The HPLC-eluted fractions (6–7.3 min) were reduced with 1 M NaBD4 in 2 M NH4OH QM/MM MD simulations. Classical MD snapshots of two enzyme-product for 18 h at 4 °C and acetylated with acetic anhydride. Myo-inositol (20 ng) was complexes were submitted to QM/MM MD simulations to re-equilibrate the added as an internal standard. Acetylated alditols were analysed on a low polarity system within the QM/MM Hamiltonian. The method by Laio and co-workers81 25 m × 0.25 mm (internal diameter) Chrompack Capillary Column CP-Sil 5 CB was used, which combines Car-Parrinello MD82 with force field-based MD. Only LB/MS (Varian Inc., CA, USA) using He as a gas flow with a Hewlett-Packard the atoms of Glc were included in the QM region (24 atoms), whereas the 6890 Series GC System and a Hewlett-Packard 5973 mass selective detector remainder of the protein and solvent were in the MM region. The MM system was (Agilent Technologies)72. treated with the Amber force field, for consistency with the previous classical MD simulations. The QM system was treated with Density Functional Theory (DFT), using a plane wave (PW) basis set with a kinetic energy cut-off of 70 Ry, Martins- NMR spectroscopy of Glc bound to HvExoI. NMR spectra were acquired at 283 K Troullier pseudopotentials83 and the Perdew-Burke-Ernzerhoff (PBE) exchange- on the Bruker AVANCE III 600 MHz and 800 MHz spectrometers (Bruker, MA, correlation functional84; this is at the same level theory as that described for sugar USA). For Saturation Transfer Difference (STD) experiments, a stock solution conformational Free Energy Landscapes (FELs)33,34. The QM system was enclosed of recombinant HvExoI deglycosylated by endoglycosidase H25 in milliQ water in an isolated orthorhombic box of size 14.0 Å × 13.2 Å × 13.5 Å (−1 subsite) or (143 µM, 20 mM sodium acetate buffer, pH 5.25, 159 mM NaCl) was diluted to 12.0 Å × 15.3 Å × 13.2 Å (+1 subsite). The electronic fictitious was set at 700 au and 40 µM in the same buffer in D2O, and 60 equivalents (2.4 mM) of Glc from a stock the simulation time step was at 0.12 fs. Both complexes were equilibrated at 300 K solution in D2O were added. For transferred Nuclear Overhauser Effect Spectro- for around 5 ps before starting metadynamics simulations of Glc conformational scopy (trNOESY) experiments, the same recombinant enzyme stock solution was FEL. A similar procedure was used in previous studies of conformational diluted to 57 µM with a deuterated 20 mM sodium acetate buffer, pH 5.25, to preferences of sugars in carbohydrate-active enzymes34. which five equivalents of Glc (0.285 mM) were added followed by 300 msec mixing time. In a separate experiment, to the 57 µM recombinant enzyme stock solution, Ab initio metadynamics simulations. QM/MM metadynamics85,86 simulations five equivalents of Glc and three equivalents of thiocellobiose (0.171 mM) were were performed to explore the conformational FEL of Glc bound in the −1 or added, both from a stock solution in D2O, followed by 300 msec mixing time. STD the +1 subsites. We employed the metadynamics driver provided by the Plumed2 spectra were acquired with 2048 scans, using 30 ms Gaussian shaped pulses for plugin, combined with the well-tempered metadynamics approach87. Ring- selective protein saturation at 0.65 ppm; the reference spectrum employed off- puckering coordinates qx and qy were used as collective variables, as described resonant saturation at 100 ppm. Different saturation times (0.5, 1, 2, 3 and 5 s) were for isolated Glc33. Metadynamics parameters were selected from the analysis of used. Prior to Free Induction Decay acquisition, an excitation sculpting double evolution of the puckering coordinates during the unbiased QM/MM MD simu- echo for water suppression and a T2 filter (100 ms XY16 with 500 μs echo times lation and from expected energy barrier heights based on our previous work to minimise dephasing from JHH coupling evolution) protein signal suppression on both isolated Glc and a β-D-glucoside bound in the active site of 1,3;1,4-β-D- was inserted. STD spectra for free Glc were acquired under the same conditions glucanase34. The height and the width of Gaussian terms were initially set at 1 kcal to show that no direct ligand saturation occurred. trNOESY spectra were acquired mol−1 and 0.1 Å, respectively, and a bias factor of ɣ = 10 was used for the well- with 64 scans, 4 K points in the direct dimension, 256 points in the indirect tempered approach. The deposition time was set at 24 fs (200 MD steps). The same dimension and 300 ms mixing time. For water suppression we employed a set up was used for both enzyme-product complexes. Simulations were stopped WATERGATE scheme using 3–9–19 binomial pulses. 1H signal assignments after adding 7,668 Gaussians (enzyme-product complex with Glc bound in the of ligands were derived from two-dimensional TOtal Correlation SpectroscopY −1 subsite) and 7725 Gaussians (Glc bound in +1 subsite). Analyses of evolution (TOCSY), NOESY and 13C-Heteronuclear Single Quantum Coherence (HSQC) of puckering coordinates and FEL showed that the metadynamics simulation spectra. The NOESY spectrum of free thiocellobiose was acquired at 313 K with reached a diffusive regime and that the energy differences converged. 700 ms mixing time to obtain positive NOEs. A NOESY spectrum acquired under the conditions used for the trNOESY spectrum (800 MHz, 283 K, 300 ms mixing time) yielded no observable NOEs, indicating NOE zero-crossing for thiocellobiose Computational methods to explore the Glc exit path in HvExoI. MD simula- (for τc = 1.12/ω0). tions in complex with Glc trapped in the −1 subsite (PDB 3WLH) proceeded as follows. Hydrogen atoms were added and protonation states of titrable residues at pH 7.0 were assigned with PROPKA3.188 through PDB2PQR version 2.0.089. Classical MD simulations. The structure of native HvExoI in complex with Glc In this complex, Glu491 was modelled in a non-ionised form; this is the required (defined as β-D-glucopyranose) (PDB 3WLH) with occupancy 0.5 for −1 and protonation state for the acid/base catalyst in substrate/product complexes of +1 subsites, suggests that the two subsites are alternatively occupied. We set up retaining GH3 enzymes that catalyse hydrolysis by the generally accepted double- two enzyme-product complexes, with Glc bound in the −1 subsite, whereas the displacement reaction mechanism22. The entire system was solvated in a water box +1 subsite was occupied by water molecules, and vice versa. The protonation state of 103.0 × 88.0 × 101.0 Å3 volume, which ensured that there was a 5 Å thick layer of of the titratable residues was selected based on their hydrogen bonding environ- water molecules in each direction from any atom. Sodium ions were added to ment. Namely, all Arg, Lys, Asp and Glu residues were in ionised form, except for neutralise the total charge of the system forming a fully solvated model of 83,428 the acid-base residue Glu491, which was protonated (i.e. non-ionised) to represent atoms. The system was energy-minimised to remove steric clashes and heated the catalytic acid/base form required in the first step of the reaction, and the through 4 ns equilibration (integration steps of 2 fs) until reaching the temperature predominant form at the optimal pH of the enzyme. Histidine residues 98, 111, of 300 K at 1 Bar. This process also involved the application of positional restraints 234, 331, 419, 487 and 586 were considered as Nε-protonated and His207, His262 on the protein and ligand atoms that were gradually released until the whole and His377 as Nδ-protonated. Assigned titration states agreed with pKa values system was fully equilibrated. Equilibration of the system was stopped when root- predicted by PROPKA3.1 and with the optimised hydrogen placement algorithm of mean-square deviation (RMSD) values of the main chain reached the stable value MolProbity. Disulfide bridges were considered between Cys151 and Cys159, and lower than 1.60 Å. Two independent simulations of 100 ns each were run under between Cys513 and Cys518. Three sodium ions were included to achieve neu- NPT conditions, whereas temperature and pressure were controlled through trality, in addition to 150 mM NaCl to mimic the physiological ionic strength. The Langevin dynamics and the Nosé-Hoover algorithm combined with the Langevin system size consisted of approximately 150,000 atoms. The initial structures were piston method79,80, respectively. Periodic boundary conditions were applied, and equilibrated by classical MD simulations. Protein and Glc were described via the the SHAKE algorithm78 was used to adjust O-H separations of water molecules. parm99SB73 and the Glycam0674 force fields, respectively, whereas the TIP3P The cut-off value of 12 Å was used for non-bonded interactions. All calculations model70 was used for water molecules and Å qvist parameters for the sodium and were carried out with NAMD2.976,77 using the CHARMM22 force field84,85, for chloride ions. The classical MD simulations were performed with NAMD2.976,77. the protein and TIP3P for water molecules75, while the CHARMM36 all-atom The SHAKE/RATTLE78 algorithm was used to constrain the bonds involving carbohydrate force field90 was used for Glc. hydrogen atoms in both the solvent and the solute. Multiple time step integration was carried out in r-RESPA; the production simulations used a base time step of 2 fs and a secondary time step of 4 fs for long-range interactions. Non-bonded Docking of Glc and disaccharides in HvExoI. Molecular docking calculations of interactions were computed at every step with a cut-off separation of 9 Å. Periodic Glc and β-D-glucopyranosyl-(1,2)-D-glucose (G2OG) or β-D-glucopyranosyl-(1,3)- boundary conditions and the particle mesh Ewald method were employed to D-glucose (G3OG) to HvExoI for selected protein structures were performed with evaluate long-range electrostatic interactions. Simulations were carried out in the Gold5.291 applying a search cavity of the 20 Å radius and using the Goldscore NPT (N, number of particles; P, system pressure; T, temperature constant) scoring function. Side-chain flexibility for Glu491 was allowed, whereas Trp434 ensemble, with a Nosé-Hoover thermostat and a Langevin piston Nosé-Hoover and Trp286 were flexible as specified below. Complex 1 was generated by docking barostat79,80 to maintain the temperature at 300 K and the pressure at 1 atm, Glc to the HvExoI:G2OG complex derived from the HvExoI:G2SG-OMe crystal respectively. Each of the two enzyme-product complex simulations was run for structure (PDB 6MD6). Complex 2 was generated by docking G2OG to the 40 ns without restraints, ensuring that the system reached equilibrium. Analyses HvExoI:Glc complex derived from the native HvExoI crystal structure (PDB of trajectories used VMD and in-house scripts. One representative snapshot from 3WLH). Both complexes contain Glc bound at the −1 subsite and G2OG at +1 the last time window (≈5 ns) of each simulation was taken as the starting structure and putative +2 subsites. Complexes 3–6 were obtained by docking G2OG or for further QM/MM MD simulations. Crystallographically observed 4C1 Glc G3OG to the HvExoI:Glc structure (PDB 3WLH) subjected to MD simulations, 14 NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z ARTICLE with the major conformation of Trp434 (CA-CB-CG-CD1 dihedral angle of References around −50°). MD simulations were carried out under periodic boundary condi- 1. Warshel, A. et al. Electrostatic basis for enzyme catalysis. Chem. Rev. 10, tions of NPT ensembles, as described for MD simulations of HvExoI in complex 3210–3235 (2006). with Glc. Docking of G2OG, G3OG and β-D-glucopyranosyl-(1,6)-D-glucose 2. Adamczyk, A. J., Cao, J., Kamerlin, S. C. & Warshel, A. Catalysis by (G6OG) lacking the Glc product in the active site was carried out as dihydrofolate reductase and other enzymes arises from electrostatic described above. preorganization, not conformational motions. Proc. Natl Acad. Sci. USA 108, 14115–14120 (2011). GaudiMM calculations. Ternary complexes 1–3 before MD simulation (as 3. Fried, S. D. & Boxer, S. G. Thermodynamic framework for identifying free described above) and complex 3 after MD simulation with substrates bound in a energy inventories of enzyme catalytic cycles. Proc. Natl Acad. Sci. USA 110, productive mode, were used as starting points to investigate Glc displacement 12271–12276 (2013). pathways from the −1 subsite of HvExoI. The protein contained Glc at the 4. Marc, A. & Engasser, J. M. Influence of substrate and product diffusion on −1 subsite, and G2OG or G3OG attached at +1 and putative +2 subsites. heterogenous kinetics of enzymic reversible reactions. J. Theor. Biol. 94, GaudiMM20 (https://github.com/insilichem/gaudi), a recently developed modular 179–189 (1982). multi-objective genetic algorithm platform allows conformational exploration of 5. Elber, R. Long-timescale simulation methods. Curr. Opin. Struct. Biol. 15, defined genes (defined below) with multiple evaluation operators of fitness or 151–156 (2005). objectives. The GPathFinder extension was used, defining as genes protein and 6. Varghese, J. N., Hrmova, M. & Fincher, G. B. Three-dimensional structure of substrates molecules, rotational bonds of ligands and rotamers of the protein, and a barley β-D-glucan exohydrolase, a family 3 glycosyl hydrolase. Structure 7, as objectives, minimisation of van der Waals contacts and maximisation of 179–190 (1999). separations between Glc geometric centres at the −1 subsite. This afforded the low- 7. Hrmova, M. et al. Catalytic mechanisms and reaction intermediates along cost computational pre-identification of putative exit channels. The ProDy (Protein the hydrolytic pathway of plant β-D-glucan glucohydrolase. Structure 9, Dynamics)92 Normal Mode Algorithm (NMA) was used to generate alternative 1005–1016 (2001). starting structures for the GPathFinder. Normal mode calculations were performed 8. Hrmova, M. & Fincher, G. B. Barley β-D-glucan exohydrolases substrate with Tangram implementation on the UCSF Chimera suite93, using residues specificity and kinetic properties. Carbohydr. Res. 305, 209–221 (1998). grouped in clusters of 100 (to avoid protein distortion and achieve acceptable 9. Hrmova, M. & Fincher, G. B. Dissecting the catalytic mechanism of a plant conformations), generating 10 modes with a cut-off of 15.0 Å and a Gamma LJ of beta-D-glucan glucohydrolase through structural biology using inhibitors 1.0. In the case of the crystal structure (PDB 3WLH), the second frame of the mode and substrate analogues. Carbohydr. Res. 342, 1613–1623 (2007). with frequency 31.2 cm−1 was chosen, while the second frame of the mode with 10. Yoshida, E. et al. Role of a PA14 domain in determining substrate specificity of frequency 28.7 cm−1 was selected for complex 3. In both cases, the displacement a glycoside hydrolase family 3 β-glucosidase from Kluyveromyces marxianus. parameter in the Tangram interface was set to 150. Biochem. J. 431, 39–49 (2010). 11. McAndrew, R. P. et al. From soil to structure, a novel dimeric β-glucosidase PELE calculations. The PELE software (https://pele.bsc.es)21,94 was used to belonging to glycoside hydrolase family 3 isolated from compost using simulate Glc displacement and to identify the Glc migration route from the metagenomic analysis. J. Biol. Chem. 288, 14985–14992 (2013). −1 subsite. PELE is a Monte Carlo based algorithm that generates configurations 12. Bacik, J.-P., Whitworth, G. E., Stubbs, K. A., Vocadlo, D. J. & Mark, B. L. along a path through a sequential ligand and protein geometric perturbation Active site plasticity within the glycoside hydrolase NagZ underlies a scheme, side-chain conformational prediction and minimisation steps. The ligand dynamic mechanism of substrate distortion. Chem. Biol. 19, 1471–1482 (Glc) is perturbed by random translations and rotations, and by internal rotations (2012). of protein rotatable bonds, while Cα carbons of the protein are displaced following 13. Suzuki, K. et al. Crystal structures of glycoside hydrolase family 3 β- one of the six lowest Anisotropic Normal Modes95. The algorithm optimally glucosidase 1 from Aspergillus aculeatus. Biochem. J. 452, 211–221 (2013). arranges sidechains surrounding perturbed atoms using a rotameric library side- 14. Pachl, P. et al. Crystal structure of native α-L-rhamnosidase from Aspergillus chain optimisation and energy minimisation with the OPLS-AA force field96. The terreus. Acta Crystallogr. D74, 1078–1084 (2018). move of the system is accepted, defining a minimum, or rejected according to a 15. Hrmova, M. et al. Three-dimensional structure of the barley β-D-glucan Metropolis criterion for a given temperature. The PELE scheme performs an glucohydrolase in complex with a transition state mimic. J. Biol. Chem. 279, effective exploration of the protein energy landscape97–99. Here, the exit path is 4970–4980 (2004). defined by separations between the centre of masses of Glc molecules at the 16. Hrmova, M. et al. Structural rationale for low nanomolar binding of transition beginning (−1 subsite) and at the end of simulations (in solvent). Complexes 1–3 state mimics to a family GH3 β-d-glucan glucohydrolase from barley. with Glc at the −1 subsite and G2OG or G3OG at +1 and putative +2 subsites Biochemistry 44, 16529–16539 (2005). were selected from MD simulations to identify potential exit pathways. One 17. Hrmova, M. et al. Structural basis for a broad specificity in higher plant representative structure from each MD simulation of complexes 1–3 was selected β-D-glucan glucohydrolases. Plant Cell 14, 1–22 (2002). as these complexes show G2OG or G3OG bound in the productive modes for 18. Lombard, V., Golaconda Ramulu, H., Drula, E., Coutinho, P. M. & Henrissat, subsequent catalysis. PELE simulations were terminated, when respective separa- B. The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic Acids tions between centres of masses and C4-OH groups of Glc molecules in initial and final structures were above 15 Å. Res. 42, D490–D495 (2014). 19. Jones, G., Willett, P., Glen, R. C., Leach, A. R. & Taylor, R. Development and validation of a genetic algorithm for flexible docking. J. Mol. Biol. 267, Reporting summary. Further information on experimental design is available 727–748 (1997). in the Nature Research Reporting Summary linked to this article. 20. Rodríguez-Guerra Pedregal, J., Sciortino, G., Guasp, J., Municoy, M. & Maréchal, J.-D. GaudiMM: a modular multi-objective platform for molecular Data availability modeling. J. Comput. Chem. 38, 2118–2126 (2017). The atomic coordinates and structure factors were deposited in the Protein Data Bank 21. Borrelli, K. W., Vitalis, A., Alcantara, R. & Guallar, V. PELE: protein energy (www.pdb.org) with the PDB accessions: native HvExoI, and in complex with 3dGlc, landscape exploration. a novel Monte Carlo based technique. J. Chem. Theory. 4dGlc, octyl-O-Glc, octyl-S-Glc and PEG400 are 3WLH, 3WLJ, 3WLK, 3WLM, 3WLN Comput. 1, 1304–1311 (2005). and 3WLL, respectively. The PDB accessions of WT recombinant ligand-free HvExoI 22. Hrmova, M. et al. Barley β-D-glucan exohydrolases with β-D-glucosidase and in complex with Glc, G2SG-OMe and G6SG-OMe are 3WLI, 3WLO, 6MD6 and activity. Purification and determination of primary structure from a cDNA 3WLP, respectively. The PDB accessions of recombinant R158A/E161A ligand-free clone. J. Biol. Chem. 271, 5277–5286 (1996). HvExoI and in complex with Glc and G6SG-OMe, are 3WLQ, 3WLR and 6MI1, 23. Hrmova, M., Varghese, J. N., Høj, P. B. & Fincher, G. B. Crystallization respectively. A reporting summary for this Article is available as a Supplementary and preliminary X-ray analysis of β-glucan exohydrolase isoenzyme ExoI Information file. All other data generated and/or analysed during the current study are from barley (Hordeum vulgare). Acta Crystallogr. D54, 6887–689 (1998). available from the corresponding author upon reasonable request. 24. Luang, S., Hrmova, M. & Ketudat Cairns, J. R. High-level expression of barley β-D-glucan exohydrolase HvExoI from a codon-optimized cDNA in Pichia pastoris. Prot. Exp. Purif 73, 90–98 (2010). Code availability 25. Luang, S., Ketudat Cairns, J. R., Streltsov, V. A. & Hrmova, M. Crystallisation Software applications described in Methods with associated references were used without of wild-type and variant forms of a recombinant β-D-glucan glucohydrolase code modifications. For analyses of geometrical parameters along MD simulations, we from barley (Hordeum vulgare L.) by macroseeding with wild-type native have used an in-house script, which is available upon request from the corresponding microcrystals and preliminary X-ray analysis. Int. J. Mol. Sci. 11, 2759–2769 author. (2010). 26. Cremer, D. & Pople, J. A. A general definition of ring puckering coordinates. J. Am. Chem. Soc. 97, 1354–1358 (1975). Received: 10 December 2018 Accepted: 22 March 2019 27. Quiocho, F. A. Carbohydrate-binding proteins: tertiary structures and protein-sugar interactions. Ann. Rev. Biochem 55, 287–315 (1986). NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications 15 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z 28. Fernandez-Alonso, M. C., Canada, F. J., Jimenez-Barbero, J. & Cuevas, G. insights into the catalytic machinery of family I glycosidases. J. Mol. Biol. 371, Molecular recognition of saccharides by proteins. Insights on the origin of the 1204–1218 (2006). carbohydrate-aromatic interactions. J. Am. Chem. Soc. 127, 7379–7386 (2005). 56. Horton, D. 1-Thio-D-glucose. Methods Carbohydr. Chem 2, 433–437 (1963). 29. Fersht, A. Structure and Mechanism in Protein Science, pp. 1-631, (W. H. 57. Classon, B., Liu, Z. & Samuelsson, B. New halogenation reagent systems Freeman and Co, New York, 1999). useful for the mild one-step conversion of alcohols into iodides or bromides. 30. Mayer, M. & Meyer, B. Characterization of ligand binding by saturation J. Org. Chem. 53, 6126–6130 (1998). transfer difference NMR spectroscopy. Angew. Chem. Int. Ed. Engl. 38, 58. Dong, H., Pei, Z., Angelin, M., Byström, S. & Ramström, O. Efficient synthesis 1784–1788 (1999). of β-D-mannosides and β-D-talosides by double parallel or double serial 31. Clore, G. M. & Gronenborn, A. M. Theory and applications of the transferred inversion. J. Org. Chem. 72, 3694–3701 (2007). Nuclear Overhauser Effect to the study of the conformations of small ligands 59. Horton, D., Wolfrom, M. L. & Thiosugars., I. Synthesis of derivatives of bound to proteins. J. Magn. Reson. 48, 402–417 (1982). 2-amino-2-deoxy-1-thio-D-glucose. J. Org. Chem. 27, 1794–1800 (1962). 32. Montero, E. et al. NMR studies of the conformation of thiocellobiose bound 60. Johnston, B. D. & Pinto, B. M. Synthesis of thio-linked disaccharides by 1 → 2 to a beta-glucosidase from Streptomyces sp. FEBS Lett. 421, 243–248 (1998). intramolecular thioglycosyl migration: oxacarbenium versus episulfonium ion 33. Biarnés, X., Ardèvol, A., Iglesias-Fernández, J., Planas, A. & Rovira, C. intermediates. J. Org. Chem. 65, 4607–4617 (2000). Catalytic itinerary in 1,3-1,4-β-glucanase unraveled by QM/MM 61. Otwinowski, Z. & Minor, W. A processing of X-ray diffraction data collected metadynamics. Charge is not yet fully developed at the oxocarbenium ion-like in oscillation mode. Methods Enzymol. 276, 307–326 (1996). transition state. J. Am. Chem. Soc. 133, 20301–20309 (2011). 62. Marshudov, G. N. P. et al. REFMAC5 for the refinement of macromolecular 34. Ardèvol, A. & Rovira, C. Reaction mechanisms in carbohydrate-active crystal structures. Acta Crystallogr. D67, 355–367 (2011). enzymes: glycoside hydrolases and glycosyltransferases. Insights from ab initio 63. Adams, P. D. et al. PHENIX: a comprehensive Python-based system for quantum mechanics/molecular mechanics dynamic simulations. J. Am. Chem. macromolecular structure solution. Acta Crystallogr. D66, 213–221 (2010). Soc. 137, 7528–7547 (2015). 64. McRee, D. E. XtalView/Xfit - a versatile program for manipulating atomic 35. Lee, P. I. & Peppas, N. A. Prediction of polymer dissolution in swellable coordinates and electron density. J. Struct. Biol. 125, 156–165 (1996). controlled release systems. J. Control Release 6, 207–215 (1987). 65. Emsley, P., Lohkamp, B., Scott, W. G. & Cowtan, K. Features and development 36. Mills, J. E. & Dean, P. M. Three-dimensional hydrogen-bond geometry and of coot. Acta Crystallogr. D66, 486–501 (2010). probability information from a crystal survey. J. Comput. Aided. Mol. Des. 10, 66. Winn, M. D., Isupov, M. N. & Murshudov, G. N. Use of TLS parameters 607–622 (1966). to model anisotropic displacements in macromolecular refinement. 37. Nashine, V. C., Hammes-Schiffer, D. & Benkovic, S. J. Coupled motions Acta Crystallogr. D57, 122–133 (2001). in enzyme catalysis. Curr. Opin. Chem. Biol. 14, 644–651 (2010). 67. Painter, J. & Merritt, E. A. A molecular viewer for the analysis of TLS rigid- 38. Buyong, M. & Nussinov, R. Enzyme dynamics point to stepwise body motion in macromolecules. Acta Crystallogr. D61, 465–471 (2005). conformational selection in catalysis. Curr. Opin. Chem. Biol. 14, 652–659 68. Brünger, A. T. Free R value: a novel statistical quantity for assessing the (2010). accuracy of crystal structures. Nature 355, 472–475 (1992). 39. Lampe, J. N., Brandman, R., Sivaramakrishnan, S. & de Montellano, P. R. 69. Brünger, A. T. X-PLOR Version 3.851: a system for crystallography and NMR. Two-dimensional NMR and all-atom molecular dynamics of cytochrome (Yale University Press, New Haven, CT, 1992). P450 CYP119 reveal hidden conformational substates. J. Biol. Chem. 285, 70. Ramachandran, G. N., Ramakrishnan, C. & Sasisekharan, V. Stereochemistry 9594–9603 (2010). of polypeptide chain configurations. J. Mol. Biol. 7, 95–99 (1963). 40. Pieraccini, S., Sironi, M. & Colombo, G. A molecular simulation analysis 71. Laskowski, R. A., MacArthur, M. W., Moss, D. S. & Thornton, J. M. of the origins of regioselectivity. Chem. Physics Lett 418, 373–376 (2006). PROCHECK—a program to check the stereochemical quality of protein 41. Sundaresan, V., Chartron, J., Yamaguchi, M. & Stout, C. D. Conformational structures. J Appl. Cryst. 26, 283–291 (1993). diversity in NAD(H) and interacting transhydrogenase nicotinamide 72. Lau, E. & Bacic, A. Capillary gas chromatography of partially methylated nucleotide binding domains. J. Mol. Biol. 346, 617–629 (2004). alditol acetates on a high polarity, cross-linked, fused silica BPX70. Column. 42. Yang, L.-W. & Bahar, I. Dynamic coupling between the SH2 and SH3 domains J. Chrom. 637, 100–103 (1993). of c-Src and Hck underlies their inactivation by C-terminal tyrosine 73. Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A. & Case, D. A. phosphorylation. Structure 13, 893–904 (2005). Development and testing of a general amber force field. J. Comput. Chem. 25, 43. Gutteridge, A. & Thornton, J. Conformational changes observed in enzyme 1157–1174 (2004). crystal structures upon substrate binding. J. Mol. Biol. 346, 21–28 (2005). 74. Kirschner, K. N. et al. GLYCAM06: a generalizable biomolecular force field. 44. Hammes-Schiffer, S. & Benkovic, S. J. Relating protein motion to catalysis. Carbohydrates. J. Comput. Chem. 29, 622–655 (2008). Annu. Rev. Biochem. 75, 519–541 (2006). 75. Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. & 45. Díaz, M. D., Fernández-Alonso, M. C., Cuevas, G., Cañada, F. J. & Jiménez- Klein, M. L. Comparison of simple potential functions for simulating Barbero, J. On the role of aromatic-sugar interactions in the molecular liquid water. J. Chem. Phys. 79, 926–935 (1983). recognition of carbohydrates: A 3D view by using NMR. Pure Appl. Chem. 80, 76. Kale, L. et al. NAMD2: greater scalability for parallel molecular dynamics. 1827–1835 (2008). Comput. Phys. 151, 283–312 (1999). 46. Hrmova, M. et al. Barley (1,3)-β-D-glucan endohydrolase mutants synthesise 77. Phillips, J. C. et al. Scalable molecular dynamics with NAMD. Comput. Chem. crystalline (1,3)-β-D-glucans. J. Biol. Chem. 277, 30102–30111 (2002). 26, 1781–1802 (2005). 47. Davies, G. & Henrissat, B. Structures and mechanisms of glycosyl hydrolases. 78. Ryckaert, J. P., Ciccotti, G. & Berendsen, H. J. C. Numerical integration of Structure 3, 853–859 (1995). the Cartesian equations of motion of a system with constraints: molecular 48. Sørensen, A., Lübeck, M., Lübeck, P. S. & Ahring, B. K. Fungal beta- dynamics of n-alkanes. Comput. Phys 23, 327–341 (1977). glucosidases: a bottleneck in industrial use of lignocellulosic materials. 79. Nosé, S. J. A unified formulation of the constant temperature molecular Biomolecules 3, 612–631 (2013). dynamics methods. Chem. Phys. 81, 511–519 (1984). 49. van Dongen, S. F. M., Elemans, J. A., Rowan, A. E. & Nolte, R. J. Processive 80. Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. catalysis. Angew. Chem. Int. Ed. 52, 11420–11428 (2014). Phys. Rev. A. 31, 1695–1697 (1985). 50. Sørlie, M., Zakariassen, H., Norberg, A. L. & Eijsink, V. G. H. Processivity and 81. Laio, A., VandeVondele, J. & Rothlisberger, U. A. Hamiltonian electrostatic substrate-binding in family 18 chitinases. Biocatal. Biotransform. 30, 353–365 coupling scheme for hybrid Car-Parrinello molecular dynamics simulations. (2012). J. Chem. Phys. 116, 6941–6947 (2002). 51. Hrmova, M. & Fincher, G. B. Structure-function relationships of β-D-glucan 82. Car, R. & Parrinello, M. Unified approach for molecular dynamics and endo- and exohydrolases from higher plants. Plant Mol. Biol. 47, 73–91 density-functional theory. Phys. Rev. Lett. 55, 471–2474 (1985). (2001). 83. Troullier, N. & Martins, J. L. Efficient pseudopotentials for plane-wave 52. Varghese, J. N. et al. Three-dimensional structures of two plant β-glucan calculations. Phys. Rev. B 43, 1993–2006 (1991). endohydrolases with distinct substrate specificities. Proc. Natl Acad. Sci. USA 84. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation 91, 2785–2789 (1994). made simple. Phys. Rev. Lett. 77, 3865–3868 (1996). 53. Kleywegt, G., Zou, J.-Y., Divne, C. & Davies, G. J. & Sinning, I. et al. The 85. Barducci, A., Bonomi, M. & Parrinello, M. Metadynamics. WIREs Comput. crystal structure of the catalytic core domain of endoglucanase I from Mol. Sci. 1, 826–843 (2011). Trichoderma reesei at 3.6 Å resolution, and a comparison with related 86. Laio, A. & Parrinello, M. Escaping free-energy minima. Proc.Natl Acad. Sci. enzymes. J. Mol. Biol. 273, 1–15 (1997). USA 99, 12562–12566 (2002). 54. Divne, C., Stahlberg, J., Teeri, T. T. & Jones, T. A. High-resolution crystal 87. Barducci, A., Bussi, G. & Parrinello, M. Well-tempered metadynamics: a structures reveal how a cellulose chain is bound in the 50 Å long tunnel of smoothly converging and tunable free-energy method. Phys. Rev. Lett. 100, cellobiohydrolase I from Trichoderma reesei. J. Mol. Biol. 275, 309–325 (1998). 020603 (2008). 55. Isorna, P. et al. Crystal structures of Paenibacillus polymyxa β-glucosidase B 88. MacKerell, A. D. Jr. et al. All-atom empirical potential for molecular modeling complexes reveal the molecular basis of substrate specificity and give new and dynamics studies of proteins. J. Phys. Chem. B 102, 3586–3616 (1998). 16 NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09691-z ARTICLE 89. MacKerell, A. D. Jr., Feig, M. & Brooks, I. I. I. C. L. Extending the treatment 72130118) and J.R.-G. to Generalitat de Catalunya and the European Social Fund of backbone energetics in protein force fields: limitations of gas-phase (2017FI_B2_00168). S.F. was supported by Glyco@Alps (ANR-15-IDEX-02), Carnot quantum mechanics in reproducing protein conformational distributions Institut PolyNat, Labex ARCANE and CBH-EUR-GS (ANR-17-EURE-0003), and thanks in molecular dynamics simulations. J. Comput. Chem. 25, 1400–1415 to Chimie Moléculaire de Grenoble for access to facilities. M.A.-P. and C.R. acknowledge (2004). MareNostrum and Minotauro and BSC-CNS (RES-QCM-2017–2–001) resources. Use 90. Guvench, O. et al. CHARMM additive all-atom force field for carbohydrate of synchrotrons was supported by the Australian Synchrotron Research Program, which derivatives and their utility in polysaccharide and carbohydrate-protein is funded by the Commonwealth of Australia under the Major National Research modeling. J. Chem. Theor. Comp 7, 3162–3180 (2011). Facilities Program. 91. Jones, G., Willett, P., Glen, R. C., Leach, A. R. & Taylor, R. Development and validation of a genetic algorithm for flexible docking. J. Mol. Biol. 267, 727–748 (1997). Author contributions S.F. synthesised thio-oligosaccharides; M.Hrmova. performed enzyme inhibition kinetics; 92. Bakan, A., Meireles, L. M. & Bahar, I. ProDy: protein dynamics inferred M.Hijnen. conducted SPR; S.L. and J.R.K.C. carried out mutagenesis; S.L. and M.H. from theory and experiments. Bioinformatics 27, 1575–1577 (2011). 93. Pettersen, E. F. et al. UCSF Chimera—a visualization system for exploratory performed crystallisation; V.A.S., J.N.V. and M.Hrmova. collected X-ray data; V.A.S., research and analysis. J. Comput. Chem. 25, 1605–1612 (2004). S.L., A.P. and M.Hrmova. refined crystal structures; M.Hrmova. executed GC/MS and HPLC; A.A. and J.J-B. performed NMR spectroscopy; I.T. conducted initial MD 94. Madadkar-Sobhani, A. & Guallar, V. PELE web server: atomistic study of calculations; M.A.-P. and C.R. performed ab initio QM/MM metadynamics and biomolecular systems at your fingertips. Nucleic Acids Res. 41, W322–W328 (2013). conformational FEL calculations; F.M., L.T-S., J.-E.S.-A., J.R-G., J.M.L., J.-D.M. and L.M. 95. Cossins, B. P., Hosseini, A. & Guallar, V. Exploration of protein implemented docking, classical MD simulations, GaudiMM and PELE calculations; M.Hrmova. wrote the manuscript with contributions from co-authors; M.Hrmova. conformational change with PELE and metadynamics. J. Chem. Theory. oversaw the project. Comput. 8, 959–965 (2012). 96. Kaminski, G. A., Friesner, R. A., Tirado-Rives, J. & Jorgensen, W. L. Evaluation and reparametrization of the OPLS-AA force field for proteins via comparison with accurate quantum chemical calculations on peptides. Additional information Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- J. Phys. Chem. B 105, 6474–6487 (2001). 019-09691-z. 97. Hernández-Ortega, A. et al. Substrate diffusion and oxidation in GMC oxidoreductases: an experimental and computational study on fungal Competing interests: The authors declare no competing interests. aryl-alcohol oxidase. Biochem. J. 436, 341–350 (2011). 98. Lucas, F. et al. Molecular determinants for selective C25-hydroxylation of Reprints and permission information is available online at http://npg.nature.com/ vitamins D2 and D3 by fungal peroxygenases. Catal. Sci. Technol 6, 288–295 reprintsandpermissions/ (2016). 99. Gygli, G., Lucas, M. F., Guallar, V. & van Berkel, W. J. H. The ins and outs Journal peer review information: Nature Communications thanks the anonymous of vanillyl alcohol oxidase: Identification of ligand migration paths. reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports PLoS Comput. Biol. 13, e1005787 (2017). are available. Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in Acknowledgements published maps and institutional affiliations. We thank M. Raab (Slovak Academy of Sciences), B.J. Smith (La Trobe University) and G.B. Fincher (University of Adelaide) for interest in this research. H. Tong (Advanced Photon Source), N. Matsugaki and S. Wakatsuki (Photon Factory), and F. Pettolino (University of Melbourne) are thanked for advice and S. Pradeau for technical assistance. Open Access This article is licensed under a Creative Commons Funding was supported by the Huaiyin Normal University and the Australian Research Attribution 4.0 International License, which permits use, sharing, Council (DP120100900) to M.Hrmova., the Australian Synchrotron (MX1 and MX2 adaptation, distribution and reproduction in any medium or format, as long as you give beamlines) to M.Hrmova. and V.A.S., the Advanced Photon Source (14-ID-B beamline) appropriate credit to the original author(s) and the source, provide a link to the Creative to J.N.V. and M.Hrmova., the Photon Factory (BL5 beamline) to V.A.S, the Suranaree Commons license, and indicate if changes were made. The images or other third party University of Technology and the Thailand Research Fund (BRG5980015) to J.R.K.C. material in this article are included in the article’s Creative Commons license, unless and S.L., Generalitat de Catalunya (Commissioner for Universities and Research, indicated otherwise in a credit line to the material. If material is not included in the Department of Innovation, Universities and Enterprise) and the European Union article’s Creative Commons license and your intended use is not permitted by statutory through a Beatriu de Pinoś fellowship (BP-B 2013) to M.A.-P., Generalitat de Catalunya regulation or exceeds the permitted use, you will need to obtain permission directly from (SGR2017–1189), MICINN (CTQ2017–85496-P) and Spanish Structures of Excellence the copyright holder. To view a copy of this license, visit http://creativecommons.org/ María de Maeztu (MDM-2017–0767) to C.R., Spanish Ministerio de Ciencia e Innova- licenses/by/4.0/. ción to J.M.L. (CTQ2017–83745-P) and J.-D.M. (CTQ2017–87889-P), and Generalitat de Catalunya (SGR2017–13234) to J.-D.M. L.M. thanks to the Universitat Autònoma de Barcelona Talent Program, F.M. to CONICYT (PFCHA/Doctorado Becas CHILE/2012 © Crown 2019 NATURE COMMUNICATIONS | (2019)10:2222 | https://doi.org/10.1038/s41467-019-09691-z | www.nature.com/naturecommunications 17