By means of QM(DFT)/MM metadynamics we have unraveled the hydrolytic reaction mechanism of Neisseria polysaccharea amylosucrase (NpAS), a member of GH13 family. Our results provide an atomistic picture of the active site reorganization along the catalytic double-displacement reaction, clarifying whether the glycosyl-enzyme reaction intermediate features an α-glucosyl unit in an undistorted ⁴C₁ conformation, as inferred from structural studies, or a distorted ¹S₃-like conformation, as expected from mechanistic analysis of glycoside hydrolases (GHs). We show that, even though the first step of the reaction (glycosylation) results in a ⁴C₁ conformation, the α-glucosyl unit undergoes an easy conformational change toward a distorted conformation as the active site preorganizes for the forthcoming reaction step (deglycosylation), in which an acceptor molecule, i.e., a water molecule for the hydrolytic reaction, performs a nucleophilic attack on the anomeric carbon. The two conformations (⁴C₁ ad E₃) can be viewed as two different states of the glycosyl-enzyme intermediate (GEI), but only the E₃ state is preactivated for catalysis. These results are consistent with the general conformational itinerary observed for α-glucosidases.

Introduction

Glycoside hydrolases (GHs) are enzymes responsible for the degradation or hydrolysis of glycosidic bonds in carbohydrates. GHs are of fundamental interest in glycobiology and glycomics (Gamblin et al., 2009) since they are responsible for the modification of polysaccharides and glycoconjugates involved in numerous biological processes such as cell-cell recognition and polysaccharide degradation for biofuel processing (Pauly and Keegstra, 2008). These enzymes also provide paradigms for enzymatic catalysis that extend beyond the bounds of carbohydrate chemistry (Koshland, 1953; Davies et al., 2012). GHs are systematically classified in 156 families (by March 2019) according to their sequence similarity (http://www.cazy.org) (Cantarel et al., 2009). Enzymes from the same family usually share the same catalytic mechanism, either with retention or inversion of the anomeric configuration (Figure S1), in which an oxocarbenium ion-like transition state is formed (Rye and Withers, 2000; Zechel and Withers, 2000). Enzymes from the same family also share the same conformational itinerary of the saccharide unit at the -1 subsite (hereafter named -1 sugar), in which this saccharide unit changes conformation following a straight line or a equatorial one in Stoddart's diagram (Figure S2), either the one representing the Northern of Southern projection of the puckering sphere (Iglesias-Fernandez et al., 2015). Because there are a few conformations that satisfy the stereoelectronic requirements for a stable oxocarbenium ion (⁴H₃,³H₄, B_2,5,^2,5B, ³E, E₃, ⁴E, and E₄) the number of itineraries is limited (Speciale et al., 2014; Ardèvol and Rovira, 2015). However, there are several GH families for which the itinerary has not been proven yet or it is controversial.

Members of GH13 family are retaining GHs that catalyze the cleavage of the α-glucosidic linkages, such as α-amylase, responsible for starch and amylose degradation. This family also gathers unique enzymes called amylosucrases that catalyze from sole sucrose the synthesis of α-D-glucopyranosyl homopolymers and oligomers accompanied with limited sucrose hydrolysis (Albenne et al., 2004). These enzymes are thus considered to be both glycoside hydrolases and transglucosylases. The substrate conformational itinerary from the Michaelis complex (MC, Figure S1) to the covalent glycosyl-enzyme intermediate (GEI), is expected to be identical.

Several X-Ray structures on GH13 enzymes have demonstrated that the conformation of the -1 sugar ring in the MC is ⁴C₁ (Fujimoto et al., 1998; Mirza et al., 2001; Skov et al., 2002). The lack of distortion in the conformation, unlike many GHs, is because the α-stereochemistry at the anomeric carbon assures that the leaving group (+1 sugar) is in a reactive orientation for the S_N2 reaction to take place, unlike e.g., β-glucosidases (Biarnés et al., 2006, 2011). In fact, α-glucosidases have been suggested to follow a ⁴C₁ → [⁴H₃]^‡ → ¹S₃ itinerary, which is the opposite as the ¹S₃ → [⁴H₃]^‡ → ⁴C₁ for β-glucosidases (Davies et al., 2012). These differences are being exploited in the field of inhibitor design (Beenakker et al., 2017).

Concerning the glycosyl-enzyme covalent intermediate (GEI) of the reaction, structural analyses of GH13 amylosucrase show that the -1 sugar bears a ⁴C₁ conformation (Jensen et al., 2004). In the case of α-amylase, the GEI trapped using sugar analogs exhibits a similar ⁴C₁ conformation (Zhang et al., 2009; Caner et al., 2016). QM/MM calculations of human pancreatic α-amylase captured the relevant role of the catalytic residues during the hydrolysis mechanism (Pinto et al., 2015), but the conformational itinerary of the -1 sugar was not investigated. On the other hand, the GEI of Bifidobacterium adolescentis sucrose phosphorylase, a GH13 enzyme, shows the -1 sugar distorted in a ¹S₃ conformation, suggesting a completely different conformational pathway (⁴C₁ → TS → ¹S₃) (Mirza et al., 2006). Being both enzymes from the same family and acting over the same substrate, this discrepancy is puzzling. The ⁴C₁ conformation of the GEI is also in contrast with the irreversible cyclosulphate inhibitors of α-glucosidases from other GH families, which show an unambiguous ¹S₃ conformation of the GEI analog (Artola et al., 2017). To solve this conundrum, we here uncover the conformational itinerary of the α-glucosyl unit of Neisseria polysaccharea amylosucrase (NpAS), a member of family GH13 for which a structure of the Michaelis complex with the natural substrate (sucrose) is available (Figure 1), during catalysis. Our simulations, performed by means of QM/MM metadynamics methods, show that the GEI can exhibit both a relaxed ⁴C₁ conformation or a distorted E₃-like conformation of the -1 sugar, depending on whether the catalytic water is properly placed and oriented for catalysis or it is on its way in.

FIGURE 1

Figure 1. NpAS in complex with sucrose. (A) Crystal structure from PDB 1JGI. (B) QM/MM optimized structure of the active site (the Glu328Gln mutation of the crystal structure was reverted).

Results

Michaelis Complex Structure

To model the Michaelis complex of NpAS, we reverted the Glu328Gln mutation of the X-ray structure of the NpAS-sucrose complex and performed classical MD simulations. The conformation of the α-glucosyl ring at the -1 sugar changed from ⁴C₁ to B_3,O during the first 6 ns of the MD simulation, due to an increase in the distance between the nucleophile and the anomeric carbon, but returned to ⁴C₁ during the remaining 8 ns, indicating that it is probably the most stable state. The ⁴C₁ conformation was maintained during the subsequent QM/MM MD equilibration.

To further ascertain which is the most favored conformation, we computed the conformational free energy landscape (FEL) of the α-glucosyl ring at the -1 subsite by QM/MM metadynamics, using the Cremer-Pople puckering coordinates θ/φ as collective variables. This is a well-tested approach that we have used with success to analyse the conformation of carbohydrates in isolation or in the active sites of GHs (Biarnés et al., 2007; Alonso-Gil et al., 2017). The computed FEL, shown in Figure 2, confirms that the most stable conformation is ⁴C₁, with a secondary minimum (4.7 kcal/mol above in free energy) corresponding to a B_3,O distorted conformation. The most relevant hydrogen bond interactions involving the -1 sugar for the most stable ⁴C₁ conformer are listed in Table 1 (see atom labeling in Figure 1B). This conformer shows a longer glycosidic bond distance (C1-O1 = 1.44 Å) compared to B_3,O, (1.42 Å) and the C1-O5 bond is shorter by 0.03 Å. Notably, the distance between the nucleophile oxygen (O_Asp393) and the anomeric carbon (C1) is shorter when the sugar is not distorted, with values of 3.10 Å (⁴C₁) and 3.30 Å (B_3,O), which would facilitate the nucleophilic attack by Asp286. Clearly, only the ⁴C₁ conformation of the substrate is preactivated for catalysis.

FIGURE 2

Figure 2. Conformational free energy landscape (FEL, in a Mercator representation) of the -1 sugar (α-glucosyl) of NpAS in complex with sucrose. Contour lines at 1 kcal/mol.

TABLE 1

Table 1. Calculated values of the most important interactions between the -1 α-glucosyl residue and the enzyme in the Michaelis complex (energies in kcal/mol, distances in Å).

Glycosylation Reaction

The simulation of the first step of the double displacement reaction (glycosylation) was initiated from a snapshot of the global minimum (⁴C₁ conformation) of the conformational FEL of the Michaelis complex. Three collective variables were used to drive the glycosylation reaction (Figure S3 and Methods section), representing the proton transfer (CV1), the nucleophilic attack (CV2) and the glycosidic bond cleavage (CV3).

The free energy landscape reconstructed from the metadynamics simulation (reaction FEL) is shown in Figure 3A (for further detail, Figure 3B shows a two-dimensional projection). The FEL exhibits two clear minima in opposite regions: reactants (Michaelis complex, MC) and products (glycosyl-enzyme intermediate, GEI), separated by a transition state (TS1). The rate-limiting step of NpAS is not known but we can assume it is the glycosylation step, as in most GHs acting on substrates with a sugar aglycone (Li et al., 2001). Under this assumption, the reaction free energy barrier (17.2 kcal/mol) is in very good agreement with the experimental value of 17.9 kcal/mol estimated from the room temperature rate constant (Potocki de Montalk et al., 2000).

FIGURE 3

Figure 3. (A) Free energy landscape of the glycosylation reaction, obtained from QM/MM metadynamics simulations with three CVs. (B) Two-dimensional projection of the FEL into the collective variables CV1 and CV3.

Figure 4 and Table 2 show the evolution of the main catalytic distances along the minimum energy pathway and the structures of the active site at the stationary states are represented in Figure 5 (top panels). The glycosylation reaction begins with the approach of Glu328 to the glycosidic oxygen. The distance between the Glu328 acid/base proton and the glycosidic oxygen (O1-H_Glu328) quickly decreases from 2.5 to 1.2 Å. Afterwards (1 in Figure 4), the acid/base residue transfers the proton to the glycosidic oxygen and the glycosidic bond (C1-O1) starts to increase. Protonation of the glycosidic oxygen takes place in the anti-configuration with respect to the C1-O5 bond, as common in α-glycosidases (Alonso-Gil et al., 2017). At the reaction transition state (TS1), the glycosidic bond is partially broken (1.98 Å) but the nucleophile has not yet started to attack the anomeric carbon (2.98 Å). Therefore, the reaction follows a dissociative mechanism and can be described as D_NA_N (Guthrie and Jencks, 1989; Schramm and Shi, 2001). Afterwards, the distance between the nucleophilic oxygen of Asp286 and the anomeric carbon decreases and the glycosyl-enzyme covalent bond forms. As it was found for endo-β-glucanase (Biarnés et al., 2011), the maximum oxocarbenium ion character does not occur at the TS but later on the reaction pathway (2 in Figure 4). At this point, both the nucleophile and the leaving group are well separated from the C1 atom and its charge is higher than the one at the TS (by 0.05 electrons, Figure S4). The calculations thus reveal that the glycosylation reaction features an early TS with respect to charge development. The glycosyl-enzyme intermediate is almost formed at 3 (C1-O_Asp286 = 2.14 Å). Finally, the hydrogen bond interaction between the free oxygen of Asp286 (O') and the 6-OH breaks and Asp286 rotates around the covalent bond with the -1 sugar, which collapses to an undistorted ⁴C₁ conformation (GEI).

FIGURE 4

Figure 4. Evolution of relevant active site distances during the glycosylation reaction catalyzed by NpAS.

TABLE 2

Table 2. Calculated values of the most relevant catalytic distances (in Å) and puckering coordinates (in degrees) and their standard deviations along the glycosylation minimum energy pathway.

FIGURE 5

Figure 5. Reaction mechanism of NpAS obtained in this work (MC, Michaelis complex; TS, transition state; GEI, glycosyl-enzyme intermediate; P, reaction products). Hydrogen atoms have been also omitted, except the ones attached to heteroatoms. For the sake of clarity, the leaving β-fructose molecule has been omitted from the deglycosylation reaction (bottom panels).

It is interesting to analyse in detail the conformation of the -1 sugar during the reaction, in relation with the hydrogen bond between the nucleophile residue (Asp286) and the 6-OH. The evolution of the sugar conformation during the glycosylation reaction, defined by the θ puckering coordinate, is shown in Figure 6. From the MC to the TS, the conformation evolves from ⁴C₁ to ⁴H₃/⁴C₁. Afterwards, the sugar adopts an E₃ envelope conformation (2), followed by a ⁴H₃ half-chair (3). At this point, a sudden conformational change toward ⁴C₁ takes place (θ changes from 58° to 20°), concomitant with rotation of Asp286 around the glycosyl-enzyme bond, and the disruption of the 6-OH hydrogen bond interaction. The covalent intermediate is already formed before the rotation and the two configurations of the GEI (with and without hydrogen bond with the 6-OH) are close in energy. However, only the one without the hydrogen bond, and thus with a ⁴C₁ conformation, corresponds to a minimum on the FEL of Figure 3. Thus, the conformation of the GEI after the first step of the double displacement reaction is ⁴C₁. In other words, the conformational itinerary is in principle cyclic, i.e., it starts and ends up in the same conformation (⁴C₁ → [⁴H₃/E₃]^‡ → ⁴C₁). This is in agreement with the conformations of the experimental structures of the MC and GEI of NpAS (both in ⁴C₁ conformation) (Mirza et al., 2001; Skov et al., 2002; Jensen et al., 2004). However, the computed itinerary differs from the expected one for α-glucosidases, for which a distorted GEI is expected (⁴C₁ → [⁴H₃]^‡ → ¹S₃) (Davies et al., 2012; Speciale et al., 2014). As we will see in the next section, both views can be reconciled if we consider the reorganization of the active site required to start the deglycosylation reaction.

FIGURE 6

Figure 6. Evolution of the conformation of the -1 sugar (α-glucosyl residue) during the glycosylation reaction.

Deglycosylation Reaction

As shown above, the glycosylation reaction leads to an active site configuration in which the -1 sugar is in a ⁴C₁ conformation. This configuration is not suitable for the second step of the double displacement reaction for two main reasons. First, the leaving group (C1-O_Asp286) is not in axially oriented, which makes the S_N2 reaction unfavorable. Second, there is no water molecule in the vicinity of the anomeric carbon. The nearest water molecule remains at ~6 Å from the anomeric carbon in the simulation (Figure 7, left panel), forming hydrogen bond interactions with the acid/base residue (Glu328) and an aspartate residue (Asp393). This water molecule is the best candidate to act as a nucleophile in the deglycosylation process, but it is still too far and not well oriented to attack the anomeric carbon. Most likely, there is an energy barrier to bring the water molecule to a reactive configuration.

FIGURE 7

Figure 7. Reorganization of the active site of NpAS in the GEI state as the catalytic water approaches the anomeric carbon. The leaving β-fructose (product of the glycosylation reaction) is not shown for clarity.

To obtain the reactive configuration for the deglycosylation reaction, we selected the above water molecule and approached it to the anomeric carbon with the metadynamics algorithm using the O_water-C1 distance as collective variable. In order to allow the C1-Asp286 bond to adopt an axial orientation (necessary for the water nucleophilic attack), a second collective variable that accounts for the hydrogen bond between Asp286 and the 6-OH was used. This simulation facilitated the identification of a configuration of the glycosyl-enzyme intermediate in which the water molecule is at ≈ 3.5 Å from C1 (Figure 7, right panel). The position of the water molecule is stabilized by hydrogen bond interactions with the 2-OH of the -1 sugar, Asp393, and the Glu328 catalytic residue. The new configuration, hereafter named GEI*, is separated from the initial GEI by a small free energy barrier (4.7 kcal/mol; Figure S5). Remarkably, the -1 sugar at GEI* is preactivated for catalysis, as it exhibits a distorted conformation (E₃) with a pseudo-axial orientation of the leaving group. Moreover, the water molecule is well oriented for nucleophilic attack, as it has the oxygen atom lone pairs pointing toward the anomeric carbon. Therefore, this configuration is the one we chose to start the modeling of the deglycosylation reaction. Interestingly, a distorted conformation of the GEI of lysozyme, a family 22 β-GH, was recently described by the Mulholland group based on QM(SCC-DFTB)/MM calculations (Limb et al., 2019). However, the GEI* states found here differs from the one of (Limb et al., 2019) in that the nucleophile carboxylate group exhibits the usual syn configuration with respect to the C1-O_Asp286 bond (τ_{C1−O−C−O} ≈ 0°). Differences in the substrate, active site and reaction stereochemistry might be the reason of the discrepancy.

The deglycosylation reaction was modeled using three collective variables that take into account the cleavage of the glycosyl-enzyme covalent bond (C1-O_Asp286), the attack of the water oxygen to the anomeric carbon (C1-O_wat), and the deprototonation of the water molecule by Glu328. The FEL reconstructed from the metadynamics simulation, shown in Figure 8, shows two minima corresponding to the activated glycosyl-enzyme intermediate (GEI*) and the hydrolysis products (P), which are 7.6 kcal/mol more stable. The transition state (TS2) is 13.3 kcal/mol higher in energy with respect to the GEI* state, consistent with deglycosylation not being rate limiting.

FIGURE 8

Figure 8. (A) Free energy landscape (FEL) of the deglycosylation reaction, obtained from QM/MM metadynamics simulations with three CVs. (B) Two-dimensional projection of the FEL into the collective variables CV1 and CV2.

The structure of the active site and the -1 sugar conformation along the deglycosylation reaction pathway is shown in Figure 5 (bottom panels) and Figure 9, whereas Table 3 lists the evolution of the most important distances. The deglycosylation reaction begins by cleavage of the glycosyl-Asp286 bond, followed by attack of the water molecule on the anomeric carbon, while the water forms a tight hydrogen bond with the acid/base residue (Glu328). The system overcomes the transition state (TS2) and the conformation of the -1 sugar changes to an E₃ envelope. Afterwards, the Glu328 abstracts a proton from the catalytic water, the covalent bond between the water molecule and the anomeric carbon forms and the -1 sugar adopts a ⁴C₁ conformation. Therefore, the conformational itinerary of the deglycosylation reaction of NpAS can be described as E₃ → [E₃]^‡ → ⁴C₁.

FIGURE 9

Figure 9. Evolution of the conformation of the -1 sugar (α-glucosyl residue) during the deglycosylation reaction.

TABLE 3

Table 3. Calculated values of the most relevant catalytic distances (in Å) and their standard deviations along the deglycosylation minimum energy pathway.

Summary and Conclusions

By means of QM(DFT)/MM metadynamics we have unraveled the reaction mechanism of α-amylosucrase during hydrolysis, in particular the conformational dynamics of the glycosyl-enzyme intermediate (GEI), for which structural studies have assigned a ⁴C₁ conformation. Our results show that the glycosylation reaction, assisted by the catalytic residues Asp286 (general base) and Glu328 (acid/base) features a dissociative transition state, in which the α-glucosyl residue at the -1 subsite adopts a ⁴H₃ conformation. The GEI intermediate adopts a ⁴C₁ conformation, thus being not preactivated for catalysis. However, the active site is very dynamic and it can easily evolve toward another configuration, in which the -1 sugar adopts a reactive E₃ conformation, with a pseudo-axial orientation of the leaving group, as the catalytic water enters the active site. The two conformations can be viewed as two different states of the GEI intermediate, dry or wet. We conclude that the catalytic itinerary of amylosucrase for the Michaelis complex → TS → intermediate enzymatic half reaction could be either described indistinctly as ⁴C₁ → [⁴H₃]^‡ → ⁴C₁ (cyclic itinerary) or ⁴C₁ → [⁴H₃]^‡ → E₃ (linear itinerary). This reconciles the results obtained by X-ray, with an undistorted GEI (⁴C₁), with the expected conformation of α-glucosidases (distorted GEI with a conformation near ¹S₃).

Computational Details

Model Building

The initial structure of the NpAS-sucrose complex was taken from the Michaelis complex structure of the Asp328Asn mutant (PDB 1JGI) (Mirza et al., 2001). The protonation states and hydrogen atom positions of all histidine residues were selected based on their hydrogen bond network (specifically, histidines 39, 173, 233, 332, 370, 392, 414, 463, 512, and 540 were considered neutral with the proton at N_ε, histidines 192, 306, 565 and 591 were neutral with a proton at N_δ and histidines 377, 382, and 601 were protonated). All Asp and Glu residues were taken as deprotonated (i.e., negatively charged) except Glu328 (the acid/base residue). The system was solvated with 20725 water molecules and 16 sodium ions were added to achieve neutrality of the protein structure, forming a rectangular box with dimensions of 89.1 Å × 106.4 Å × 81.1 Å.

Classical MD Simulations

Classical MD simulations of the NpAS-sucrose complex at room temperature were performed using the AMBER11 software (Case et al., 2010). The force-fields ff99SB (Hornak et al., 2006), GLYCAM06 (Kirschner et al., 2008) and TIP3P (Jorgensen et al., 1983) were used for the protein, sucrose substrate and solvent water, respectively. The system was equilibrated in several steps. First, all water molecules were relaxed with a gradient minimizer, holding the protein, and substrate fixed. Next, the whole system was allowed to relax. To gradually reach the desired temperature of 300 K, spatial constraints were initially added to the interactions between the protein and the substrate, while water molecules and sodium ions were allowed to move freely. The constraints were then removed and the whole system was allowed to reach the desired temperature. During all the process, the acid/base residue (Glu328) was very mobile, with the proton position evolving from cis to trans conformations and vice versa, thus its relative position with respect to the sugar was controlled with a smooth restraint (force constant of 5 kcal/mol/Ų). The simulation was pursued for 100 ps at constant pressure, allowing the cell volume to evolve, until the density stabilized (~1.06 g/cm³). The MD simulation was extended to 15 ns at constant volume without restraints until the system had reached equilibrium. A timestep of 1 fs was used, increasing it to 2 fs during the last 14 ns (using SHAKE). The −1 sugar conformation evolved from ⁴C₁ to B_3,O and returned to the original ⁴C₁ conformation during the last 7 ns of the simulation. Analysis of the trajectories was carried out by using standard tools of AMBER and VMD (Humphrey et al., 1996). The root-mean-square-deviation (RMSD) of the protein backbone atoms with respect to the crystal structure was stabilized around 1.4 Å in the equilibrated structure. One snapshot from the last 0.5 ns of simulation was taken as starting structure for the subsequent QM/MM MD simulations.

QM/MM MD Simulations

QM/MM MD simulations were performed with the CPMD software (Car and Parrinello, 1985), using the QM/MM interface developed by Laio et al. (2002). The QM region was considered as follows: (i) the whole sucrose molecule for the simulations of the conformational free energy landscape (FEL) of the -1 sugar (α-glucosyl); (ii) the sucrose, the acid/base residue (Glu328, capped at the C_β) and the nucleophile residue (Asp286, capped at the C_α) for the simulation of the glycosylation reaction; (iii) same as (ii) plus the catalytic water, for the simulation of the deglycosylation reaction. In all cases, the frontier atoms between QM and MM region were described using pseudopotential carbon link atoms. The fictitious electronic mass of the Car-Parrinello Lagrangian was taken as 600 au and the timestep was set at 0.12 fs in all CPMD simulations. All systems were enclosed in an isolated cubic box of 12.0 Å × 12.0 Å × 12.0 Å, using a fictitious electron mass of 700 au and a time step of 0.12 fs. The Kohn-Sham orbitals were expanded in a plane wave (PW) basis set with a kinetic energy cutoff of 70 Ry. Ab initio pseudopotentials generated within the Troullier-Martins scheme were employed (Troullier and Martins, 1991). The Perdew, Burke and Ernzerhoff generalized gradient-corrected approximation (PBE) (Perdew et al., 1996) was selected in view of its good performance in previous work on isolated sugars (Biarnés et al., 2007; Marianski et al., 2016), glycosidases (Jin et al., 2016) and glycosyltransferases (Bilyard et al., 2018).

QM/MM Metadynamics Simulations

QM/MM metadynamics (Laio and Parrinello, 2002; Barducci et al., 2011) simulations were performed to characterize the conformational FEL of the α-glucosyl residue of sucrose in the active site of NpAS and to simulate the different steps of the enzymatic reaction. The following collective variables were used: (i) Conformational FEL: the Cremer-Pople puckering coordinates (Cremer and Pople, 1975) phi and theta (φ, θ) of the α-glucosyl unit, following the methodology previously used in our group to rationalize and predict catalytic itineraries of GHs (Biarnés et al., 2007; Ardèvol and Rovira, 2015; Iglesias-Fernandez et al., 2015): (ii) glycosylation reaction: three collective variables representing the proton transfer (CV1), the nucleophilic attack (CV2) and the glycosidic bond cleavage (CV3) (Figure S3); (iii) deglycosylation reaction: three collective variables representing the covalent enzyme-substrate interaction (CV1), the water attack (CV2) and the proton transfer (CV3) were used (Figure S3). The metadynamics algorithm (Laio and Parrinello, 2002; Barducci et al., 2011), provided by the Plumed 2 plugin (Tribello et al., 2014), was used to explore the conformational free energy landscape of the systems. The height/width of the Gaussian terms was tested according to the oscillations of the CVs in free dynamics. In the case of the conformational FEL simulation, the height/width of the Gaussian terms was set at 0.75 kcal/mol/0.10 Å and a new Gaussian-like potential was added every 400 MD steps. The simulation of the conformational FEL was stopped once energy differences between the two wells (⁴C₁ and B_3,O) were maintained (1971 Gaussian terms were added), which was further tested by a time-independent free energy estimator (Tiwary and Parrinello, 2015). The error in the region of the two relevant minima and the pathway interconnecting them was ≤ 1.2 kcal/mol (Figure S6).

In the case of the glycosylation reaction, the height/width of the Gaussian terms was set at 1 kcal·mol⁻¹/0.20 Å (CV1) and 1 kcal/mol/0.10 Å (CV2 and CV3) and a new Gaussian-like potential was added every 300 MD steps. Walls for each CV at appropriate distances were used to reduce the FEL space to the chemical event. For the deglycosylation reaction, values of 1 kcal/mol/ 0.10 Å (height/width) and 250 MD steps (deposition time) were used. Metadynamics simulations were stopped after one crossing over the transition state (Figure S7), as recommended for chemical reactions (Ensing et al., 2005). Previous work on carbohydrate-active enzymes shows that the error associated to the metadynamics is <1 kcal/mol using this criteria (Raich et al., 2016). The total number of Gaussian terms added was 2284 (glycosylation reaction) and 4947 (deglycosylation reaction). The reaction coordinate was taken from the minimum free energy pathway, computed according to the intrinsic reaction coordinate method (Fukui, 1981). Structures at a given point along the reaction coordinate were taken from averages over a small region defined by CV1 ± 0.2, CV2 ± 0.2, CV3 ± 0.2 Å and were used for analysis.

Author Contributions

CR and IA designed calculations. SA-G performed and analyzed the simulations. JC helped with the computational work and analysis of the results. CR wrote the manuscript with the help of all authors.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by grants from the Spanish Minister of Science, Innovation and Universities (MICINN) (CTQ2017-85496-P to CR), the Agency for Management of University and Research Grants of Generalitat de Catalunya (AGAUR) (2017SGR-1189 to CR), the French National Research Agency (ANR) (CarbUniVax ANR-15-CE07-0019-01 to IA) and the Spanish Structures of Excellence María de Maeztu (MDM-2017-0767 to CR). The authors gratefully acknowledge the computer resources at and the technical support provided by the Barcelona Supercomputing Center (BSC-CNS, Barcelona, Spain), the TGCC-Curie supercomputer (Paris, France) and the mesocenter of Région Midi-Pyrénéées (CALMIP, Toulouse, France). SA-G and JC acknowledge MICINN for predoctoral fellowships (FPI-BES-2012-051782 and FPI-BES-2015-072055, respectively).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2019.00269/full#supplementary-material

References

Albenne, C., Skov, L. K., Mirza, O., Gajhede, M., Feller, G., D'amico, S., et al. (2004). Molecular basis of the amylose-like polymer formation catalyzed by Neisseria polysaccharea amylosucrase. J. Biol. Chem. 279, 726–734. doi: 10.1074/jbc.M309891200

PubMed Abstract | CrossRef Full Text | Google Scholar

Alonso-Gil, S., Males, A., Fernandes, P. Z., Williams, S. J., Davies, G. J., and Rovira, C. (2017). Computational design of experiment unveils the conformational reaction coordinate of GH125 alpha-mannosidases. J. Am. Chem. Soc. 139, 1085–1088. doi: 10.1021/jacs.6b11247

PubMed Abstract | CrossRef Full Text | Google Scholar

Ardèvol, A., and Rovira, C. (2015). Reaction mechanisms in carbohydrate-active enzymes: glycoside hydrolases and glycosyltransferases. Insights from ab initio quantum mechanics/molecular mechanics dynamic simulations. J. Am. Chem. Soc. 137, 7528–7547. doi: 10.1021/jacs.5b01156

PubMed Abstract | CrossRef Full Text | Google Scholar

Artola, M., Wu, L., Ferraz, M. J., Kuo, C. L., Raich, L., Breen, I. Z., et al. (2017). 1,6-cyclophellitol cyclosulfates: a new class of irreversible glycosidase inhibitor. ACS Cent. Sci. 3, 784–793. doi: 10.1021/acscentsci.7b00214

PubMed Abstract | CrossRef Full Text | Google Scholar

Barducci, A., Bonomi, M., and Parrinello, M. (2011). Metadynamics. Wiley Interdiscip. Rev. Comput. Mol. Sci. 1, 826–843. doi: 10.1002/wcms.31

CrossRef Full Text | Google Scholar

Beenakker, T. J. M., Wander, D. P. A., Offen, W. A., Artola, M., Raich, L., Ferraz, M. J., et al. (2017). Carba-cyclophellitols are neutral retaining-glucosidase inhibitors. J. Am. Chem. Soc. 139, 6534–6537. doi: 10.1021/jacs.7b01773

PubMed Abstract | CrossRef Full Text | Google Scholar

Biarnés, X., Ardèvol, A., Iglesias-Fernández, J., Planas, A., and Rovira, C. (2011). Catalytic itinerary in 1,3-1,4-beta-glucanase unraveled by QM/MM metadynamics. Charge is not yet fully developed at the oxocarbenium ion-like transition state. J. Am. Chem. Soc. 133, 20301–20309. doi: 10.1021/ja207113e

CrossRef Full Text | Google Scholar

Biarnés, X., Ardèvol, A., Planas, A., Rovira, C., Laio, A., and Parrinello, M. (2007). The conformational free energy landscape of beta-D-glucopyranose. Implications for substrate preactivation in beta-glucoside hydrolases. J. Am. Chem. Soc. 129, 10686–10693. doi: 10.1021/ja068411o

PubMed Abstract | CrossRef Full Text | Google Scholar

Biarnés, X., Nieto, J., Planas, A., and Rovira, C. (2006). Substrate distortion in the Michaelis complex of Bacillus 1,3-1,4-beta-glucanase. Insight from first principles molecular dynamics simulations. J. Biol. Chem. 281, 1432–1441. doi: 10.1074/jbc.M507643200

PubMed Abstract | CrossRef Full Text | Google Scholar

Bilyard, M. K., Bailey, H. J., Raich, L., Gafitescu, M. A., Machida, T., Iglesias-Fernandez, J., et al. (2018). Palladium-mediated enzyme activation suggests multiphase initiation of glycogenesis. Nature 563, 235–240. doi: 10.1038/s41586-018-0644-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Caner, S., Zhang, X., Jiang, J., Chen, H. M., Nguyen, N. T., Overkleeft, H., et al. (2016). Glucosyl epi-cyclophellitol allows mechanism-based inactivation and structural analysis of human pancreatic alpha-amylase. FEBS Lett. 590, 1143–1151. doi: 10.1002/1873-3468.12143

PubMed Abstract | CrossRef Full Text | Google Scholar

Cantarel, B. L., Coutinho, P. M., Rancurel, C., Bernard, T., Lombard, V., and Henrissat, B. (2009). The Carbohydrate-Active EnZymes database (CAZy): an expert resource for glycogenomics. Nucleic Acids Res. 37, D233–238. doi: 10.1093/nar/gkn663

PubMed Abstract | CrossRef Full Text | Google Scholar

Car, R., and Parrinello, M. (1985). Unified approach for molecular dynamics and density-functional theory. Phys. Rev. Lett. 55, 2471–2474. doi: 10.1103/PhysRevLett.55.2471

PubMed Abstract | CrossRef Full Text | Google Scholar

Case, D. A., Darden, T. A., Cheatham, T. E., Simmerling, C., Wang, J., Duke, R., et al. (2010). AMBER 11. San Francisco, CA: University of California.

Google Scholar

Cremer, D., and Pople, J. A. (1975). General definition of ring puckering coordinates. J. Am. Chem. Soc. 97, 1354–1358. doi: 10.1021/ja00839a011

CrossRef Full Text | Google Scholar

Davies, G. J., Planas, A., and Rovira, C. (2012). Conformational analyses of the reaction coordinate of glycosidases. Acc. Chem. Res. 45, 308–316. doi: 10.1021/ar2001765

PubMed Abstract | CrossRef Full Text | Google Scholar

Ensing, B., Laio, A., Parrinello, M., and Klein, M. L. (2005). A recipe for the computation of the free energy barrier and the lowest free energy path of concerted reactions. J. Phys. Chem. B 109, 6676–6687. doi: 10.1021/jp045571i

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujimoto, Z., Takase, K., Doui, N., Momma, M., Matsumoto, T., and Mizuno, H. (1998). Crystal structure of a catalytic-site mutant alpha-amylase from Bacillus subtilis complexed with maltopentaose. J. Mol. Biol. 277, 393–407. doi: 10.1006/jmbi.1997.1599

PubMed Abstract | CrossRef Full Text | Google Scholar

Fukui, K. (1981). The path of chemical reactions - the IRC approach. Acc. Chem. Res. 14, 363–368. doi: 10.1021/ar00072a001

CrossRef Full Text | Google Scholar

Gamblin, D. P., Scanlan, E. M., and Davis, B. G. (2009). Glycoprotein synthesis: an update. Chem. Rev. 109, 131–163. doi: 10.1021/cr078291i

PubMed Abstract | CrossRef Full Text | Google Scholar

Guthrie, R. D., and Jencks, W. P. (1989). IUPAC recommendations for the representation of reaction mechanisms. Acc. Chem. Res. 22, 343–349. doi: 10.1021/ar00166a001

CrossRef Full Text | Google Scholar

Hornak, V., Abel, R., Okur, A., Strockbine, B., Roitberg, A., and Simmerling, C. (2006). Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins 65, 712–725. doi: 10.1002/prot.21123

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: Visual molecular dynamics. J. Mol. Graph. 14, 33–38. doi: 10.1016/0263-7855(96)00018-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Iglesias-Fernandez, J., Raich, L., Ardèvol, A., and Rovira, C. (2015). The complete conformational free-energy landscape of β-xylose reveals a two-fold catalytic itinerary for β-xylanases. Chem. Sci. 6, 1167–1177. doi: 10.1039/C4SC02240H

PubMed Abstract | CrossRef Full Text | Google Scholar

Jensen, M. H., Mirza, O., Albenne, C., Remaud-Simeon, M., Monsan, P., Gajhede, M., et al. (2004). Crystal structure of the covalent intermediate of amylosucrase from Neisseria polysaccharea. Biochemistry 43, 3104–3110. doi: 10.1021/bi0357762

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Y., Petricevic, M., John, A., Raich, L., Jenkins, H., Portela De Souza, L., et al. (2016). A beta-mannanase with a lysozyme-like fold and a novel molecular catalytic mechanism. ACS Cent. Sci. 2, 896–903. doi: 10.1021/acscentsci.6b00232

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. doi: 10.1063/1.445869

CrossRef Full Text | Google Scholar

Kirschner, K. N., Yongye, A. B., Tschampel, S. M., Gonzalez-Outeirino, J., Daniels, C. R., Foley, B. L., et al. (2008). GLYCAM06: a generalizable biomolecular force field. Carbohydrates. J. Comput. Chem. 29, 622–655. doi: 10.1002/jcc.20820

PubMed Abstract | CrossRef Full Text | Google Scholar

Koshland, D. E. (1953). Stereochemistry and the mechanism of enzymatic reactions. Biol. Rev. Camb. Philos. Soc. 28, 416–436. doi: 10.1111/j.1469-185X.1953.tb01386.x

CrossRef Full Text | Google Scholar

Laio, A., and Parrinello, M. (2002). Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 99, 12562–12566. doi: 10.1073/pnas.202427399

PubMed Abstract | CrossRef Full Text | Google Scholar

Laio, A., Vandevondele, J., and Rothlisberger, U. (2002). A hamiltonian electrostatic coupling scheme for hybrid Car-Parrinello molecular dynamics simulations. J. Chem. Phys. 116, 6941–6947. doi: 10.1063/1.1462041

CrossRef Full Text | Google Scholar

Li, Y. K., Chir, J., and Chen, F. Y. (2001). Catalytic mechanism of a family 3 beta-glucosidase and mutagenesis study on residue Asp-247. Biochem. J. 355, 835–840. doi: 10.1042/bj3550835

PubMed Abstract | CrossRef Full Text | Google Scholar

Limb, M.a,.L, Suardiaz, R., Grant, I. M., and Mulholland, A. J. (2019). Quantum mechanics/molecular mechanics simulations show saccharide distortion is required for reaction in hen egg-white lysozyme. Chem. Eur. J. 25, 764–768. doi: 10.1002/chem.201805250

PubMed Abstract | CrossRef Full Text | Google Scholar

Marianski, M., Supady, A., Ingram, T., Schneider, M., and Baldauf, C. (2016). Assessing the accuracy of across-the-scale methods for predicting carbohydrate conformational energies for the examples of glucose and alpha-maltose. J. Chem. Theory Comput. 12, 6157–6168. doi: 10.1021/acs.jctc.6b00876

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirza, O., Skov, L. K., Remaud-Simeon, M., Potocki De Montalk, G., Albenne, C., Monsan, P., et al. (2001). Crystal structures of amylosucrase from Neisseria polysaccharea in complex with D-glucose and the active site mutant Glu328Gln in complex with the natural substrate sucrose. Biochemistry 40, 9032–9039. doi: 10.1021/bi010706l

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirza, O., Skov, L. K., Sprogoe, D., Van Den Broek, L. A., Beldman, G., Kastrup, J. S., et al. (2006). Structural rearrangements of sucrose phosphorylase from Bifidobacterium adolescentis during sucrose conversion. J. Biol. Chem. 281, 35576–35584. doi: 10.1074/jbc.M605611200

PubMed Abstract | CrossRef Full Text | Google Scholar

Pauly, M., and Keegstra, K. (2008). Cell-wall carbohydrates and their modification as a resource for biofuels. Plant J. 54, 559–568. doi: 10.1111/j.1365-313X.2008.03463.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Perdew, J. P., Burke, K., and Ernzerhof, M. (1996). Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868. doi: 10.1103/PhysRevLett.77.3865

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinto, G. P., Bras, N. F., Perez, M. A., Fernandes, P. A., Russo, N., Ramos, M. J., et al. (2015). Establishing the catalytic mechanism of human pancreatic alpha-amylase with QM/MM methods. J. Chem. Theory Comput. 11, 2508–2516. doi: 10.1021/acs.jctc.5b00222

PubMed Abstract | CrossRef Full Text | Google Scholar

Potocki de Montalk, G., Remaud-Simeon, M., Willemot, R. M., Sarcabal, P., Planchot, V., and Monsan, P. (2000). Amylosucrase from Neisseria polysaccharea: novel catalytic properties. FEBS Lett. 471, 219–223. doi: 10.1016/S0014-5793(00)01406-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Raich, L., Nin-Hill, A., Ardevol, A., and Rovira, C. (2016). Enzymatic cleavage of glycosidic bonds: strategies on how to set up and control a QM/MM metadynamics simulation. Methods Enzymol. 577, 159–183. doi: 10.1016/bs.mie.2016.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Rye, C. S., and Withers, S. G. (2000). Glycosidase mechanisms. Curr. Opin. Chem. Biol. 4, 573–580. doi: 10.1016/S1367-5931(00)00135-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Schramm, V. L., and Shi, W. (2001). Atomic motion in enzymatic reaction coordinates. Curr. Opin. Struct. Biol. 11, 657–665. doi: 10.1016/S0959-440X(01)00269-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Skov, L. K., Mirza, O., Sprogoe, D., Dar, I., Remaud-Simeon, M., Albenne, C., et al. (2002). Oligosaccharide and sucrose complexes of amylosucrase. Structural implications for the polymerase activity. J. Biol. Chem. 277, 47741–47747. doi: 10.1074/jbc.M207860200

PubMed Abstract | CrossRef Full Text | Google Scholar

Speciale, G., Thompson, A. J., Davies, G. J., and Williams, S. J. (2014). Dissecting conformational contributions to glycosidase catalysis and inhibition. Curr. Opin. Struct. Biol. 28C, 1–13. doi: 10.1016/j.sbi.2014.06.003

CrossRef Full Text | Google Scholar

Tiwary, P., and Parrinello, M. (2015). A time-independent free energy estimator for metadynamics. J. Phys. Chem. B 119, 736–742. doi: 10.1021/jp504920s

PubMed Abstract | CrossRef Full Text | Google Scholar

Tribello, G. A., Bonomi, M., Branduardi, D., Camilloni, C., and Bussi, G. (2014). PLUMED 2: new feathers for an old bird. Comput. Phys. Commun. 185, 604–613. doi: 10.1016/j.cpc.2013.09.018

CrossRef Full Text | Google Scholar

Troullier, N., and Martins, J. L. (1991). Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B Condens. Matter. 43, 1993–2006. doi: 10.1103/PhysRevB.43.1993

CrossRef Full Text | Google Scholar

Zechel, D. L., and Withers, S. G. (2000). Glycosidase mechanisms: anatomy of a finely tuned catalyst. Acc. Chem. Res. 33, 11–18. doi: 10.1021/ar970172

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, R., Li, C., Williams, L. K., Rempel, B. P., Brayer, G. D., and Withers, S. G. (2009). Directed “in situ” inhibitor elongation as a strategy to structurally characterize the covalent glycosyl-enzyme intermediate of human pancreatic alpha-amylase. Biochemistry 48, 10752–10764. doi: 10.1021/bi901400p

PubMed Abstract | CrossRef Full Text | Google Scholar