Article https://doi.org/10.1038/s41467-022-33180-5 The evolutionary advantage of an aromatic clamp in plant family 3 glycoside exo-hydrolases Received: 25 May 2022 Sukanya Luang1,10, Xavier Fernández-Luengo2,10, Alba Nin-Hill 3,10, Victor A. Streltsov 4, Julian G. Schwerdt1, Santiago Alonso-Gil3, Accepted: 3 September 2022 James R. Ketudat Cairns 5, Stéphanie Pradeau 6, Sébastien Fort 6 , Jean-Didier Maréchal 2, Laura Masgrau 2,7, Carme Rovira 3,8 & Maria Hrmova 1,9 Check for updates 1234567890():,; 1234567890():,; In the barley β-D-glucan glucohydrolase, a glycoside hydrolase family 3 (GH3) enzyme, the Trp286/Trp434 clamp ensures β-D-glucosides binding, which is fundamental for substrate hydrolysis during plant growth and development. We employ mutagenesis, high-resolution X-ray crystallography, and multi- scale molecular modelling methods to examine the binding and conforma- tional behaviour of isomeric β-D-glucosides during substrate-product assis- ted processive catalysis that operates in GH3 hydrolases. Enzyme kinetics reveals that the W434H mutant retains broad specificity, while W434A behaves as a strict (1,3)-β-D-glucosidase. Investigations of reactant move- ments on the nanoscale reveal that processivity is sensitive to mutation- specific alterations of the tryptophan clamp. While wild-type and W434H utilise a lateral cavity for glucose displacement and sliding of (1,3)-linked hydrolytic products through the catalytic site without dissociation, con- sistent with their high hydrolytic rates, W434A does not adopt processive catalysis. Phylogenomic analyses of GH3 hydrolases disclose the evolu- tionary advantage of the tryptophan clamp that confers broad specificity, high catalytic efficiency, and processivity. The barley (Hordeum vulgare L.) β-D-glucan glucohydrolase contains the catalytic nucleophile Asp285 and the acid/base catalyst HvExoI is the only plant family 3 glycoside hydrolase (GH3) for Glu491, while the +1 subsite is shaped in a clamp with a coin slot which the atomic structure is available1–3, although the formed by the Trp286 and Trp434 residues that hold β-D-glucosides Carbohydrate-Active Enzymes CAZy database4 lists nearly entering the active site in a productive mode for hydrolysis1–3,5 48,400 GH3 entries (August 2022) from archaea, bacteria, fungi, (Fig. 1a, b). HvExoI catalyses the hydrolytic removal of nonreducing and eukaryotes, including economically important plants Triti- glucose (Glc) moiety from a broad spectrum of poly- and oligomeric β- cum, Oryza, Glycine, Zea, Brassica, Medicago, Phaseolus, and D-glucosides. This hydrolytic reaction underpins a vital role of HvExoI Solanum. during plant growth and development, to participate during e.g. pri- HvExoI folds into a two-domain structure consisting of an (α/β)8 mary root extension and the recycling of mixed-linkage (1,4;1,3)-β-D- barrel (domain 1) and an (α/β)6 sandwich (domain 2) with a 13 Å-deep glucans during germination and root elongation1,3,6. The capacity of active site pocket embracing the subsites −1 and +1, positioned at the HvExoI to hydrolyse a panel of β-D-glucosides was rationalised from the interface between the two domains. At the −1 subsite, the pocket observation that non-reducing Glc at the −1 subsite is locked into a A full list of affiliations appears at the end of the paper. e-mail: maria.hrmova@adelaide.edu.au Nature Communications | (2022)13:5577 1 Article https://doi.org/10.1038/s41467-022-33180-5 Asp285 Trp286 a Trp286 b Asp285 Tyr253 C1 Tyr253 C1 His207 His207 C4 C4 Glc Glc -1 subsite -1 subsite Asp95 Asp95 Arg158 Glu491 Trp434 Glu491 Trp434 Arg158 β-D-Glucan glucohydrolases c Subclade 2 Mp bglu III Cp bglu Pp bgl u Sf bg Subclade 1 Mp bbglu I glu II Trp286/Tyr286 Subclade 3 and Rs u xo IIII Trp286 S r bglu At E xo o I lu Ps b lu Gln434/Tyr434/ At t E bglu Trp286 lu I Mp and Ps m bg bg glu II A a E npro and Trp434 Trp434 C A Gly434/His434/Lys434/Trp434 x La I Fv u n p I 1 o bg r o Ex o II 0.8 0.6 lu Es Hv Ex 1 bg 1 lu Zm I pro II 0.6 un 0.5 0.9 At b At unpro 1 glu At 0.3 1 0.8 H3 Cf G 1 II 1 Hv Exo 0. 8 0.2 Cb b 0.9 0.5 0.5 1 Nt Exo glu 1 0.9 Tm bglu Trp286 7 0. and 1 Xc bg lu Trp434 1 Trp286 1 and Trp434 0.9 1 0.9 β-D-Glucosidases/ 0.3 β-D-Glucosidases β-N-Acetylhexosaminidases 0.9 1 5 0. 1 Hj GH 0.5 1 Bs bglu Dc bg Gg ave3 0.6 0.9 0.6 Ni bg lu 0.4 1 0.3 lu Ci bglu 0.9 1 0.8 0.4 Wa bglu Cr G 0.5 1 lu H3 Sc bg 0.7 1 S Pst ac-g 9 0. a l bglu C p h c-g c-ami P 1 Km r bglu3 s e lc- G xo am Cs A GH lu 1 g H i Pb n b 3 bg 1 T lu 1 3 1 H lu G bg 1 Se 1 β-N-Acetylhexosaminidases 0.6 Cg glc ln- D d c bg lu g E b bg Mp b bxyl II V T xyl-ara Te Tb I b g lu Rf GH Mp bxyl III lI Hj bxyl Nc bglu Mp bxyl Hv bxyl lu xy 2.0 Mp 3 β-D-Xylosidases/ β-D-Xylosidases α-L-Arabinofuranosidases Fig. 1 | Exo-hydrolytic enzymes of the GH3 family with an entrapped Glc within 3.3 Å. The positions of the −1 subsites in panels a and b are indicated. molecule and a ligand-free form, and their phylogeny. a Hordeum vulgare c Unrooted phylogenetic tree of representative GH3 family members (accession HvExoI with entrapped Glc (PDB 3WLH; atomic colours), Bacteroides ovatus β- numbers listed in Supplementary Table 1) revealed the presence of four major glucosidase with perfused Glc (PDB 5JP0; green cpk), and a ligand-free form of clades with five activities [β-D-glucan glucohydrolases (EC 3.2.1.-), β-D-glucosidases Pseudoalteromonas sp. exo-(1,4;1,3)-β-glucanase (PDB 3UT0; yellow cpk). Distances (EC 3.2.1.21), β-N-acetylhexosaminidases (EC 3.2.1.52), β-D-xylosidases (EC 3.2.1.37), ached of HvExoI residues to Glc in the −1 subsite (black dashes) are within 3.3 Å. and α-L-arabinofuranosidases (EC 3.2.1.55)]. In β-D-glucan glucohydrolases, three HvExoI and structurally equivalent residues in 3UT0 and 5JP0 superpose with subclades were identified with residue variability in the aromatic clamp (residue respective RMSD values of 0.7 Å and 1.2 Å. b Surface morphology of the active numbering corresponding to HvExoI; PDB 1EX1). The percentage of replicate trees site pocket of HvExoI coloured by an electrostatic potential: white, neutral; red, with taxa clustering together in the bootstrap test (1000 replicates) are shown at −5 kT·e−1. Separations (dashed lines) between Glc and key interacting residues are node bipartitions. fixed position through a network of hydrogen bonds (H-bonds), while was estimated that the strength of binding through aromatic residues the reducing Glc moiety at the +1 subsite is constrained through could be in the range of 2-5 kcal/mol7. stacking interactions, and thus the flexibility of attached substrates Since the HvExoI structure was elucidated1–3, many prokaryotic, and their projection is unencumbered allowing HvExoI to accom- fungal and unclassified metagenome GH3 structures have become modate isomeric (1,2)-, (1,3)-, (1,4)- and (1,6)-linked β-D-glucosides2,3. It available but provided little information on the substrate specificity of Nature Communications | (2022)13:5577 2 Article https://doi.org/10.1038/s41467-022-33180-5 this class of enzymes. These include two-domain bacterial β-N-acet- metadynamics, indicate that in this GH3 exo-hydrolase the aromatic ylglucosaminidases, for example, from Bacillus subtilis8, and a β-N- clamp fulfils the role of a recognition platform ensuring effective acetylhexosaminidase from Salmonella typhimurium9, while a β- carbohydrate-binding due to cooperative rotational and translational hexosaminidase from Burkholderia cenocepacia10 folds into a single- movements. This knowledge expands on the applicability of bioengi- domain (α/β)8 sandwich structure. Several three-domain architectures neering of carbohydrate-processing enzymes to improve catalytic consisting of an (α/β)8 barrel, an (α/β)6 sandwich, and another β- rates, and stability, and to understand product inhibition (important in sandwich domain were revealed in β-glucosidases from Thermotoga drug discovery), as the traditional industrial approaches are being neapolitana11, Hypocrea jecorina12 and Aspergillus aculeatus13, in a β-N- increasingly replaced by enzyme-based technologies. acetylglucosaminidase from Thermotoga neapolitana14, in a (1,4;1,3)-β- exo-glucanase from Pseudoalteromonas sp. BB115, and in a β- Results glucosidase JMB19063 isolated from metagenome16. These additional Mutants altered in Trp286 and Trp434 residues β-sandwich domains are described as fibronectin type III or Ig-like β- The productive binding of carbohydrates in HvExoI involves the par- sandwich folds and have similar architectures. Four-domains struc- ticipation of Trp286 and Trp434 residues that form an aromatic clamp tures consisting of (α/β)8 barrel, an (α/β)6 sandwich, and two β- at the +1-subsite, while the neighbouring subsite contains the Asp285 sandwich domain folds (PA14 and fibronectin type III) were identified and Glu491 catalysts1–3 (Fig. 1a). In HvExoI, after the hydrolytic cycle is in the β-glucosidases from Kluyveromyces marxianus17, Aspergillus completed, both subsites contain entrapped Glc products that are fumigatus, A. oryzae18 and Streptomyces venezuelae19. The precise roles displaced after a new substrate is attached22. The Glc molecule of these additional β-sandwich domains have not been fully clarified. entrapped in the −1 subsite (Fig. 1b) is held at the bottom of the pocket However, in some GH3 enzymes these domains are thought to be by at least 10 residues (Fig. 1a, b). The residues equivalent to the HvExoI linked to substrate binding and specificity17, protein structure 286 and 434 positions overlap in the homologous GH3 exo- stability14,20, and oligomerisation, meaning that dimerisation is hydrolases14,15,32, form similar aromatic clamps, and that of Bacter- required for catalytic activity16. oides β-glucosidase (PDB 5JP0) harbours the Glc molecule after in Although the details of the catalytic mechanism in HvExoI were crystallo perfusion (Fig. 1a). To resolve phylogenetic relationships of reported1–3,21, the precise information on the binding of positionally HvExoI within the GH3 family4 a phylogenetic tree of selected entries isomeric β-D-glucosides at the atomic level remains limited. Recently, was built revealing that the entries segregated in four clusters we envisaged that β-D-glucan glucohydrolases, such as HvExoI, adopt harbouring five enzyme activities: β-D-glucan glucohydrolases, so-called substrate-product assisted processive catalysis, where after β-D-glucosidases, β-N-acetylhexosaminidases, β-D-xylosidases, and α-L- an incoming substrate binds in a productive mode, the Glc product arabinofuranosidases33,34 (Fig. 1c; Supplementary Table 1). The evokes the formation of a transient lateral cavity, which serves as a inspection of the β-D-glucan glucohydrolase cluster revealed that conduit for its egress to allow for the next catalytic cycle22. To our Trp286 and Trp434 remained conserved in subclade 1, while in the knowledge, this type of processive catalysis (or non-dissociative subclades 2 and 3 enzymes, we identified variations in both Trp resi- sequential degradation) discovered in HvExoI22 has not been descri- dues – Trp286/Tyr286 and Gln434/Tyr434/Trp434 in subclade 2, and bed in hydrolases with a closed pocket-shaped catalytic site. In doing Gly434/His434/Lys434/Trp434 in subclade 3 (Fig. 1c), which may so22, we defined processive catalysis (which we termed substrate- impact substrate binding and processivity in these enzymes. product assisted processivity) in broader terms than those used for To examine the precise mechanistic roles of the Trp286 and the definitions of the tunnel-, groove-, ridge-, or cleft-like active sites Trp434 residues in the binding of isomeric (1,2)-, (1,3)-, (1,4)- and of GHs23–25. Considering that HvExoI operates on plant cell wall (1,6)-linked β-D-glucosides in HvExoI, we methodically varied these (1,3;1,4)-β-D-glucan polysaccharides, processive catalysis (as a form of residues to Ala, His, Phe, and Tyr to generate single (and in some efficient catalysis) in this enzyme is expected as rapid substrate instances) double mutants, and investigated their substrate specifi- hydrolysis is critical to embryo development during seed city, steady-state and inhibitions kinetics and substrate binding in germination5,26. Processive catalysis has been defined in (1,3)-β-D- atomic structures (Figs. 2–4; Supplementary Figs. 2–4; Supplemen- glucan endohydrolases27, (1,4)-β-D-glucan endohydrolases28, tary Tables 2–6). 23 24,25,29 cellobiohydrolases , and chitobiohydrolases . Conversely, exo- acting hydrolases21,22 and β-D-glucosidases30 with closed funnel-, Substrate specificity of WT, and Trp286 and Trp434 mutants crater-, or pocket-shaped active sites have been deemed to be We compared relative hydrolytic rates of β-linked poly- and oligo- nonprocessive30. The work with processive enzymes suggested the saccharides, and aryl-glycosides under optimal physico-chemical role of aromatic and polar amino acid residues to assist with stacking conditions (Supplementary Fig. 1; Supplementary Table 2; Supple- and H-bond interactions for polysaccharide sliding motions through mentary Note 1). Of the tested mutants, W434H was the only mutant the catalytic sites22–25,27–29,31. that exhibited broad substrate specificity similar to WT, contrary to In this context, our aim was to provide direct structural evidence W286F and W286Y, and W434F, W434Y, W286F/W434F, and W286F/ for the binding of β-D-glucosides and the unbinding of their hydrolytic W434A that lacked activities of at least with one carbohydrate sub- products and investigate the relationship between reactant move- strate (Supplementary Table 2). The absence of hydrolytic activities ments and processivity. was observed with (1,4;1,3)-β-D-glucans (barley glucan and lichenin) Here, we capitalise on a unique feature of HvExoI to bind a panel and 4NP-cellobioside in all mutants except W434H and a decrease in of isomeric β-D-glucosides through the aromatic clamp of two Trp hydrolytic rates with oligomeric substrates and aryl β-D-glucosides residues, and elucidate, how this precisely facilitates substrate binding, except in W286Y with 4-nitrophenyl β-D-glucopyranoside (4NP-Glc) catalytic efficiency, and substrate-product assisted processivity. We (Supplementary Table 2). Remarkably, W434A lost all hydrolytic present atomic details of HvExoI wild-type (WT) and its Ala, His, Phe, activities with poly- and oligosaccharides, except (1,3)-β-D-linked β-D- and Tyr mutants of the Trp286/Trp434 clamp in the complex with glucosides, turning this mutant into a strict (1,3)-β-D-glucosidase, positionally isomeric thio-analogues. Based on these structures, we albeit with low specific activities (Supplementary Table 2). Our find- evaluate the reactant movements on the nanoscale in HvExoI via ings suggest that although the mutations of Trp286 to Ala and His led classical molecular dynamics (cMD), Gaussian accelerated MD (GaMD) to folded but inactive enzymes, enzyme activity retention is more simulations and GPathFinder calculations. These findings combined sensitive to alterations of Trp286 than Trp434, likely due to the with modelling of the conformational behaviour of reactants via cMD proximity of Trp286 to Asp285, the catalytic nucleophile (Supple- and quantum mechanics/molecular mechanics (QM/MM) MD mentary Table 2). Nature Communications | (2022)13:5577 3 Article https://doi.org/10.1038/s41467-022-33180-5 a WT G2SG-OMe b WT 4NP-G3SG Glu287 Asp285 Trp286 Tyr253 Trp286 His207 4NP-G3SG G2SG-OMe PEG Lys206 Glu491 Gol Asp95 Gly57 Gly57 Arg158 Trp434 Glu491 Trp434 -1 subsite +1 subsite +2 subsite -1 subsite +1 subsite +2 subsite c WT G4SG4OG4SG4OG4SG d WT G6SG-OMe Arg291 Asp285 Trp286 Asp285 Tyr253 Trp286 His207 G6SG-OMe His207 G4SG4OG4S Lys206 Lys206 Asp95 Gly57 Asp95 Gly57 Glu491 Trp434 Glu491 Trp434 Arg158 Arg158 -1 subsite +1 subsite +2 subsite -1 subsite +1 subsite Fig. 2 | WT HvExoI in complex with thio-saccharide analogues. a Thio- as insets rotated by around −90 via x-axes relative to main panels) contoured at 1.0 o sophoroside (G2SG-OMe; pink); b thio-laminaribioside (4NP-G3SG; magenta); σ levels. Glycerol (Gol) or polyethylene glycol (PEG) and water molecules are shown c thio-cellotrioside (G4SG4OGS; cyan); and d thio-gentiobioside (G6SG-OMe; in cpk orange sticks and red spheres, respectively. Distances at separations within green) moieties are bound across the −1 to +2 subsites. Panels a to d contain derived 3.3 Å are shown in dashed lines. Ligand designations above structural images difference 2 m|Fo|-D|Fc| electron density maps (blue) of thio-analogues (also shown indicate thio-analogues that were perfused in crystals. Steady-state kinetic constants of WT, and Trp286 and Trp434 Inhibition constants Ki of WT, and Trp286 and Trp434 mutants mutants To further examine binding modes of β-D-glucosides, we used the thio- Consistent with our predictions, catalytic rate kcat (first-order rate) analogues with (1,2)-, (1,3)-, (1,4)- and (1,6)-linkages to quantify appar- constants and in some instances, Michaelis-Menten KM values of single ent inhibition constants (Ki) of the HvExoI:inhibitor complexes W286A, W286H, W286F, W434A, W434H, W434F, and W434Y, and (Supplementary Table 4). In WT HvExoI, the (1,3)-β-D-linked double W286F/W434F and W286F/W434A mutants with poly- and (4NP-G3SG), (1,4)-β-D-linked (G4SG-OMe), and (1,6)-β-D-linked (G6SG- oligomeric substrates decreased compared to WT (Supplementary OMe) thio-analogues showed a competitive inhibition with Ki values in Table 3). This trend continued with the catalytic efficiency kcat/KM low to medium micromolar ranges, except G2SG-OMe, which was (second-order rate) constants, where these values corresponded to weakly inhibitory (Supplementary Table 4). The W434A mutant was those of relative hydrolytic rates (Supplementary Table 2). These sensitive to 4NP-G3SG and G4SG-OMe, and the data with 4NP-G3SG values confirmed the strict (1,3)-β-oligosaccharide specificity of correlated with substrate specificity and catalytic efficiency (Supple- W434A, with little change in the apparent strength of binding of mentary Tables 2 and 3). As expected, the (1,3)- and (1,4)-linked thio- laminaribiose, but by an about 21-fold decrease of kcat (compared to inhibitors were less potent with tested mutants, while the most effec- WT), consequently yielding a low kcat/KM value. A similar profile was tive thio-inhibitor was G6SG-OMe, except that it did not bind to noted for W286F/W434A, although in this case, the kinetic parameters W434A, and compared to WT, it showed 8- and 6-fold higher Ki con- decreased even more, pointing to a negative effect of the non- stant values with W434H and W434F, respectively (Supplementary conservative (radical) W434A mutation on substrate binding and Table 4). These data pointed out that the replacements of Trp434 with hydrolysis (Supplementary Table 3). Notable was the behaviour of Ala, His, or Phe residues in the +1 subsite, disturbed the optimal binding W286Y that showed increased catalytic efficiency values with (1,3;1,6)- of thio-analogues, due to the removal or weakening stacking interac- β-D-glucan (laminarin), laminaribiose (G3OG), and 4NP-Glc, while these tions crucial for the correct positioning of inhibitors in the active site. values in W434Y were considerably lower (Supplementary Table 3). In W286F/W434F, the catalytic efficiency values were also drastically WT and mutant HvExoI:thio-β-D-glucoside complexes reveal the reduced compared to WT, whereas the W434H and W434F mutants details of binding at the atomic level showed higher catalytic efficiency than the double mutants, with We observed that the thio-analogue Ki constants (except G2SG-OMe) W434F being more active except with laminaribiose. These data signal measured with the Trp286 and Trp434 mutants, correlated with their that although the mutations of Trp434 affected the catalytic properties capacity to hydrolyse (1,3)-, (1,4)- and (1,6)-linked β-D-glucosides. To of HvExoI more substantially compared to Trp286, both Trp residues decipher their precise binding modes at the atomic levels, we pursued are required for binding the isomeric β-D-glucosides to achieve pro- structural studies (Figs. 2–4; Supplementary Figs. 2–4; Supplementary ductive binding modes required for high catalytic efficiency. Tables 5 and 6). Nature Communications | (2022)13:5577 4 Article https://doi.org/10.1038/s41467-022-33180-5 a W434A G2SG-OMe b W434A 4NP-G3SG3OG Arg291 Arg291 PEG Asp285 Tyr253 Trp286 Trp286 Tyr253 His207 G2SG-OMe His207 G3SG Lys206 Lys206 PEG Asp95 Gly57 Glu491 Asp95 Gly57 Glu491 Arg158 W434A W434A Arg158 -1 subsite +1 subsite -1 subsite +1 subsite c W434A G4SG4OG4SG4OG4SG d W434A G6SG-OMe Trp286 Asp285 Trp286 Asp285 Tyr253 G6SG-OMe His207 G4SG4OG4S His207 Lys206 Lys206 Gly57 Glu491 Asp95 Glu491 Gly57 Asp95 Lys493 W434A Arg158 Arg158 W434A -1 subsite +1 subsite +2 subsite -1 subsite +1 subsite Fig. 3 | Thio-saccharide analogues are bound in the active site of the HvExoI thio-analogues (also shown as insets rotated by around −90o via x-axes relative to W434A mutant. a Thio-sophoroside (G2SG-OMe, pink); b thio-laminaribioside main panels) contoured at 1.0 σ levels. Polyethylene glycol (PEG) and water mole- (4NP-G3SG, magenta); c thio-cellotrioside (G4SG4OGS; cyan); and d thio- cules are shown in cpk orange sticks and red spheres, respectively. Distances at gentiobioside (G6SG-OMe, green) moieties are bound across the −1 to +2 subsites. separations within 3.3 Å are shown in dashed lines. Ligand designations above Panels a–d contain derived difference 2 m|Fo|-D|Fc| electron density maps (blue) of structural images indicate thio-analogues that were perfused in crystals. G2SG-OMe binding. The G2SG-OMe thio-analogue bound weakly to native and WT recombinant HvExoI with the respective Ki of 2 located at the +1 to +2 subsites, with the 4NP group clearly discernible ×10−3 M and 4 ×10−3 M (Supplementary Table 4), reflecting its ability at the +2 subsite (Fig. 2b). However, in W286A, and W286Y, 4NP-G3SG to attach at the +1 and +2 subsites instead of the −1 and +1 subsites was not observed, and in W286F, the reducing-end Glc moiety of 4NP- (Fig. 2a)22, which was most likely due to the rigidity of this thio- G3SG at the +1 subsite remained disordered (Supplementary Fig. 2d). ligand. On the other hand, in W286F the electron density map of This emphasised that the Trp286 mutations to Ala, Tyr, and Phe dis- G2SG-OMe was defined at the −1 and +1 subsites (Supplementary turbed the binding of (1,3)-linked thio-β-D-glucosides and that this Fig. 2c), where this thio-ligand formed H-bonds with Asp95, Lys206, residue was necessary for inhibitor binding. His207, Arg158, Tyr253, Asp285, Arg291, and Glu491, and stacking interactions with W286F and Trp434, whereas G2SG-OMe did not G4SG-OMe and G4SG4OG4SG4OG4SG binding. G4SG-OMe bound bind to W286A and W286Y. The dispositions of G2SG-OMe at the −1 to native and WT recombinant HvExoI with the Ki values of 1.0 ×10−3 M and +1 subsites in the W434A, W434F, W434H, and W434Y mutants (Supplementary Table 4). This compares favourably to the hex- (Figs. 3a, 4a; Supplementary Figs. 3a and 4a) revealed similar asaccharide G4SG4OG4SG4OG4SG, where we previously observed the structural features of the ligand binding to those seen in W286F Ki constant of 0.6 ×10−3 M with native HvExoI2. These binding con- (Supplementary Fig. 2c). stants suggested that the (1,4)-linked thio-analogues should locate in the enzyme’s active sites. This occurred with WT, and the W434A, W434F, W434H, and W434Y mutants, while no binding was observed 4NP-G3SG and 4NP-G3SG3OG binding. 4NP-G3SG bound tightly to in the W286A, W284F, and W286Y mutants. The well-defined electron native and WT recombinant HvExoI with the Ki constant of 0.7 ×10−3 M density maps in W434A (Fig. 3c) and W434F (Supplementary Fig. 3c) (Supplementary Table 4); the Ki value of 0.2 ×10−3 M was previously revealed three saccharide moieties at the −1 to +2 subsites with the observed with 4NP-G3SG3OG in WT HvExoI3. These relatively tight reducing thio-saccharide moiety protruding beyond the environment binding constants suggested that these thio-analogues should be seen of the active site pocket at the +2 subsite. In W434A, we observed that in the active sites, as we observed with WT, and the W286F, W434A, the tri-thio-cellobioside moiety was somewhat misplaced (compared W434F, W434H, and W434Y mutants. The well-defined electron den- to WT) and that the W434A substitution unexpectedly allowed Lys493 sity maps of the (1,3)-linked thio-ligands illustrated that the saccharide to be rotated outward and become engaged in the binding of this thio- moieties were attached through the networks of H-bonds formed ligand; here, 11 residues participated in forming H-bonds with tri-thio- between nine to ten residues in each case (Figs. 2b, 3d, 4b; Supple- cellobioside moiety (Fig. 3c). To generate the W434H and W434Y mentary Figs. 2d, 3b, 4b). The 4NP-G3SG thio-analogue in WT was complexes, we perfused the G4SG-OMe ligand, consequently, only two Nature Communications | (2022)13:5577 5 Article https://doi.org/10.1038/s41467-022-33180-5 a W434H G2SG-OMe b W434H 4NP-G3SG Arg291 Arg291 Asp285 Asp285 Tyr253 Trp286 Trp286 Tyr253 His207 His207 G3SG G2SG-OMe Lys206 PEG Lys206 Gol Gly57 Gly57 Asp95 Asp95 Glu491 Glu491 W434H W434H Arg158 Arg158 -1 subsite +1 subsite -1 subsite +1 subsite c W434H G4SG-OMe d W434H G6SG-OMe Arg291 Asp285 Trp286 Tyr253 Asp285 Trp286 Tyr253 His207 G6SG-OMe G4SG Lys206 His207 Lys206 Asp95 Gly57 Gly57 Asp95 Gol Glu491 W434H Glu491 W434H Arg158 Arg158 -1 subsite +1 subsite -1 subsite +1 subsite Fig. 4 | Thio-saccharide analogues are bound in the active site of the HvExoI analogues (also shown as insets rotated by around −90o via x-axes relative to main W434H mutant. a Thio-sophoroside (G2SG-OMe, pink); b thio-laminaribioside panels) contoured at 1.0 σ levels. Polyethylene glycol (PEG) and water molecules are (4NP-G3SG, magenta); c thio-cellobioside (G4SG; cyan); and d thio-gentiobioside shown in cpk orange sticks and red spheres, respectively. Distances at separations (G6SG-OMe, green) moieties are bound across the −1 to +2 subsites. Panels a to within 3.3 Å are shown in dashed lines. Ligand designations above structural images d contain derived difference 2 m|Fo|-D|Fc| electron density maps (blue) of thio- indicate thio-analogues that were perfused in crystals. thio-saccharide moieties were observed in active sites (Fig. 4c; Sup- W286Y (Supplementary Fig. 2b) but that these maps were well-formed plementary Fig. 4c). But with W286H, and W286Y, we could not see the in W286A and W286F (Supplementary Fig. 2a, e), suggesting that in the binding of G4SG-OMe, which correlated with their weaker Ki values W286Y mutant the process of ligand-binding in (or near) the active site compared to WT (Supplementary Table 4), and suggested that the may have been perturbed. These findings reinforce the conclusion that conservative Trp286 mutation to Phe and Tyr affected the binding of G6SG-OMe binds in HvExoI regardless of mutations (particularly into (1,4)-linked thio-β-D-glucosides. aromatic residues) of Trp286 or Trp434. Although these mutations impact the precise orientation of the glucopyranoside moieties, as G6SG-OMe binding. G6SG-OMe bound tightly to native and WT observed in W434A, they permit G6SG-OMe binding. recombinant HvExoI with the respective Ki values of 0.1 ×10−3 M and 0.2 ×10−3 M (Supplementary Table 4), hinting that this thio-analogue Conformational behaviour of S- and O-linked β-D-glucosides should be observed in the active sites of WT and mutants (Figs. 2d, 3d, To determine the conformational landscape of β-D-glucosides bound 4d; Supplementary Figs. 2a, b, e, 3d, 4d). This is what we indeed in HvExoI, we analysed enzyme-thio-ligand-binding modes in the observed in WT, and the W286A, W286Y, W286F, W434A, W434F, crystal structures before and after cMD simulations. As described W434H, and W434Y mutants, where well-defined electron density above, the β-D-glucopyranose rings of thio-saccharide analogues at the maps of G6SG-OMe and tight H-bond networks that were formed −1 subsites in the WT crystal structure complexes were engaged in through eight to ten residues (Asp95, Lys206, His207, Arg158, Tyr253, H-bond networks and observed in the 4C1 conformation for G2SG-OMe Asp285, Arg291, Glu491, Trp286, and Trp434). In W434A we noticed and G6SG-OMe, the 4H3 conformation for 4NP-G3SG, and the 4E/4H3 approximately 120o rotation (through the x-axis) of the reducing-end conformations for G4SG4OG (Supplementary Figs. 5 and 6; Supple- moiety of G6SG-OMe compared to its disposition in WT (or in-solution mentary Table 7), while the β-D-glucopyranose moieties of thio-ligands structure of G6OG or gentiobiose), due to the presence of a rotatable at the +1 subsites adopted the 4C1 conformation. To analyse the C1-O-C6 glycosidic bond of G6SG-OMe; this did not occur in WT or dynamic behaviour of β-D-glucosides, a total of eight models were other mutants. To find out, if the Glc molecule alone adopts similar or generated for G2(S or O)G, G3(S or O)G, G4(S or O)G, and G6(S or O)G different dispositions at the −1 and +1 subsites, we perfused Glc into systems, wherein we reverted mutations (if present) in HvExoI, WT, W434A, and W286Y, and compared their positions, which at the removed ions and other molecules from thio-analogues, and for the −1 subsite were identical in all instances, but at the +1 subsite of W434A O-linked systems, replaced S-atoms with O-atoms. While cMD simula- the Glc molecule was disordered (data not shown). We noted that the tions revealed that β-D-glucopyranose rings in thio-analogues adopted quality of the G6SG-OMe electron density map slightly declined in predominantly the 4C1 conformations in both subsites, alternate Nature Communications | (2022)13:5577 6 Article https://doi.org/10.1038/s41467-022-33180-5 a G2OG b G3OG -1 +1 -1 +1 Asp285 Asp285 C2 C3 C1 E  C1 E  Glu491 Glu491 c G4OG4OG d G6OG -1 +1 +2 -1 +1 Asp285 Asp285  C1 E C4 C1 E  C6 Glu491 Glu491 Fig. 5 | Conformational behaviour of O-linked oligosaccharides bound in the 200 ns in Mercator projections are indicated). Complexes of WT with G3SG, active site of HvExoI calculated by cMD simulations plotted as a function of θ G4SG4OG, G6SG, and W434A in complex with G2SG were used to predict binding and φ puckering coordinates. a Sophorose (G2OG; pink); b laminaribiose (G3OG; modes of O-linked oligosaccharides. The β-D-glucopyranose moieties at the −1 magenta); c cellotriose (G4OG4OG; cyan); and d gentiobiose (G6OG; green) ligands to +2 subsites, and dispositions of Asp285 and Glu491 catalytic residues are are bound at the −1 to +2 subsites (colour gradients from purple at 0 ns to yellow at also shown. conformations also participated, such as 4H3/4E/4H5 or 1S3/1,4B, simulations of this complex showed that the Glc moieties of G4OG3OG depending on a type of the glycosidic linkage (Fig. 5; Supplementary adopted the stable 4C1 conformations at the +1 and +2 subsites (Sup- Fig. 6; Supplementary Table 7). Additionally, we reconstructed plementary Fig. 10). However, Glc at the −1 subsite adopted a state, that O-derivatives in the active sites based on the thio-analogue complexes was nearly identical to the Glc moiety in G3OG, and not entirely a and subjected these to cMD (Fig. 5; Supplementary Fig. 7; Supple- perfect chair, thus we defined it as a pseudo-chair (Supplementary mentary Table 7). Compared to the analyses of the thio-analogue Fig. 10b). Again, this is consistent with the conformational FEL of complexes discussed above, these data indicated that the β-D-gluco- G4OG3OG at −1 subsite, computed by QM/MM MD metadynamics pyranose rings of G2OG, G3OG, G4OG4OG, and G6OG adopted mainly (Fig. 6a; Supplementary Fig. 11), which shows a wide free energy 4 C1 conformations in both subsites, although other states contributed minimum that extends from 4C1 towards 4H3. to binding in the −1 subsite (32% in the 4H5/4E region for G3OG, and Additionally, we investigated the binding of the G3OG product respective 5%, 10% and 27% in 4H3, 4E, 4H5 regions for G2OG, G6OG, and originating from the G4OG3OG substrate in WT and the W434A G4OG4OG). In the +1 subsite, the G3OG β-D-glucopyranose ring was mutant, after the non-reducing Glc moiety was cleaved off, using cMD located at respective 66% and 28% in the 4C1 and 2SO regions, while the simulations (Figs. 6b, c). The cluster analysis of the trajectory of the β-D-glucopyranose rings of G2OG, G6OG, and G4OG4OG occupied the G3OG product in WT indicated that the Glc moieties adopted at least 4 C1 regions in more than 90% of occurrences (Supplementary Table 7). two poses: ~96% of the populations located in the +1 and +2 subsites, Notably, the Mercator projection of G3OG (Supplementary Fig. 8) and ~4% to the intermediary poses at the −1 and +1 subsites (Fig. 6b). concurred with its conformational Free Energy Landscape (FEL) map Conversely, the cluster analysis of the G3OG product in W434A indi- (Supplementary Fig. 9), computed by QM/MM MD metadynamics. As cated a dramatic shift in poses to ~15% of the populations located in the recently found for other exo-acting GH enzymes35,36, these findings +1 and +2 subsites, and ~74% to the intermediary (midway) poses at the predicted that HvExoI likely exhibits a two-fold hydrolytic itinerary −1 and +1 subsites (Fig. 6c). These observations imply that Trp434 plays with a positionally isomeric disaccharide such as G3OG, where the a key role in the sliding of the (1,3)-β-D-linked hydrolytic products hydrolytic reaction starts either from a distorted 4H3/4E/4H5 region or through the active site during catalysis and that Trp434 prevents the from an undistorted 4C1 β-D-glucopyranosyl ring, which would evolve sliding of mixed-linkage substrates. towards a transition state with a 4H3/4E conformation in both cases. Binding and displacements of laminaritriose (G3OG3OG) and Binding and hydrolysis of the G4OG3OG natural substrate were laminarihexaose (G3OG3OG3OG3OG3OG) hydrolytic products explored via cMD and QM/MM MD metadynamics simulations These studies were performed in the WT, W434H, and W434A enzymes To understand the binding of one of the natural substrates of HvExoI, with active sites accommodating hydrolytic products of laminaritriose, we selected the kinked G4OG3OG oligosaccharide (Fig. 6; Supple- i.e. Glc in the −1 subsite and laminaribiose (G3OG) in the +1 and mentary Figs. 10 and 11), which is the key hydrolytic product origi- +2 subsites, or laminarihexaose, i.e. Glc at the −1 subsite and laminar- nating from (1,4;1,3)-β-D-glucans; these polymers fulfil the role of ipentaose (G3OG3OG3OG3OG) in the +1 and +2 subsites, and exposed structural cell wall components in cereals, such as barley5. The com- to bulk solvent. The selection of these mutants was based on relative plex with G4OG3OG was generated based on the WT HvExoI:G4SG4OG hydrolytic rates of poly- and oligosaccharide substrates, where these complex (Fig. 2c), upon replacing G4 with G3. The analysis of cMD mutants either retain (W434H) or lose (W434A) their broad substrate Nature Communications | (2022)13:5577 7 Article https://doi.org/10.1038/s41467-022-33180-5 a Conformational FEL map of Glc at -1 subsite of G4OG3OG in WT Glc at -1 subsite -1 +1 +2 Asp285 Trp286 C1 E  E  C4 C1 C3 Glu491 Trp434 b Behaviour of the G3OG product after hydrolysis from G4OG3OG in WT -1 +1 +2 Glc at -1/+1 subsite Asp285 Trp286 C1 E  C3 95.5% Glu491 Trp434 -1 +1 +2 Asp285 Trp286 E 4.1% C1 C3 Glu491 Trp434 c Behaviour of the G3OG product after hydrolysis from G4OG3OG in W434A -1 +1 +2 Glc at -1/+1 subsite Asp285 Trp286 C3 14.6% C1  E Glu491 W434A -1 +1 +2 Asp285 Trp286 73.7% C1 E C3 Glu491 W434A Fig. 6 | Binding of the G4OG3OG natural substrate in the HvExoI active site (Mercator projection) of the G3OG product after the −1 Glc moiety is hydrolysed calculated via cMD and QM/MM MD metadynamics simulations, plotted as a from G4OG3OG and removed from the −1 subsite in W434A. Colour gradients from function of θ and φ puckering coordinates. a Conformational FEL map of the Glc purple at 0 ns to yellow at 150 ns in Mercator projections are shown. In moiety at the −1 subsite in WT; Glc moieties of G4OG3OG in the subsites −1 to +2 are b, c, respective cluster analyses of trajectories (populations in %) indicate two indicated along with the key residues; b conformational behaviour (Mercator modes of G3OG binding in WT and W434A, where each saccharide binds in the +1 projection) of the G3OG product after the −1 Glc moiety is hydrolysed from and +2 subsites or Glc moieties bind in intermediary positions between those G4OG3OG and removed from the −1 subsite in WT; c conformational behaviour subsites and the −1 subsite. specificity (Supplementary Table 2). To assess the effect of saccharides and laminarihexaose hydrolytic products showed contrasting on HvExoI behaviour, cMD simulations were also performed with the behaviour compared to WT (Supplementary Fig. 12). cMD simu- WT apo-form, lacking bound Glc or β-D-glucosides. lations of the mutants with the Glc and the G3OG products (run The cMD and GaMD simulations with laminaritriose and for 1000 ns) revealed that the protein backbones remained stable laminarihexaose hydrolytic products indicated that the WT (respective RMSD values: 1.7 Å for W434H, and 1.25 Å for W434A) structure was stable (Supplementary Fig. 12), as indicated by and that Glc was attached at the −1 subsite, but, after around backbone Root-Mean-Square-Deviation (RMSD) values, although 100 ns in W434H and 200 ns in W434A, G3OG moved to bulk high mobilities in some loops were detected. Conversely, cMD solvent, indicating that the mutants were ineffective in retaining and GaMD simulations of W434H and W434A with laminaritriose the G3OG product (Supplementary Fig. 13). Such behaviour was Nature Communications | (2022)13:5577 8 Article https://doi.org/10.1038/s41467-022-33180-5 not observed with laminarihexaose hydrolytic products, within For the W434A mutant (Glc at the −1 subsite and laminaripentaose the time-frames of cMD simulations (600 ns), whereby these in the +1 and +2 subsites, and bulk solvent), GPathFinder identified a products remained attached in WT and both mutants (Supple- different Glc exit path (path 2) (Figs. 7c, 8; Supplementary Movie 3). In mentary Fig. 13). These simulations indicate that the dynam- this mutant, Glc exited through the preformed active site opening after ics of ternary complexes of mutants presented variabilities the non-conservative mutation of W434 into Ala, while the Gly210- compared to WT, where these differences depend on the length Met220 loop adopted a position that was seen in WT. This was sup- of hydrolytic products bound, i.e. longer products were more ported by the higher flexibility of the region containing W434A, which stable. adopted an altered conformation resulting from the displacement of To resolve variability with the length of hydrolytic products in Trp430 that could otherwise block path 2. In this Glc egress path, the mutants, GaMD simulations (800 ns) with enhanced sampling were toll-like Arg158-Asp285-Glu491 barrier, Tyr253, and the Arg291-Glu220 conducted. We found that the protein scaffold of the W434A mutant salt bridge were not involved. was more flexible (RMSD value 2.5 Å) than that of WT or W434H (RMSD When laminaripentaose was removed from the +1 and +2 subsites values between 1.0 Å − 1.5 Å) (Supplementary Fig. 12), and that the (and bulk solvent) in WT, W434H, and W434A, we identified path 3 (we Gly210-Met230 loop in the vicinity of the catalytic site exhibited Root- termed it vertical) for Glc displacement (Fig. 8), where Glc utilised the Mean-Square-Fluctuation (RMSF) values of up to 7 Å in both mutants space vacated by laminaripentaose and exited via the Trp286/Trp434 compared to WT (RMSF value ≤ 2 Å). While Glc remained highly aromatic clamp, and through the +2 subsite into bulk solvent. This restrained in the −1 subsite in WT and W434H, it was mobile in W434A, path could only be available when no products or substrates are hinting that the overall flexibility of the catalytic site was little affected bound. This scenario for Glc egress is, however, only a hypothetical in W434H, but strongly affected in W434A. During GaMD simulations, option, also challenging our previous study22, where in native HvExoI, the G3OG and laminaripentaose products remained bound in WT, but two Glc molecules were bound at the −1 and +1 subsites, effectively G3OG vacated the catalytic site in both mutants, in contrast to lami- blocking path 3. naripentaose, which initially remained bound (~25% and ~50% of In summary, GPathFinder simulations revealed that path 1 is the simulation time for W434A and W434H, respectively), but eventually, it preferred route in WT and W343H with bound (1,3)-linked hydrolytic traversed to bulk solvent. These findings suggested that the replace- products, and is independent of a substrate length, contrary to W434A, ments of Trp434 by His and specifically by Ala led to a loss of G3OG and which utilises path 2, where the essential structural attributes of laminaripentaose retention, and thus, a potential unbinding of these substrate-product assisted processivity (toll-like barrier, salt bridge, products could occur before the next catalytic cycle begins. This also loop motions, lateral cavity) do not participate. implies that in W434H and W434A, but not in WT, the (1,3)-linked products would need to re-bind and cannot slide through the active Discussion site uninterrupted. In nature, processive (non-dissociative) catalysis is omnipresent in enzymes that synthesise, modify, and degrade biopolymers such as Glc displacement route (egress) or substrate-product assisted nucleic acids, proteins, and polysaccharides, and has been described in processivity in WT nucleases, peptidases, proteases, proteasomes, and endo- and exo- The modelling of the Glc displacement from the −1 subsite performed acting polysaccharide hydrolases22–25,28–31,38. Processive polysaccharide with GPathFinder37, commenced from the most populated cluster of hydrolases have certain common characteristics that stem from the structures based on cMD, to further explore dynamic traits of Glc architectures of their catalytic sites that dictate spatial arrangements unbinding. These calculations were executed in two arrangements governing substrate specificity and binding (including directionality) with (i) Glc at the −1 subsite and laminaripentaose in the +1 and and the dispositions of catalytic residues to ensure a high efficacy of +2 subsites and bulk solvent (Figs. 7, 8; Supplementary Movies 1–3), or catalysis. Processive endo- and exo-acting cellulases23,28,39,40 and exo- Glc at the −1 subsite and G3OG in the +1 and +2 subsites (Supple- acting chitinases24,25,29,31 are known to be driven by processive hydro- mentary Fig. 14); and (ii) Glc in the −1 subsite and the di- and penta- lysis or substrate-assisted catalysis, and the GH3 exo-acting β-D-glucan saccharide products removed to explore the Glc egress route without glycohydrolases with a pocket-shaped active site (such as HvExoI) by products being bound. substrate-product assisted processive catalysis22. The difference Initial simulations performed with WT (Glc at the −1 subsite and between these two types is in the role of a hydrolytic product, that in laminaripentaose in the +1 and +2 subsites and bulk solvent) yielded the latter group evokes the formation of a lateral cavity, which serves as the Glc egress path (path 1) (Figs. 7a, 8; Supplementary Movie 1), a conduit for Glc egress to allow for the next catalytic cycle using the which was facilitated by the re-orientation of a toll-like Arg158- retaining catalytic mechanism (Fig. 9). This type of catalysis enables Asp285-Glu491 barrier (including the side-chain of Asp285) and the multiple events that proceed without an enzyme losing contact with motions of the loop carrying Tyr253. Next, Glc migration induced a oligo- or polymeric substrates, contrary to distributive catalysis, where rotation of Tyr253 that opened the lateral cavity due to the dis- the catalyst follows a conventional association-conversion-dissociation ruption of the Arg291-Glu220 salt bridge – these re-orientations and pattern. motions collectively opened the lateral cavity for Glc to exit. During The crucial role of substrate binding in HvExoI implies a pair of these simulations, CH–π interactions between the non-reducing- the aromatic Trp286 and Trp434 residues that form a clamp at the end Glc moiety of laminaripentaose and the Trp286/Trp434 aro- +1 subsite, which is positioned next to a pocket-shaped site at the matic clamp remained stationary, holding the Glc moiety in the −1 subsite, housing the Asp285 and Glu491 catalysts1–3. These resi- +1 subsite. dues, typically around 5.5-6.5 Å apart locked in a rigid configuration GPathFinder simulations with W434H (Glc at the −1 subsite and via a dense H-bond network, are responsible for the hydrolysis of laminaripentaose in the +1 and +2 subsites and bulk solvent), revealed oligo- and polymeric β-D-glucosides to release Glc from non-reducing that the Glc egress path was identical to that of WT (path 1) (Figs. 7b, 8; ends of substrates (Figs. 1a, b, and 9). To understand the mechanism Supplementary Movie 2). Here, a toll-like Arg158-Asp285-Glu491 bar- that drives the unidirectional motion of the Glc product in HvExoI rier also re-oriented, Asp285 changed its rotameric state, and as Glc during substrate-product assisted processive catalysis, and specifi- advanced through the lateral cavity, the side-chain of Tyr253 re- cally the roles of Trp286 and Trp434, we analysed substrate specifi- orientated, and the Arg291–Glu220 salt bridge was broken – these city in the context of elementary steps of Glc movement. We created events enlarged the lateral cavity volume and facilitated Glc egress HvExoI mutants replacing Trp286 and Trp434, with Ala, Phe, Tyr, and from the −1 subsite to bulk solvent. His residues, and compared their substrate specificity, steady-state Nature Communications | (2022)13:5577 9 Article https://doi.org/10.1038/s41467-022-33180-5 a Glc displacement route in WT with laminaripentaose bound in +1 and +2 subsites and bulk solvent Step 1 Step 2 Step 3 Step 4 Glc displacement: 0.4 Å Glc displacement: 4.8 Å Glc displacement: 9.0 Å G3OG3OG3OG3OG Trp286 Asp285 Arg291 Glc Glu220 Trp434 Glu491 Thr496 Arg158 -1/+1 +2 b Glc displacement route in W434H with laminaripentaose bound in +1 and +2 subsites and bulk solvent Step 1 Step 2 Step 3 Step 4 Glc displacement: 6.9 Å Glc displacement: 9.1 Å Glc displacement: 9.9 Å G3OG3OG3OG3OG Arg291 Asp285 Trp286 Glu220 Glc Glu491 W434H Arg158 +2 +1 -1 c Glc displacement route in W434A with laminaripentaose bound in +1 and +2 subsites and bulk solvent Step 1 Step 2 Step 3 Step 4 Glc displacement: 1.6 Å Glc displacement: 6.0 Å Glc displacement: 9.4 Å G3OG3OG3OG3OG Arg291 Trp286 Asp285 Glc Glu220 W434A Glu491 Arg158 +2 +1 -1 Fig. 7 | Glc displacement routes in WT, and the W434H, and W434A mutants of the participation of the Arg158-Asp285-Glu491 toll-like barrier (triangles in dashed HvExoI with Glc and laminaripentaose, calculated through cMD and GPath- lines; distances indicated in Å). Surface morphologies of structures are coloured by Finder. a WT; b W434H; c W434A – Four steps along Glc displacement routes based electrostatic potentials: white, neutral; blue, +5 kT·e−1; red, −5 kT·e−1. Separations in on converged structures in complex with Glc (−1 subsite; carbons in cpk yellow Å between the positions of C1 carbons of Glc molecules, as they move from the lines, and spheres) and laminaripentaose (+1 and +2 subsites and bulk solvent; −1 subsites in initial and final structures are indicated. The lateral cavity (a, b; carbons: magenta sticks) obtained by docking. Selected residues (carbons: atomic indicated in black dotted ellipsoids) forms transiently, is partly exposed to bulk sticks), and positions of Glc, and laminaripentaose are indicated. Glc in the solvent, and facilitates Glc displacement. In panel c, Glc is displaced to bulk solvent −1 subsite is separated from the lateral cavity that is evolving during Glc egress by through the opening formed by the W434A mutation. kinetics, and inhibition parameters relative to WT. We conducted Dictyostelium, and have emerged with the evolution of land plants – atomic resolution (1.5 Å − 2.5 Å) X-ray crystallographic analyses of these enzymes are present in liverworts and mosses (Embryophyte), 25 structures of WT and mutant forms in complex with positional ferns and conifers (Tracheophyte), a basal Angiosperm Amborella isomers of thio-analogues, revealing the details of β-D-glucoside trichopoda, and dicotyledons (Eudicots) and monocotyledons binding. We examined the conformational behaviour of thio- (including grasses) (Fig. 1c). The β-D-glucan glucohydrolases with the analogues and O-linked β-D-glucosides, and explored the binding Trp286/Trp434 signature are also richly represented in bacteria, e.g. and hydrolysis of the G4OG3OG natural substrate via cMD and QM/ the sludge Chloroflexi bacterium fermenting carbohydrates41, and a MM MD metadynamics simulations. To define the atomistic details plant-pathogenic Xanthomonas (Fig. 1c), similarly to other bacteria of and dynamics of reactant movements, we examined the unidirec- the GH3 family4 (Bacteroides ovatus, Niveispirillum irakense, Pseu- tional motion of Glc and its egress, and substrate sliding through the doalteromonas sp. BB1), although β-D-glucan glucohydrolases from catalytic site or binding and displacements of laminaritriose and Rhodothermus marinus and Saccharopolyspora erythraea have the laminarihexaose hydrolytic products in the W434H and W434A Tyr286-Trp434 signature34. mutants through cMD, GaMD, and GPathFinder simulations. These During the investigations of physicochemical properties and data were corroborated with phylogenomic analyses of GH3 family substrate specificity in HvExoI and its mutants, we detected that the members, to trace the evolutionary origin of substrate-product impact of Trp286 and Trp434 mutations on the hydrolysis of 4NP-Glc assisted processivity, where β-D-glucan glucohydrolases are expec- under a variety of pH and temperature conditions was relatively small ted to retain an aromatic Trp286/Trp434 clamp. In the unrooted compared to WT, except for W286H with its pH optimum shifted to the phylogenetic tree (Fig. 1c; Supplementary Table 1), we observed the alkaline region (from 5.25 for WT to 6.7 for W286H) (Supplementary clustering of plant β-D-glucan glucohydrolases to subclades 1-3 with Fig. 1) due to a proximity effect to Asp285. However, the relative both Trp residues conserved or diversified. We detected that through hydrolytic rates and specificities of mutants with (1,2)-, (1,3)-, (1,4)- and the evolutionary history of plants, β-D-glucan glucohydrolases with (1,6)-linked poly-, and oligosaccharides, and aryl-glycosides, were the aromatic clamp are absent in red alga Porphyra and unicellular substantially altered. The W434A mutant lost its hydrolytic activities green algae Cyanobacterium and Chlamydomonas, and amoeba with polymers and oligomers, except (1,3)-linked oligosaccharides, Nature Communications | (2022)13:5577 10 Article https://doi.org/10.1038/s41467-022-33180-5 Laminaripentaose Laminaripentaose Trp286 Trp286 Glc Glc Path 2 Path 2 Path 3 Trp434 Path 1 Path 3 Trp434 Path 1 Fig. 8 | Glc egress paths 1-3, illustrated on the top of WT HvExoI with bound Glc with conforming coordinates for Glc in each frame. Glc egress paths were identi- (carbons in cpk yellow lines, and spheres) and laminaripentaose (cpk magenta fied by cMD simulations and GPathFinder. Surface morphology is coloured by sticks). Path 1 (WT and W434H), path 2 (W434A), and path 3 (vertical) are depicted electrostatic potentials: white, neutral; blue, +5 kT·e−1; red, −5 kT·e−1. Positions of in orange, cyan, and magenta spheres, respectively, and aligned next to each other. Trp286 and Trp434 delineating the +1 subsite, are indicated in grey sticks. The right Path 3 is found when no ligand is bound at the +1 and +2 subsites. Illustrated panel depicts the detailed view of paths 1-3, with HvExoI being rotated by around trajectories represent frames that correspond to consecutive stages of Glc egress −80o via the y-axis relative to the left panel. and acted as a strict (1,3)-β-D-specific β-D-glucosidase (Supplementary conformational flexibility of glucopyranose rings was present. This Table 2). In this mutant, no activity was observed with laminarin [(1,3)- was supported by our previous observations that D-glucono−1,5-lac- β-linked polymer with (1,6)-β-branching at the 3:2 ratio42] suggesting tone (transition state analogue) with inherently a half-chair con- that the binding of branched substrates may be influenced by steric formation was an effective inhibitor (Ki = 2.1 × 10−6 M)21. This is effects. These modified hydrolytic rates in the Trp286 and Trp434 acceptable because the active site must be able not only to recognise mutants were followed by apparent first- and second-order rate kinetic the ground state of saccharides in 4C1 conformation but also dis- (Supplementary Table 3) and apparent inhibition constants (Supple- torted conformations43,44. We concluded that the electron density mentary Table 4). Changes in these parameters suggested that Trp286 maps observed in the complexes could represent snapshots of time- and Trp434 play key roles in orienting positionally isomeric β-D-glu- space averaged conformations, that are formed at different time cosides towards the formation of productive HvExoI:substrate intervals on the substrate-binding trajectory from the ground to complexes. transition states. To address the atomic details underlying substrate binding and One significant observation stemmed from the comparative ana- hydrolysis, the crystalline WT and mutant HvExoI forms were per- lyses in terms of binding of (1,2)-, (1,3)-, (1,4)- and (1,6)-linked thio- fused with (1,2)-, (1,3)-, (1,4)- and (1,6)-linked thio-β-D-glucosides fol- glucosides by HvExoI (Figs. 2–4; Supplementary Figs. 2–4). In these WT lowed by X-ray crystallography analyses (Figs. 2–4; Supplementary complexes, the Glc moiety at the −1 subsite was oriented almost Figs. 2–4; Supplementary Tables 5 and 6) to assess the contributions identically, contrary to the Glc moiety at the +1 subsite that adopted of residues that interacted with glycone/aglycone moieties and thus various orientations. In G3SG the apolar β face of Glc moiety was their involvements in catalysis. Our interpretations of the electron exposed to Trp286, while in the HvExoI:G2SG,:G4SG, and:G6SG com- density maps of these complexes (25 in total) led us to postulate that plexes, the apolar thio-analogue sides mostly faced Trp434. It is known it was primarily the aromatic clamp, that binds variously linked thio- that in these stacking interactions, the saccharide surface encounters a glucosides at the +1 subsite, and mediates broad substrate specificity side-chain aromatic moiety when an apolar face of a saccharide in HvExoI. In the bound (1,2)- and (1,6)-linked thio-glucosides, we exposes three or more C–H groups close in space45. Concerning the observed 4C1 conformations for Glc moieties at the −1 subsite, or 4H3 substrate specificity from another angle, it is important to compare the or 4H3/4E respective conformations for the (1,3)- and (1,4)-linked thio- G6SG position in enzyme forms. While the Glc moieties of G6SG at the glucosides (Supplementary Table 7), suggesting that certain −1 subsite were oriented nearly identically in WT and mutant Nature Communications | (2022)13:5577 11 Article https://doi.org/10.1038/s41467-022-33180-5 a Retaining catalytic mechanism Glu491 O O O O HO O H ROH H O O O O OR H OH O O O O O O Asp285 b Reaction kinetics k1 k2 k3 k4 E+ S E~S E~P1~Glc E~P1 E~Glc -k1 -k2 Glc Glc c Substrate-product product assisted processivity 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: Dimer binds in Dimer is hydrolysed, After new dimer binds, Dimer progresses to productive mode. Glc at +1 subsite diffuses Glc remains catalytic site. Glc exits from away, Glc at -1 subsite remains bound in -1 subsite. -1 subsite via lateral cavity bound non-covalently. formed near +1 subsite. d Glc-product egress trajectories Step 1 Step 4 Step 1 Step 4 G3OG Trp286 substrate Trp286 Glc Glc product product G3OG substrate Trp434 W434A Fig. 9 | Mechanistic view of the catalytic function of HvExoI, a family GH3 squares) bound in the −1 and +1 subsites (step 1) is hydrolysed and the reducing end enzyme. a HvExoI uses the retaining catalytic mechanism, where nucleophile Glc (aglycon) diffuses away from +1 subsite, the non-reducing end Glc (cyan square) Asp285 (magenta) and acid/base catalyst Glu491 (green) play key roles in substrate remains non-covalently trapped. After a new dimer binds (step 3) and advances to hydrolysis1–3. b Reaction kinetics indicates the formation of Michaelis complexes the catalytic site, Glc modifies its binding patterns and egresses (large arrow) via a (E-S and E~S), and enzyme-product complexes before (E~P1~Glc) and following the transient lateral cavity (cylinder in dotted lines) formed near the Trp clamp first Glc egress (E~P1), and after hydrolysis of P1 (E~Glc)1, 2. S and P1 are trimeric and (step 4)22. Curved arrows indicate dimer entries into the catalytic site. d Glc dis- dimeric molecules, respectively. First-order rate constants are shown. Entrapped placement paths in WT (left panels), and the W434A mutant (right panels) with Glc Glc product from E~Glc is displaced after a new substrate is attached. The hydro- and the G3OG dimer, calculated through cMD and GPathFinder. Large arrows lytic cycle with a (1,3)-linked substrate is repeated until it is hydrolysed to Glc indicate Glc positions at the initial (step 1) and after the final (step 4 – following Glc (indicated by vertical and horizontal lines). c Simplified mechanism of Glc dis- departure from −1 subsites) stages of catalysis (this work). placement with a disaccharide22. After the disaccharide (empty blue and filled grey structures, the position of the Glc moiety at the +1 subsite varied – it (consisting of London dispersion forces)46. The geometries of these was identical in WT and W286A, while in W434A it was rotated by 120° interactions with Trp, Phe, Tyr, and His involved face-to-face geome- through the x-axis. This comparison underscored the pivotal role of tries, although the edge-face and offset stacked geometries also occur Trp434 in correctly orienting saccharide substrates engaged in stack- in protein complexes47. These stacking interactions are mediated by ing interactions, also known as dispersion-driven C-H-π forces Trp residues in 40%, trailed by Tyr, Phe, and His48. Nature Communications | (2022)13:5577 12 Article https://doi.org/10.1038/s41467-022-33180-5 To address the structural features of β-D-glucosides binding (or (W434H), or behaved as a strict (1,3)-β-D-glucosidase (W434A). Michaelis complexes) in HvExoI at the atomic level, we compared the Summarising these findings, we revealed that: (i) the stability and conformational behaviour of thio-analogues bound in crystalline time occupancy of the laminaribiose or laminaripentaose products in forms (Supplementary Figs. 5 and 6) and after cMD simulations (Sup- the subsites +1 and +2 were affected by W434H and W434A muta- plementary Fig. 8). These analyses revealed 4C1 conformations of the tions; (ii) in WT and W434H, the lateral Glc egress path was used glucopyranose moieties in the +1 subsite (Fig. 5; Supplementary Fig. 8; during substrate-product assisted processivity (path 1; Fig. 8); (iii) in Supplementary Table 7), and the pseudo-4C1 conformation for the W434A, Glc exited through the W434A-Trp286 clamp and a gap sugar moiety bound at the −1 subsite, although previous cMD analyses formed by the non-conservative W434A mutation, and without of a fungal GH3 enzyme did not find evidence of sugar distortion47. engaging substrate-product assisted processivity (path 2; Fig. 8); (iv) These findings supported the analyses of puckering conformations of in WT and both mutants, Glc egressed through the aromatic Trp286/ glucopyranose rings in Michaelis complexes of other exo-acting GH Trp434 clamp, when no oligosaccharides are present at this site (path enzymes2,35,36, and suggested that these exo-acting enzymes differ 3; Fig. 8), which was energetically most favoured; (v) substrate- from endo-acting ones, where substrate distortions were system- product assisted processivity can only be utilised, when the atically observed49–53. unbinding of the product at the +1 subsite is slower than Glc egress Additionally, via cMD simulations and QM/MM MD metady- through the lateral cavity (Figs. 7–9; Supplementary Fig. 14; Supple- namics, we investigated the trajectory of the G3OG product origi- mentary Movies 1–3). Finally, our unconstrained MD simulations of nating from the G4OG3OG substrate, where the Glc moieties in 4C1 HvExoI complexed with laminaripentaose, and Glc removed from the conformations and FEL maps shared properties (Fig. 6; Supplemen- −1 subsite revealed that the sliding of laminaripentaose proceeded tary Figs. 10 and 11). The cluster analysis of the G3OG product in WT spontaneously, even in W434H or W434A, meaning that all enzymes and W434A showed positions in markedly different proportions position Glc moieties in the active sites correctly due to the structural (Fig. 6b). This led us to conclude that the G3OG product could not complementarity. overcome a high energy barrier to achieve the productive binding Importantly, the combined experimental and computational mode, due to an unfavourable orientation of the Glc moiety in the approaches addressed the key role of the evolutionary conserved −1 subsite. Here, the CH2OH substituent faced the acid/base Glu491 Trp286/Trp434 clamp in the GH3 family, which secures Glc unbinding residue, while the anomeric carbon was positioned away from the and through C-H-π interactions guarantees the optimal binding of Asp285 nucleophile. A change in the ring puckering of the Glc moiety hydrolytic products while preventing their attachments by reducing- could eventually cause the anomeric carbon to approach Asp285, but terminals. These approaches also revealed another important func- such a possibility was eliminated based on the calculated FEL map tion, which is to assure that when 1,3-β-D-linked polymeric substrates (Fig. 6a). In W434A the G3OG product disposition allowed more bind to the active site and are hydrolysed with Glc egressing through flexibility, but G3OG still faced a high energy barrier to achieving a an ad-hoc formed lateral cavity, the hydrolytic products could slide productive binding mode. through the catalytic site and are rectified via substrate-product These findings have profound implications for the hydrolysis of assisted processivity to warrant a next hydrolytic event. This process is mixed-linkage β-D-gluco-oligosaccharides by HvExoI, where the accompanied by a cork-like motion between domains 1 and 222, which G4OG3OG natural substrate cannot slide through the active site from associates with an active site opened and closed states (Supplementary the +2 and +1 to the +1 and −1 subsites, although it could utilise the Movies 1–3). lateral cavity exit mechanism for Glc egress. The binding and Our studies provide a blueprint for emerging studies of reactant hydrolysis of the G4OG3OG natural substrate for HvExoI proceed in movements to advance the knowledge of catalysis in hydrolases and two steps: (i) after the 1st glycosidic linkage (G4OG) from the non- demonstrate that specificity and efficiency could be engineered via reducing end is cleaved off, and Glc egresses through the lateral residue alterations that dock substrates into active sites. Particularly, cavity; (ii) the G3OG product diffuses out of the active site and re- the work through our multiscale strategy contributes to the under- rebinds in the correct orientation to form a productive enzyme:- standing of catalysis in GH3 exo-hydrolases by defining struc- substrate complex. This distributive catalytic process can occur with tural trajectories and mechanisms for the breakdown of oligo- and trimeric or longer substrates that are (1,2)-, (1,4)-, (1,6)- or (1,4;1,3)- polysaccharides (Fig. 9). We suggest that our findings could be utilised linked, but not with (1,3)-linked substrates, which utilise both assets in enzyme design, as they are correlated with the statistical energy of processivity, i.e. Glc egress and substrate single-chain sliding inferred from a maximum entropy model of homologous sequences56. movement forward through the catalytic site without dissociation. Further, this knowledge can be exploited in the biomolecular design of Notably, in HvExoI, the processive character is paramount to the artificial molecular motors57,58, and more broadly in biotechnology and catalytic efficiency and is correlated with high kcat values of the (1,3)- drug design, where enzymes are sought as environmentally friendly linked substrates (Supplementary Table 3). Similar (i-ii) steps can be and cost-effective alternatives in bioremediation, and applications in used by G4OG4OG3G, which with G4OG3OG are the key low- medicine and pharmaceutics. molecular hydrolytic products generated from mixed-linkage (1,4;1,3)-β-D-glucans – polysaccharide components in cereal cell Methods walls5,6,54. Generation of HvExoI mutants, protein expression, and We alluded to the pivotal role of Trp434 in HvExoI that through purification stacking interactions participated in the orientations of substrates to Mutant HvExoI forms were prepared by site-directed mutagenesis positions amenable for hydrolysis through substrate-product assis- using the pPICZαBNH8 vector containing an optimised cDNA fragment ted processivity. To extend our knowledge of the role of Trp434 in encoding mature HvExoI as a template22,59,60. Trp286 and Trp434 were these processes, and why the evolution converged to Trp434, we singly replaced with alanine (Ala), histidine (His), phenylalanine (Phe), applied cMD, GaMD55, and GPathFinder37 simulations, where cMD and tyrosine (Tyr). Doubly mutated W286F/W434F and W286F/W434A simulations safeguard the stability of complexes in geometries close forms were also prepared to yield a total of ten mutants. Site-directed to those in crystal structures, GaMD provide an enhanced sampling mutagenesis was carried out as described in the QuikChange® Site- of conformational spaces, and GPathFinder explores dynamic traits Directed Mutagenesis Kit (Agilent Technologies, Santa Clara, CA, USA) of reactants binding or unbinding, via backbone motions and side- using primers specified in Supplementary Table 8. DNA constructs chain optimisations of discrete rotameric states. For these studies, were transformed in Pichia pastoris strain SMD 1168H to express and we selected two mutants, which retained broad substrate specificity purify proteins as described22,59,60. Nature Communications | (2022)13:5577 13 Article https://doi.org/10.1038/s41467-022-33180-5 Enzyme stability, substrate specificity, kinetic analyses, and domains (domain 1, residues 1-357; domain 2, residues 374-559) as two inhibition independent anisotropic domains with translation-libration-screw Dependency of activities of WT and mutant HvExoI on pH (McIlvaine (TLS) motion70. Electron density maps of ligands were defined in the buffers in the pH range of 3.5–8.5 at 30 °C) and temperature (100 mM active site regions at 3σ levels in the m|Fo|-|Fc| maps (m is the figure of sodium acetate buffer, pH 5.25, temperature range 0–80 °C, 15 min merit and D is an estimated coordinate error) and the water molecules incubations) were determined using 0.2% (w/v) 4NP-Glc21,59. Relative were located automatically with CCP4 ARP at levels higher than 2σ rates of hydrolysis of β-linked poly-, oligosaccharides, and aryl- levels. Water molecules were retained if they satisfied H-bond criteria glucosides were determined using near-homogenous HvExoI at 50 and if the 2 m|Fo|-D|Fc| electron density maps were confirmed. During pmoles. Apparent first-order rate parameters (KM and kcat) were model building and refinement, 5% of the data were flagged for cross- determined with (1,3;1,6)-β-D-glucan (laminarin), barley (1,4;1,3)-β-D- validation using Rfree statistics71,72. PROCHECK73 was used to check the glucan, laminaribiose, cellobiose, and 4NP-Glc with substrate con- stereo-chemical quality of models. Ramachandran plots of structures centrations ranging between 0.3-3 times KM values in duplicate as showed that 99.8% of the residues were found in the most favourable, described59,61. Apparent inhibition constants of enzyme:inhibitor additionally allowed, and generously allowed regions of the plot, with complexes (Ki) using the 4NP-Glc substrate were determined at six the well-defined density map for Ile432, the only residue located in concentrations with G2SG-OMe, 4NP-G3SG, G4SG-OMe, and G6SG- disallowed regions. Structural images were prepared in PyMOL Mole- OMe at 0.4-3 times the Ki values in duplicate as described21,22. Kinetic cular Graphics System v2.3.3 (Schrödinger LLC, Portland, OR, USA). parameters were calculated from the analyses of Michaelis-Menten The summary of 25 structures with thio-analogues and Glc and PDB plots with Grafit 7.0.362 and processed by a proportional weighted fit, accession numbers, data collection, and refinement statistics are using nonlinear regression63. The relative rates of hydrolysis, kinetic, reported in Supplementary Table 6. The (2 m|Fo|-D|Fc|) difference and inhibition data were rounded to one or two decimal points (Sup- electron density maps were calculated using Phenix67 with the mtz and plementary Tables 2–4). PDB files after CCP4 REFMAC5 refinements. The maps were contoured at the 1.0 σ level. Ligands bound in active sites (Figs. 2–4; Supple- Protein crystallisation and ligand perfusion mentary Figs. 2–4) were visualised in PyMol Molecular Graphics Sys- Near-homogenous WT and mutated HvExoI were crystallised by tem v2.3.3. macro- and cross-seeding with native HvExoI microcrystals22,60,64. Briefly, crystallisation was carried out at 4˚C in handing drops cMD simulations of Michaelis complexes derived from thio- consisting of 4-6 μL of HvExoI (4.5−12 mg/mL) and 4 µL of pre- ligand complexes cipitant solution A (100 mM HEPES-NaOH buffer, pH 7.0 containing A total of eight models were generated for the following systems: 2.3% (w/v) polyethylene glycol 400 and 1.6 M ammonium sulfate). W286F-G2(S or O)G in complex with G2SG-OMe; PDB 6JG7; WT-G3(S Drops were placed on 22 mm siliconised circular glass coverslips, or O)- in complex with 4NP-G3SG; PDB 6JG2; WT-G4(S or O)G in into which three to four 5−10 µm crystals of native HvExoI64 were complex with G4SG4OG4S; PDB 6JG1, and WT-G6(S or O)G in com- transferred by a cat whisker. Coverslips were suspended over 1 mL plex with G6SG-OMe; PDB 3WLP. The procedure consisted of (i) of mother liquor solutions (1.7 M ammonium sulfate in 46.9 mM reverting the mutation in the W286F crystal structure; (ii) removing HEPES-NaOH buffer, pH 7.0). Crystals of mutants grew within 10-28 SO42- anions, glycerol (GOL), polyethylene glycol (PEG), N-acetyl days at 4˚C to 100-400 µm sizes in the longest dimension. Thio- glucosamine (NAG), and aglycon moieties of thio-analogues; (iii) analogue and Glc perfusion at saturating concentrations (list and neutralising the systems with three Na+ cations; and (iv) adding a conditions summarised in Supplementary Table 5) was carried out cubic water box around proteins. For O-linked saccharide systems, an as described22. Briefly, crystals were perfused with thio-analogues or additional step was taken by replacing S-atoms in thio-analogues with Glc dissolved in precipitant solution A for 5-720 min. After various O-atoms. The tleap code, included in AMBER 1174 was chosen to perfusion times at 4±2 °C or ambient temperature, crystals were parametrise systems using the ff99SB force field75 for the protein, cryo-protected with 15% (v/v) glycerol in solution A and mounted on Glycam0676 for saccharides, and the TIP3P force field77 for water goniometers. molecules. The equilibration protocol consisted of a minimisation step (8 ×103 steps), progressive heating to 300 K (700 ps), NPT Crystal structure determination equilibration under a pressure of 1 atm to ensure physiological den- The X-ray diffraction data, except the WT:4NP-G3SG complex (Sup- sity (100 ps), and a final NVT equilibration (with a time step value of plementary Table 6), were collected at the MX1 beamline of the Aus- 2 ps using the SHAKE78 algorithm to constrain bonds involving tralian Synchrotron at 100 K (Oxford Instruments, UK) with a hydrogens) until the RMSD values of protein backbones stabilised. collimating mirror, double-crystal Si(111) monochromator on 03BM1 Default parameters were used for non-bonding interactions and dipole/bending magnet, and ADSC Quantum 210r Detector. X-ray data long-range electrostatics. Simulations at NVT were extended to of the WT HvExoI-4NP-G3SG complex were collected at 100 K (Oxford 200 ns each. Instruments, UK) on an in-house Rigaku MicroMax-007HF with a second-generation true microfocus rotating anode generator, oper- cMD simulations of the Michaelis complex with a natural sub- ating at 40 kV and 50 mA, Osmic VariMax optics for Cu, Kα, and the strate G4OG3OG bound to WT and W434A, based on the R-AXIS IV++ Imaging Plate detector. All data were collected at 0.5-1° WT:G4SG4OG4S complex oscillations throughout 180-720°. Data were processed with the The natural substrate G4OG3OG was modelled, based on the WT DENZO/SCALEPACK HKL 200065. Autoindexing determined that space HvExoI crystal structure in complex with G4SG4OG4S (PDB 6JG1), groups of crystals that were consistent with a tetragonal space group where, S-atoms in G4SG4OG4S were replaced by O-atoms, the gly- P43212 in all instances. The structures were refined using CCP4 cosidic linkage was re-arranged to be 3-O-linked between the 2nd and REFMAC566 and Phenix67. The previously solved structures of native 3rd Glc moieties from the non-reducing end, and the atoms were HvExoI: PDB 1EX1, 1IEQ, 1IEV, 1IEW, 1IEX, and 1J8V) without ligands, carefully moved such that the system maintained β-anomeric con- ions, glycerol, and water molecules were used as starting models2,3,22. formation. Further, SO42- anions, GOL, PEG, and NAG were removed, The iterative model building using REFMAC566, MIFIT68, and Coot69 and the system was neutralised with three Na+ cations and solvated in allowed tracing of all the residues in electron density maps. Following a cubic water box. From this structure, two more systems repre- convergence in a standard model refinement, approximately 2% senting the product of the hydrolytic reaction of G4OG3OG were improvements in the Rwork/Rfree factor ratios were achieved by refining modelled. The Glc product at the −1 subsite was removed in both Nature Communications | (2022)13:5577 14 Article https://doi.org/10.1038/s41467-022-33180-5 systems and for one system the WT structure was retained, while for laminaripentaose at the +1 and +2 subsites and pointing to bulk the other, Trp434 was mutated to Ala. These three systems were solvent). For a comparison, models for the apo and laminaritriose parametrised using the tleap program of AmberTools2179 using the substrate-bound forms of HvExoI were also generated. All models ff19SB protein force field80, Glycam0676 for saccharides, and the were built based on the crystal structures of WT (in complex with TIP3P force field77 for water molecules. Amber2079 was used to 4NP-G3SG; PDB 6JG2), W434A (in complex with 4NP-G3SG3OG; PDB minimise and equilibrate all systems using the protocol described 6L1J), and W434H (in complex with 4NP-G3SG; PDB 6JGN). The above, and applying restraints when necessary until the RMSD values starting structures for cMD simulations were prepared using USFC of protein backbones stabilised. Simulations were extended to Chimera95. Ion, water, GOL, and PEG molecules present in crystal 150 ns each. structures were removed and proteins were protonated at pH 7, and with the catalytic acid Glu491 in a neutral form as it corre- QM/MM MD metadynamics simulations of G3OG and G4OG3OG sponds to this step of the catalytic cycle1,2. The (1,3)-linked β-D- substrates in WT glucosides (laminaribiose, laminaritriose, and laminaripentaose) For these studies or conformational Free Energy Landscape (FEL) used in simulations were generated in the carbohydrate builder simulations of Glc moieties of G3OG and G4OG3OG positioned in GLYCAM, using the legacy version of the software (https://glycam. active sites of HvExoI, we used the CPMD code (parallelised plane org/cb/legacy). To ensure that the starting structures had Glc wave/pseudopotential implementation of Density Functional The- bound in the −1 subsite and the substrates non-reducing ends ory, Copyright IBM Corp 1990-2008)81 patched with the open-source occupied the +1 and +2 sites, oligosaccharides were spatially mat- community-developed PLUMED2 library82,83. For the conformational ched to substrates present in crystal structures of WT and FEL simulation of the G3OG disaccharide, the QM box (region) mutant HvExoI using the match commands in USFC Chimera. included 45 atoms corresponding to the disaccharide in the −1 and System solvation, coordinate, and topology file acquisitions +1 subsites; the MM region included the protein and water molecules (prmtop and inpcrd) were performed in tleap, included in the (85,478 MM atoms). For the conformational FEL simulation of the AMBERTools21 suite96. Systems were solvated in TIP3P water G4OG3OG trisaccharide, the whole sugar was not included in the QM molecules and neutralised with Na+ ions, using a cubic box, with box (region) but instead, the Glc molecule in the −1 subsite and the standard periodic boundary conditions set. The system was para- half-ring of Glc in the +1 subsite (with two monovalent pseudopo- meterised using the Amber ff99SB-ILDN force fields97 for proteins, tentials to saturate the QM region) were included in the QM box, GLYCAM_06j76 for the carbohydrates, and the TIP3P77 for water corresponding to 29 QM atoms in total. The other half-ring of the Glc molecules. Simulations were carried out in AMBER 1898 with CUDA molecule in the +1 subsite and the third Glc moiety placed in the acceleration. The cMD simulation protocol consisted of: (i) initial +2 subsite were described as MM atoms, which also included protein minimisation (105 steps); (ii) heating (100 ps, with protein back- and water molecules (11,5415 MM atoms)84 (Supplementary Fig. 15). bones restrained); (iii) NVT equilibration (100 ps, with restrained The G3OG and G4OG3OG systems were enclosed in isolated super- backbones); (iv) NPT equilibration (500 ps, unrestrained); and (v) cells of 13.5 Å x 19.3 Å x 14.8 Å and 13.3 Å3 sizes. Kohn-Sham orbitals NPT production at 300 K under a pressure of 1.0 atm. Every simu- were expanded on a plane wave basis set with a kinetic energy cut-off lation was carried out with a time step value of 2 ps using the of 70 Ry. Troullier-Martins ab initio pseudopotentials85 were used for SHAKE78 algorithm to constrain bonds involving hydrogens. A all elements. The PBE functional86 in the generalised gradient- cutoff value of 14.0 Å was used for non-bonding interactions and corrected approximation of DFT with a time step of 0.12 fs and a the Particle Mesh Ewald method for long-range electrostatics99. fictitious electron mass of 600 a.u. for the Car-Parrinello Lagrangian Production simulations were run for 1000 ns in all cases except were used in both systems. The Nosé-Hoover thermostat87 was cou- laminarihexaose, where we applied 600 ns simulation time. Due to pled to the systems to reach a constant temperature of 300 K. After the size of laminaripentaose, we also conducted GaMD simulation the systems were equilibrated, metadynamics simulations were to enhance the sampling of the conformational space adopted by initiated88. To obtain conformational FEL of Glc at the −1 subsite in proteins and ligands. both systems, we used the polar Cremer and Pople puckering coordinates89, as φ and θ collective variables. The well-tempered GaMD simulations of laminaritriose and laminarihexaose metadynamics algorithm90 was used with a bias factor of 12, a tem- hydrolysis products by WT, and the W434H and W434A mutants perature of 300 K, and an initial height of Gaussians of 0.75 kcal/mol. For GaMD simulations55, AMBER 1898 was used with the same forcefield The Gaussian width was 0.1 rad with a deposition time of 400 described in the previous paragraph. The GaMD protocol consisted of molecular dynamics steps (corresponding to 48 fs). The simulations an initial equilibration stage, where the potential boost was applied were stopped after a free energy error was below 0.6 kcal/mol (ana- and boost parameters were updated, and production runs were lysed using the block error analysis91), corresponding to 2822 Gaus- updated with fixed boost parameters. A dual boost on both dihedral sians (140 ps) and 2980 Gaussians (150 ps) for G3OG and G4OG3OG, and total potential energy was applied (igamd = 3). As a starting point, respectively (cf. Supplementary Fig. 11). the structures obtained after clustering in Chimera cMD trajectories were used. For W434H and W434A, clustering was performed on the Analyses of cMD and QM/MM MD metadynamics simulations first part of the trajectory since the oligosaccharides did not remain Data analyses were performed in VMD1.9.392 using home-made bound in the +1 and +2 subsites for extended time periods, and GaMDs Python3 scripts, the matplotlib library93, and the previously descri- production simulations were run for around 800 ns. For comparison, bed script94. The open-source, community-developed PLUMED GaMD simulations of bound laminaritriose and the apo form were also library83 version 2.7.282 was used for block analyses. The cpptraj tool examined using 400 ns production runs. GaMD inputs were generated from Ambertools21 was used to obtain representative structures of the following the recommendations from the developers (Y. Miao and R. most populated clusters in simulations. Walker, Gaussian accelerated molecular dynamics tutorial; https:// miaolab.ku.edu/GaMD/lib/GaMD_Amber-Tutorial.pdf). After the simu- cMD simulations of laminaritriose or laminarihexaose hydro- lations were completed, mdcrd files were converted into a dcd format lysis products by WT and the W434H and W434A mutants using the Amber Cpptraj tool and visualised in USCF Chimera. All For each enzyme form, a model was built with bound hydrolysed systems were checked for convergence100–102 by calculating RMSD, all- products from laminaritriose (i.e. Glc at the −1 subsite and G3OG at to-all RMSD, and principal component analyses values for each set of the +1 and +2 subsites) or laminarihexaose (Glc at the −1 subsite and simulations with a Python script. Nature Communications | (2022)13:5577 15 Article https://doi.org/10.1038/s41467-022-33180-5 GpathFinder simulations to explore Glc displacement routes in 6JGE, 6LIJ, 6LBB, and 6K6V. The PDB accession of W434F in complex WT, and the W434H and W434A HvExoI mutants with G2SG-OMe, 4NP-G3SG, (G4SG4O)3, and G6SG-OMe are 6JGG, This method uses the multi-objective genetic algorithm of the Gau- 6LC5, 6JGK, and 6LBV. The PDB accession of W434H in complex with diMM engine103, wherein genes are structural operators (e.g. atomic G2SG-OMe, 4NP-G3SG, G4SG-OMe, and G6SG-OMe are 6JGL, 6JGN, positions and displacements), and objectives (fitness functions in 6JGO, and 6JGP. The PDB accession of W434Y in complex with G2sG- evolutive algorithms) are structural and energetic descriptors (e.g. OMe, 4NP-G3SG, G4SG-OMe, and G6SG-OMe are 6JGQ, 6JGR, 6JGS, and force field functions, distances, angles, etc.). The GaudiMM engine 6JGT. Source data are provided with this paper. combines full-atom representations of a ligand and a protein and allows diverse levels of structural flexibility including internal rotation, Code availability global rotation and translation of ligands, and the exploration of Software applications described in Methods with associated refer- protein motions that are in a short-range (e.g. internal rotation, global ences were used without code modifications. For analyses of geome- rotation, and translation of ligands) or a long-range (e.g. rotameric trical parameters and MD simulations, we used in-house scripts, which exploration and normal mode analysis of the protein). Typically, in are available on GitHub (https://github.com/insilichem/utils/; https:// GpathFinder runs, fast energetic evaluations are needed, which are github.com/roviralab/utils) and ZENODO repositories106,107. provided by a docking scoring function, i.e. using Autodock Vina104 and steric clashes. GpathFinder’s philosophy relies on the conforma- References tional selection hypothesis, as collective motions (provided by the 1. Varghese, J. N., Hrmova, M. & Fincher, G. B. Three-dimensional sampling of the lowest energetic modes obtained from Normal Mode structure of a barley β-D-glucan exohydrolase, a family 3 glycosyl Analysis) are responsible for conformational changes of proteins hydrolase. Structure 7, 179–190 (1999). associated with their function, including ligand binding. Ligand dif- 2. Hrmova, M. et al. Catalytic mechanisms and reaction inter- fusion of Glc from the −1 subsite was calculated with GPathfinder for mediates along the hydrolytic pathway of plant β-D-glucan glu- six systems: WT, W434H, and W434A with or without laminaribiose (i- cohydrolase. Structure 9, 1005–1016 (2001). iii) or laminaripentaose (iv-vi) hydrolysis products bound in the +1 and 3. Hrmova, M. et al. Structural basis for a broad specificity in +2 subsites and pointing to bulk solvent for laminaripentaose, since higher plant β-D-glucan glucohydrolases. Plant Cell 14, the Glc molecule occupied the −1 subsite. In each case, the starting 1–22 (2002). structure calculations were chosen from the most representative 4. Lombard, V., Golaconda Ramulu, H., Drula, E., Coutinho, P. M. & clusters of cMD simulations described above, and GpathFinder cal- Henrissat, B. The carbohydrate-active enzymes database (CAZy) in culations were performed using standard genetic algorithm 2013. Nucleic Acids Res. 42D, 490–495 (2014). parameters37. Clashes and the scoring Autodock Vina function were 5. Hrmova, M. & Fincher, G. B. Dissecting the catalytic mechanism of used as objectives, together with genes representing rotamers and a plant β-D-glucan glucohydrolase through structural biology Normal Modes Parameters calculated with ProDy105. The pool of indi- using inhibitors and substrate analogues. Carbohydr. Res. 342, viduals consisted of a population of sixteen candidates that satisfied 1613–1623 (2007). the genetic algorithm for 800 generations. Calculations were analysed 6. Nazipova, A. et al. Forgotten actors: glycoside hydrolases during using GaudiView for visual inspection, together with the movie plug-in elongation growth of maize primary root. Front. Plant Sci. 12, in UCSF Chimera to inspect trajectories. Other analyses were per- 802424 (2022). formed through in-house UCSF Chimera built-in plug-ins, such as 7. Zolotnitsky, G. et al. Mapping glycoside hydrolase substrate sub- finding clashes/contacts and distances, to locally explore interactions sites by isothermal titration calorimetry. Proc. Natl Acad. Sci. USA of Glc with amino acid residues during egress pathways. 101, 11275–11280 (2004). 8. Litzinger, S. et al. Structural and kinetic analysis of Bacillus subtilis Unconstrained MD simulations to explore sliding of substrates N-acetylglucosaminidase reveals a unique asp-his dyad mechan- in WT, and the W434H and W434A HvExoI mutants ism. J. Biol. Chem. 285, 35675–35684 (2010). These simulations proceeded with laminaripentaose in the +1 and 9. Bacik, J.-P., Whitworth, G. E., Stubbs, K. A., Vocadlo, D. J. & Mark, +2 subsites and pointing to bulk solvent, and Glc removed from the B. L. Active site plasticity within the glycoside hydrolase NagZ −1 subsite. The systems were generated from the structures obtained underlies a dynamic mechanism of substrate distortion. Chem. from GaMD simulations and were run as described above for ternary Biol. 19, 1471–1482 (2012). complexes, with production runs of 200 ns each. 10. Davies, D. R., Abendroth, J., Staker, B. & Stewart, L. Seattle Additional methods are available in Supplementary Methods Structural Genomics Center for Infectious Disease (SSGCID). (Materials, Phylogenetic analysis, Organic Syntheses of 4NP-G3SG and Crystal structure of β-hexosaminidase 1 from Burkholderia cen- G4SG-OMe). ocepacia J2315. PDB accession 4G6C (2012). 11. Pozzo, T., Linares-Pasten, J., Karlsson, E. N. & Logan, D. T. Reporting summary Structural and functional analyses of β-glucosidase 3B from Further information on research design is available in the Nature Thermotoga neapolitana: a thermostable three-domain repre- Research Reporting Summary linked to this article. sentative of glycoside hydrolase 3. J. Mol. Biol. 397, 724–739 (2010). Data availability 12. Karkehabadi, S. et al. Biochemical characterization and crystal The data that support this study are available from the corresponding structures of a fungal family 3 β-glucosidase, Cel3A from Hypo- author upon request. The atomic coordinates and structure factors crea jecorina. J. Biol. Chem. 289, 31624–31637 (2014). have been deposited in the Protein Data Bank (www.pdb.org) with the 13. Suzuki, K. et al. Crystal structures of glycoside hydrolase family 3 following PDB accessions: WT HvExoI, and in complex with 4NP-G3SG, β-glucosidase 1 from Aspergillus aculeatus. Biochem. J. 452, and (G4SG4O)3 are 6JG2, and 6JG1, respectively. The PDB accession of 211–221 (2013). W286A in complex with G6SG-OMe is 6JG6, and W286F in complex 14. Kim, J. S., Yoon, B. Y., Ahn, J., Cha, J. & Ha, N. C. Crystal structure of with G2OG-OMe, 4NP-G3SG, and G6SG-OMe are 6JG7, 6JGA, and 6JGB. β-N-acetylglucosaminidase CbsA from Thermotoga neapolitana. The PDB accession of W286Y in complex with Glc, and G6SG-OMe are Biochem. Biophys. Res. Commun. 464, 869–874 (2015). 6JGC and 6JGD. The PDB accession of W434A in complex with Glc, 15. Nakatani, Y., Cutfield, S. M., Cowieson, N. P. & Cutfield, J. F. G2SG-OMe, 4NP-G3SG3OG, (G4SG4O)3, and G6SG-OMe are 6KUF, Structure and activity of exo-1,3/1,4-β-glucanase from marine Nature Communications | (2022)13:5577 16 Article https://doi.org/10.1038/s41467-022-33180-5 bacterium Pseudoalteromonas sp. BB1 showing a novel C-terminal 35. Nin-Hill, A. & Rovira, C. The catalytic reaction mechanism of the β- domain. FEBS J. 279, 464–478 (2012). galactocerebrosidase enzyme deficient in Krabbe disease. ACS 16. McAndrew, R. P. et al. From soil to structure, a novel dimeric β- Cat. 10, 12091–12097 (2020). glucosidase belonging to glycoside hydrolase family 3 isolated 36. Morais, M. A. B. et al. Two distinct catalytic pathways for GH43 from compost using metagenomic analysis. J. Biol. Chem. 288, xylanolytic enzymes unveiled by X-ray and QM/MM simulations. 14985–14992 (2013). Nat. Commun. 12, 367 (2021). 17. Yoshida, E. et al. Role of a PA14 domain in determining substrate 37. Sánchez-Aparicio, J. E. et al. GPathFinder: Identification of ligand- specificity of a glycoside hydrolase family 3 β-glucosidase from binding pathways by a multi-objective genetic algorithm. Int. J. Kluyveromyces marxianus. Biochem. J. 431, 39–49 (2010). Mol. Sci. 20, 3155 (2019). 18. Agirre, J. et al. Three-dimensional structures of two heavily 38. Akopian, T. N., Kisselev, A. F. & Goldberg, A. L. Processive degra- N-glycosylated Aspergillus sp. family GH3 β-D-glucosidases. Acta dation of proteins and other catalytic properties of the protea- Crystallogr. D72, 254–265 (2016). some from Thermoplasma acidophilum. J. Biol. Chem. 272, 19. Zmudka, M. W., Thoden, J. B. & Holden, H. M. The structure of DesR 1791–1798 (1997). from Streptomyces venezuelae, a β-glucosidase involved in mac- 39. Wilson, D. B. & Kostylev, M. Cellulase processivity. Methods Mol. rolide self-resistance. Protein Sci. 22, 883–892 (2013). Biol. 908, 93–99 (2012). 20. Mohsin, I., Poudel, N., Li, D. C. & Papageorgiou, A. C. Crystal 40. Vermaas, J. V. et al. The dissociation mechanism of processive structure of a GH3 β-glucosidase from the thermophilic cellulases. Proc. Natl Acad. Sci. USA 116, 23061–23067 (2019). fungus Chaetomium thermophilum. Int. J. Mol. Sci. 20, 41. Speirs, L. B. M., Rice, D. T. F., Petrovski, S. & Seviour, R. J. The 5962 (2019). phylogeny, biodiversity, and ecology of the chloroflexi in acti- 21. Hrmova, M. et al. Barley β-D-glucan exohydrolases with β-D-glu- vated sludge. Front. Microbiol. 10, 2015 (2019). cosidase activity. Purification and determination of primary 42. Hrmova, M. & Fincher, G. B. Purification and properties of three structure from a cDNA clone. J. Biol. Chem. 271, 5277–5286 (1996). (1,3)-β-d-glucanase isoenzymes from young leaves of barley 22. Streltsov, V. A. et al. Discovery of processive catalysis by an exo- (Hordeum vulgare). Biochem. J. 289, 453–461 (1993). hydrolase with a pocket-shaped active site. Nat. Commun. 10, 43. Tews, I. et al. Bacterial chitobiase structure provides insight into 2222 (2019). catalytic mechanism and the basis of Tay-Sachs disease. Nat. 23. Igarashi, K. et al. Traffic jams reduce hydrolytic efficiency of cel- Struct. Biol. 3, 638–648 (1996). lulase on cellulose surface. Science 333, 1279–1282 (2011). 44. Sulzenbacher, G., Schülein, M. & Davies, G. J. Structure of the 24. Horn, S. J. et al. Costs and benefits of processivity in enzymatic endoglucanase I from Fusarium oxysporum: native, cellobiose, degradation of recalcitrant polysaccharides. Proc. Natl Acad. Sci. and 3,4-epoxybutyl beta-d-cellobioside-inhibited forms, at 2.3 Å USA 103, 18089–18094 (2006). resolution. Biochemistry 36, 5902–5911 (1997). 25. Zakariassen, H. et al. Aromatic residues in the catalytic center of 45. Dolores Díaz, M., Fernández-Alonso, M. C., Cuevas, G., Cañada, F. chitinase A from Serratia marcescens affect processivity, enzyme J. & Jiménez-Barbero, J. On the role of aromatic-sugar interactions activity, and biomass converting efficiency. J. Biol. Chem. 284, in the molecular recognition of carbohydrates. A 3D view by using 10610–10617 (2009). NMR. Pure Appl. Chem. 80, 1827–1835 (2008). 26. Hrmova, M. & Fincher, G. B. Structure-function relationships of β- 46. Raju, R. K., Ramraj, A., Vincent, M. A., Hillier, I. H. & Burton, N. A. D-glucan endo- and exohydrolases from higher plants. Plant Mol. Carbohydrate-protein recognition probed by density functional Biol. 47, 73–91 (2001). theory and ab initio calculations including dispersive interactions. 27. Varghese, J. N. et al. Three-dimensional structures of two plant β- Phys. Chem. Chem. Phys. 10, 6500–6508 (2008). glucan endohydrolases with distinct substrate specificities. Proc. 47. Waters, M. L. Aromatic interactions in model systems. Curr. Opin. Natl Acad. Sci. USA 91, 2785–2789 (1994). Chem. Biol. 6, 736–741 (2002). 28. Parsiegla, G., Reverbel, C., Tardif, C., Driguez, H. & Haser, R. 48. Houser, J. et al. The CH-π Interaction in protein-carbohydrate Structures of mutants of cellulase Cel48F of Clostridium cellulo- binding: bioinformatics and in vitro quantification. Chemistry 26, lyticum in complex with long hemithiocellooligosaccharides give 10769–10780 (2020). rise to a new view of the substrate pathway during processive 49. Geronimo, I., Payne, C. M. & Sandgren, M. Hydrolysis and trans- action. J. Mol. Biol. 375, 499–510 (2008). glycosylation transition states of glycoside hydrolase family 3 β- 29. Nakamura, A., Okazaki, K. I., Furuta, T., Sakurai, M. & Iino, R. Pro- glucosidases differ in charge and puckering conformation. J. Phys. cessive chitinase is Brownian monorail operated by fast catalysis Chem. B. 122, 9452–9459 (2018). after peeling rail from crystalline chitin. Nat. Commun. 9, 50. Davies, G. J. et al. Snapshots along an enzymatic reaction coor- 3814 (2018). dinate: analysis of a retaining β-glycoside hydrolase. Biochemistry 30. Sørensen, A., Lübeck, M., Lübeck, P. S. & Ahring, B. K. Fungal beta- 37, 11707–11713 (1998). glucosidases: a bottleneck in industrial use of lignocellulosic 51. Biarnés, X., Nieto, J., Planas, A. & Rovira, C. Substrate distortion in materials. Biomolecules 3, 612–631 (2013). the Michaelis complex of Bacillus 1,3-1,4-β-glucanase: insight 31. Hamre, A. G., Jana, S., Reppert, N. K., Payne, C. M. & Sørlie, M. from first principles molecular dynamics simulations. J. Biol. Processivity, substrate positioning, and binding: the role of polar Chem. 281, 1432–1441 (2006). residues in a family 18 glycoside hydrolase. Biochemistry 54, 52. Biarnés, X. et al. The conformational free energy landscape of β-d- 7292–306 (2015). glucopyranose. implications for substrate preactivation in β- 32. Hemsworth, G. R. et al. Structural dissection of a complex Bac- glucoside hydrolases. J. Am. Chem. Soc. 129, teroides ovatus gene locus conferring xyloglucan metabolism in 10686–10693 (2007). the human gut. Open Biol. 6, 160142 (2016). 53. Davies, G. J., Planas, A. & Rovira, C. Conformational analyses of the 33. Harvey, A. J., Hrmova, M., De Gori, R., Varghese, J. N. & Fincher, G. reaction coordinate of glycosidases. Acc. Chem. Res. 45, B. Comparative modeling of the three-dimensional structures of 308–316 (2012). family 3 glycoside hydrolases. Proteins 41, 57–269 (2000). 54. Hrmova, M., Stratilová, B. & Stratilová, E. Broad specific xyloglu- 34. Ara, K. Z. G. et al. Characterization and diversity of the complete can:xyloglucosyl transferases are formidable players in the re- set of GH family 3 enzymes from Rhodothermus marinus DSM modelling of plant cell wall structures. Int. J. Mol. Sci. 23, 4253. Sci. Rep. 10, 1329 (2020). 1656 (2022). Nature Communications | (2022)13:5577 17 Article https://doi.org/10.1038/s41467-022-33180-5 55. Miao, Y., Feher, V. A. & McCammon, J. A. Gaussian accelerated 77. Jorgensen, W., Chandrasekhar, J., Madura, J., Impey, R. & Klein, M. molecular dynamics: unconstrained enhanced sampling and free Comparison of simple potential functions for simulating liquid energy calculation. J. Chem. Theory Comput. 11, water. J. Chem. Phys. 79, 926–935 (1983). 3584–3595 (2015). 78. Ryckaert, J. P., Ciccotti, G. & Berendsen, J. C. Numerical integra- 56. Xie, W. J., Asadi, M. & Warshel, A. Enhancing computational tion of the cartesian equations of motion of a system with con- enzyme design by a maximum entropy strategy. Proc. Natl Acad. straints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, Sci. USA 119, e2122355119 (2022). 327–341 (2017). 57. Saffarian, S., Collier, I. E., Marmer, B. L., Elson, E. L. & Goldberg, G. 79. Case, D. A. et al. Amber 2021. The University of California, San Interstitial collagenase is a Brownian ratchet driven by proteolysis Francisco (2021). of collagen. Science 306, 108–111 (2004). 80. Tian, C. et al. Ff19SB: Amino-acid-specific protein backbone 58. Kassem, S. et al. Artificial molecular motors. Chem. Soc. Rev. 46, parameters trained against quantum mechanics energy surfaces 2592–2621 (2017). in solution. J. Chem. Theory Comput. 16, 528–552 (2020). 59. Luang, S., Hrmova, M. & Ketudat Cairns, J. R. High-level expression 81. CPMD parallelized plane wave/pseudopotential implementation of barley β-d-glucan exohydrolase HvExoI from a codon- of Density Functional Theory, http://www.cpmd.org/; Copyright optimized cDNA in Pichia pastoris. Prot. Exp. Purif. 73, IBM Corp. 1990-2008. 90–98 (2010). 82. Tribello, G. A., Bonomi, M., Branduardi, D., Camilloni, C. & Bussi, G. 60. Luang, S., Ketudat Cairns, J. R., Streltsov, V. A. & Hrmova, M. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. Crystallisation of wild-type and variant forms of a recombinant β- 185, 604–613 (2014). d-glucan glucohydrolase from barley (Hordeum vulgare L.) by 83. PLUMED consortium. Promoting transparency and reproducibility macroseeding with wild-type native microcrystals and preliminary in enhanced molecular simulations. Nat. Methods 16, X-ray analysis. Int. J. Mol. Sci. 11, 2759–2769 (2010). 670–673 (2019). 61. Hrmova, M. & Fincher, G. B. Barley β-d-glucan exohydrolases. 84. Zhang, Y., Lee, T. S. & Yang, W. A pseudobond approach to Substrate specificity and kinetic properties. Carbohydr. Res. 305, combining quantum mechanical and molecular mechanical 209–221 (1998). methods. J. Chem. Phys. 110, 46 (1998). 62. Leatherbarrow, R. J. Data analysis and graphics program for the 85. Troullier, N. & Martins, J. L. Efficient pseudopotentials for plane- IBM PC. Erithacus Software Ltd. (Surrey, UK) (1998). wave calculations. Phys. Rev. B 43, 1993–2006 (1991). 63. Perella, F. EZ-FIT: a practical curve-fitting microcomputer pro- 86. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient gram for the analysis of the enzyme kinetic data on IBM-PC com- approximation made simple. Phys. Rev. Lett. 77, patible computers. Anal. Biochem. 174, 437–447 (1988). 3865–3868 (1996). 64. Hrmova, M., Varghese, J. N., Høj, P. B. & Fincher, G. B. Crystal- 87. Nosé, S. A molecular dynamics method for simulations in the lization and preliminary X-ray analysis of β-glucan exohydrolase canonical ensemble. Mol. Phys. 52, 255–268 (2006). isoenzyme ExoI from barley (Hordeum vulgare). Acta Crystallogr. 88. Laio, A. & Parrinello, M. Escaping free-energy minima. Proc. Natl D54, 687–689 (1998). Acad. Sci. USA 99, 12562–12566 (2002). 65. Otwinowski, Z. & Minor, W. A processing of X-ray diffraction data 89. Cremer, D. & Pople, J. A. A general definition of ring puckering collected in oscillation mode. Methods Enzymol. 276, coordinates. J. Am. Chem. Soc. 97, 1354–1358 (1975). 307–326 (1996). 90. Barducci, A., Bussi, G. & Parrinello, M. Well-tempered metady- 66. Marshudov, G. N. P. et al. REFMAC5 for the refinement of macro- namics: A smoothly converging and tunable free-energy method. molecular crystal structures. Acta Crystallogr. D67, Phys. Rev. Lett. 100, 020603 (2008). 355–367 (2011). 91. Trieste Tutorial: Metadynamics simulations with PLUMED; https:// 67. Adams, P. D. et al. PHENIX: a comprehensive Python-based system www.plumed.org/doc-master/user-doc/html/trieste-4.html. for macromolecular structure solution. Acta Crystallogr. D66, 92. Humphrey, W., Dalke, A. & Schulten, K. VMD: Visual molecular 213–221 (2010). dynamics. J. Mol. Graph. 14, 33–38 (1996). 68. McRee, D. E. XtalView/Xfit - a versatile program for manipulating 93. Hunter, J. D. Matplotlib: A 2D graphics environment. Comput. Sci. atomic coordinates and electron density. J. Struct. Biol. 125, Eng. 9, 90–95 (2007). 156–165 (1996). 94. Biarnés, X. F., Rovira, C. & Sanz Carrasco, F. Mecanismes de pre- 69. Emsley, P., Lohkamp, B., Scott, W. G. & Cowtan, K. Features and activació de substrat en 1,3-1,4-beta-glucanasa. Modelització development of coot. Acta Crystallogr. D66, 486–501 (2010). mitjançant dinàmica molecular de primers principis. TDX (Tesis 70. Painter, J. & Merritt, E. A. A molecular viewer for the analysis of TLS Doctorals en Xarxa). Universitat de Barcelona, Departament de rigid-body motion in macromolecules. Acta Crystallogr. D61, Química Física and Centre de Recerca en Química Teòrica (2008). 465–471 (2005). 95. Pettersen, E. F. et al. UCSF Chimera-a visualization system for 71. Brünger, A. T. Free R value: a novel statistical quantity for asses- exploratory research and analysis. J. Comput. Chem. 25, sing the accuracy of crystal structures. Nature 355, 1605–1612 (2004). 472–475 (1992). 96. Case, D. A. et al. Amber 2021, University of California, San Fran- 72. Brünger, A. T. X-PLOR Version 3.851: a system for crystallography cisco, USA (2021). and NMR. (New Haven, CT: Yale University Press) (1992). 97. Lindorff-Larsen, K. et al. Improved side-chain torsion potentials for 73. Laskowski, R. A., MacArthur, M. W., Moss, D. S. & Thornton, J. M. the Amber ff99SB protein force field. Proteins 78, PROCHECK - a program to check the stereochemical quality of 1950–1958 (2010). protein structures. J. Appl. Cryst. 26, 283–291 (1993). 98. Case, D. A. et al. AMBER 2018. (University of California, San Fran- 74. Case, D. A. et al. AMBER 11, University of California, San Fran- cisco, USA, 2018). cisco (2011). 99. Darden, T., York, D. & Pedersen, L. Particle mesh Ewald: An 75. Hornak, V. et al. Comparison of multiple Amber force fields and N⋅log(N) method for Ewald sums in large systems. J. Chem. Phys. development of improved protein backbone parameters. Proteins 98, 10089 (1993). 65, 712–725 (2006). 100. Smith, L. J., Daura, X. & van Gunsteren, W. F. Assessing equili- 76. Kirschner, K. N. et al. GLYCAM06: a generalizable biomolecular bration and convergence in biomolecular simulations. Proteins force field. Carbohydrates. J. Comput. Chem. 29, 622–655 (2008). 48, 487–96 (2002). Nature Communications | (2022)13:5577 18 Article https://doi.org/10.1038/s41467-022-33180-5 101. Grossfield, A. & Zuckerman, D. M. Quantifying uncertainty and Synchrotron Research Program, which is funded by the Commonwealth sampling quality in biomolecular simulations. Annu. Rep. Comput. of Australia under the Major National Research Facilities Program. Chem. 5, 23–48 (2009). 102. Grossfield, A. et al. Best practices for quantification of uncertainty Author contributions and sampling quality in molecular simulations [Article v1.0]. Living S.P. and S.F. synthesised thio-ligands; S.L. and J.R.K.C. implemented J. Comput. Mol. Sci. 1, 5067 (2018). mutagenesis; S.L. completed enzyme kinetics; S.L. and M.H. purified 103. Rodríguez-Guerra Pedregal, J., Sciortino, G., Guasp, J., and crystallised proteins; V.A.S. and M.H. collected X-ray data; S.L., Municoy, M. & Maréchal, J.-D. GaudiMM: a modular multi- V.A.S., and M.H. refined crystal structures; M.H. and J.G.S. performed objective platform for molecular modeling. J. Comput. Chem. phylogeny analyses; A.N.H., S.A.G., and C.R. performed cMD, and QM/ 38, 2118–2126 (2017). MM MD metadynamics calculations; X.F.L., J.-D.M., and L.M. imple- 104. Trott, O. & Olson, A. J. AutoDock Vina: improving the speed and mented cMD, GaMD, and GPathFinder simulations; X.F.L., A.N.H., and accuracy of docking with a new scoring function, efficient opti- M.H. generated animations and structural images. M.H. wrote the paper mization, and multithreading. J. Comput. Chem. 31, with contributions from S.L., S.F., X.F.L., A.N.H., J.-D. M., L.M., and C.R. 455–461 (2010). M.H. oversaw the project. 105. Bakan, A., Meireles, L. M. & Bahar, I. ProDy: protein dynamics inferred from theory and experiments. Bioinformatics 27, Competing interests 1575–1677 (2011). The authors declare no competing interests. 106. Luang, S. et al. The evolutionary advantage of an aromatic clamp in plant family 3 glycoside exo-hydrolases. Nat. Commun., Additional information ZENODO repository, https://doi.org/10.5281/zenodo. Supplementary information The online version contains 6989338 (2022). supplementary material available at 107. Luang, S. et al. The evolutionary advantage of an aromatic clamp https://doi.org/10.1038/s41467-022-33180-5. in plant family 3 glycoside exo-hydrolases. Nat. Commun., ZENODO repository, https://doi.org/10.5281/zenodo. Correspondence and requests for materials should be addressed to 6997562 (2022). Maria Hrmova. Acknowledgements Peer review information Nature Communications thanks Åsmund Funding was supported by the Australian Research Council Kjendseth Røhr and the other, anonymous, reviewer(s) for their con- (DP120100900) to M.H., and the Australian Synchrotron to M.H. and tribution to the peer review of this work. Peer reviewer reports are V.A.S. J.G.S. and M.H. thank the University of Adelaide and M.H. available. acknowledges the Huaiyin Normal University for support. J.R.K.C. and S.L. acknowledge funding from the Suranaree University of Technology, Reprints and permission information is available at the Synchrotron Light Research Institute, and the Thailand Research http://www.nature.com/reprints Fund (RSA6280073). A.N.H. and C.R. acknowledge funding from the Spanish Ministry of Science, Innovation, through MICINN/AEI/FEDER, UE, Publisher’s note Springer Nature remains neutral with regard to jur- PID2020-118893GB-100, the Spanish Structures of Excellence María de isdictional claims in published maps and institutional affiliations. Maeztu (MDM-2017-0767), and the CARBOCENTRE project, which is funded by the European Research Council under the European Union’s Open Access This article is licensed under a Creative Commons Horizon 2020 research and innovation program ERC-2020-SyG-95123. Attribution 4.0 International License, which permits use, sharing, C.R. and A.N.H. also thank the Barcelona Supercomputing Center and adaptation, distribution and reproduction in any medium or format, as Red Nacional de Supercomputación RES for computational resources at long as you give appropriate credit to the original author(s) and the Marenostrum IV and CTE-Power for technical support. S.A.-G thanks the source, provide a link to the Creative Commons license, and indicate if Spanish Ministry of Economy and Competitiveness for a Predoctoral changes were made. The images or other third party material in this Scholarship (BES-2012-051782). S.P. and S.F. recognise Glyco@Alps via article are included in the article’s Creative Commons license, unless ANR-15-IDEX-02, Carnot Institut PolyNat, Labex ARCANE, and CBH-EUR- indicated otherwise in a credit line to the material. If material is not GS through NR-17-EURE-0003, and thank NanoBio ICMG (UAR 2607) for included in the article’s Creative Commons license and your intended providing expertise for mass spectrometry and NMR analyses. J.-D.M., use is not permitted by statutory regulation or exceeds the permitted X.F.L., and L.M. acknowledge the Spanish Ministerio de Ciencia, Inno- use, you will need to obtain permission directly from the copyright vación y Universidades via the PRE2019-088412 fellowship, and the PID- holder. To view a copy of this license, visit http://creativecommons.org/ 2020-116861GB-I00 and PGC2018-098592-B-100 grants, and L.M. licenses/by/4.0/. acknowledges the Universitat Autònoma de Barcelona Talent Program. The use of the Australian synchrotron was supported by the Australian © Crown 2022 1 School of Agriculture, Food and Wine, and Waite Research Institute, University of Adelaide, Waite Research Precinct, Glen Osmond, SA, Australia. 2 Department de Química, Universitat Autònoma de Barcelona, Bellaterra, Spain. 3Department of Química Inorgànica i Orgànica and Institut de Química Teòrica i Computacional, Universitat de Barcelona, Barcelona, Spain. 4The Florey Institute, University of Melbourne, Victoria, Australia. 5School of Chemistry, Suranaree University of Technology, Nakhon Ratchasima, Thailand. 6Université Grenoble Alpes, Centre de Recherches sur les Macromolécules Végétales, Grenoble, France. 7Institut de Biotecnologia i de Biomedicina, Universitat Autònoma de Barcelona, Bellaterra, Spain. 8Institució Catalana de Recerca i Estudis Avançats, Barcelona, Spain. 9Jiangsu Collaborative Innovation Centre for Regional Modern Agriculture and Environmental Protection, School of Life Science, Huaiyin Normal University, Huai’an, China. 10These authors contributed equally: Sukanya Luang, Xavier Fernández-Luengo, Alba Nin-Hill. e-mail: maria.hrmova@adelaide.edu.au Nature Communications | (2022)13:5577 19