Calculation of substrate binding affinities for a bacterial GH78 rhamnosidase through molecular dynamics simulations.

Ram2 from Pediococcus acidilactici is a rhamnosidase from the glycoside hydrolase family 78. It shows remarkable selectivity for rutinose rather than para-nitrophenyl-alpha-l-rhamnopyranoside (p-NPR). Molecular dynamics simulations were performed using a homology model of this enzyme, in complex with both substrates. Free energy calculations lead to predicted binding affinities of -34.4 and -30.6 kJ mol-1 respectively, agreeing well with an experimentally estimated relative free energy of 5.4 kJ mol-1. Further, the most relevant binding poses could be determined. While p-NPR preferably orients its rhamnose moiety toward the active site, rutinose interacts most strongly with its glucose moiety. A detailed hydrogen bond analysis confirms previously implicated residues in the active site (Asp217, Asp222, Trp226, Asp229 and Glu488) and quantifies the importance of individual residues for the binding. The most important amino acids are Asp229 and Phe339 which are involved in many interactions during the simulations. While Phe339 was observed in more simulations, Asp229 was involved in more persistent interactions (forming an average of at least 2 hydrogen bonds during the simulation). These analyses directly suggest mutations that could be used in a further experimental characterization of the enzyme. This study shows once more the strength of computer simulations to rationalize and guide experiments at an atomic level.

Introduction

Glycosidases (glycoside hydrolases, EC 3.2.1) catalyze the hydrolysis of the terminal sugar residues from the non-reducing ends of a wide range of glycosylated compounds. Glycosidases are ubiquitous in all domains and exert diverse metabolic functions [1]. The concept of distinguishing structurally distinct glycoside hydrolase (GH) families introduced by Henrissat and Davies [2] has undoubtedly become the currently most significant classification system for glycosidases. CAZy (The Carbohydrate-Active enZYmes database) [3] provides reliable information about the putative function and reaction mechanisms of all currently available GH sequences [4]. The predicted 3-dimensional fold of a GH family (i.e., of the catalytic domain) is a reliable indicator of its glycone specificities (e.g. glucosidase, galactosidase, rhamnosidase, etc.). Unfortunately, the molecular principles determining the specificities and affinities for the aglycone, which is the reducing terminus in the case of saccharides, are at present far less understood, especially as aglycone specificities of glycosidases are generally broad and therefore a highly complex issue [1]. Accordingly, it is to date not possible to reliably predict the aglycone-specificity of a given glycosidase class based on its sequence or structure. This gap in knowledge is certainly due to the reported broad and complex functional diversity of GHs, and especially, the lack of sufficient experimental data. A key complication in this regard is that substrates that would be required for a broad biochemical characterization are often commercially unavailable (e.g. phospho-glycosidases) or too costly to allow detailed kinetic analyses.

The interest in glycosidases reaches far beyond understanding the metabolic functions of these enzymes: For example, the release of metabolites from glycosylated precursors through microorganisms involved in food fermentations or digestion is currently a subject of much debate regarding its health implications. As such, it is scarcely understood which bacterial glycosidase classes can be made responsible for increasing the bioavailability of health beneficial compounds (e.g. phenols) and/or dietary toxins masked through glycosylation, such as mycotoxins [5].

Experiments with alpha-l-rhamnosidases of Pediococcus acidilactici suggest two subfamilies of GH78 differing in their substrate affinity. Ram2 is more specific for the disaccharide rutinose while Ram presents higher affinity for aryl rhamnosides such as para-nitrophenyl-alpha-l-rhamnopyranoside (p-NPR). Experiments show that glucose is a strong inhibitor for Ram2 and rhamnose for Ram [6]. The specificity and affinity for substrates is essential to understand the differences in aglycone selectivity. The hydrolytic potential of rhamnosidase Ram2 is also not well understood since the mechanism leading to the hydrolysis of specific substrates is unclear. Generally, GH 78 rhamnosidases catalyze the hydrolysis of the glycosidic bond by the inverting mechanism [2,3], in which the aglycone is replaced by a water molecule in a single displacement step, resulting in inversion of the chirality of the anomeric carbon atom of rhamnose (Fig. 1A). The two carboxyl residues that serve as catalytic acid/base are on average 10 Å apart to allow the binding of both water and aglycone molecule between them [7]. The identity of the acid/base catalysts is not yet resolved for GH 78 rhamnosidases. Based on crystallography studies of RhaB (Bacillus sp. GL1), Cui et al. [8] proposed that Glu572 acting in combination with Asp567 or Glu841 would be the key residues. In respective order, these residues correspond to Asp222, Asp217 and Glu488 in Ram2.

Fig. 1

(A) General reaction mechanism of glycoside hydrolase family 78 α-l-rhamnosidases. R represents the aglycone moiety, which is β-d-glucopyranoside in the case of rutinose and p-nitrophenol in the case of p-NPR. (B) Chemical structure of rutinose (left) and p-NPR (right). Rutinose and p-NPR share the rhamnose moiety but differ in the second moiety which is glucose for rutinose and a nitrophenol ring for p-NPR.

Obviously, there is a strong need to obtain more detailed insights into the underlying principles of glycosidase selectivity. The crystallization of bacterial rhamnosidase revealed insight into the nature of the catalytic site by identifying several amino acids conserved within GH78 which may be involved in substrate recognition or hydrolysis. These residues are Asp217, Asp222, Trp226, Asp229 and Glu488 [8]. Computer simulations have become quite attractive to complement and support regular experiments. In silico experiments can help to rationalize and predict experimental results and have evolved into a mature research tool, complementary to experimental methods [9]. To achieve that, simulations have to be long enough such that the structural diversity is captured sufficiently. Due to the increased availability of computational tools and increasing computing capacities, molecular modeling and simulation gained considerable attractiveness over the last years. Simulation might then lead to predictions which help to focus the design of experiments [10].

The objective of the present study is to assess the substrate binding affinity of a GH78 alpha-l-rhamnosidase (Ram2) from P. acidilactici through molecular dynamics (MD) simulations. Ram2 was previously characterized experimentally, showing that this enzyme possesses remarkably distinct affinities for the disaccharide rutinose and the model substrate p-NPR [6]. MD Simulations of Ram2 in complex with rutinose and p-NPR (Fig. 1B) confirm these results and also reproduce previous findings on the function of GH78 rhamnosidases as reported by Cui et al. [8]. Accordingly, the results presented here clearly suggest that MD simulations are a useful complementary approach to estimate substrate binding affinities of glycosidases, and also to predict which amino acid residues might be involved in substrate recognition.

Methods

Since there is no crystal structure of the Ram2 protein available, a homology model of the catalytic domain of rhamnosidase 2 (GenBank accession ZP_07366943, starting at residue 192) from P. acidilactici was made. The crystal structure of alpha-rhamnosidase from Bacteroides thetaiotaomicron (PDB: 3CIH) with a sequence identity of 24% and a similarity of 42% was used as template. The homology model was made with SWISS Model [11] in the automated mode. The accuracy of the homology model was verified with the tools Verify3D [12,13] and PROCHECK [14-16]. Verify3D visualizes the quality of a structure by comparing the 3D structure of a protein with its amino acid sequence and using a collection of good structures as reference. The average 3D–1D score of all 334 residues of the Ram2 structure is greater than zero and fluctuates between 0.2 and 0.8. PROCHECK validates the stereo chemical quality of a structure, in this case by comparing it with 118 structures with a resolution of at least 2.0 Å and an R-factor greater than 20%. 99% of the residues are located in the allowed regions. Of these, 85% are in the most favorable regions, 12% in the additional allowed regions and 2% in the generously allowed regions. Only 1% of the residues have disallowed φ,ψ value combinations. PROCHECK also shows that 99.5% of the bond lengths and 95% of the bond angles are within the limits, which are considered to be less than 2.0 standard deviations from the mean. Additionally, the root mean square deviation (rmsd) between the modeled structure and the template structure was calculated using a structural alignment in the program PyMOL [17] to gain a minimized rmsd of 1.6 Å.

The molecular operating environment (MOE 2011.10) program [18] was used to dock the substrates into the active site. First, the site finder tool of MOE was used to find possible binding sites. The most prominent suggestion was in agreement with the catalytic site reported by Cui et al. [8]. Rutinose was docked into this site. The protein was rigid during the docking, while the ligand was flexible. During the docking a conformational analysis was performed which led to conformations of the ligand with favorable torsion angles of the rotatable bonds while keeping the bond lengths and bond angles the same. The placement methodology used was Triangle Matcher. This method aligns ligand atom triplets with alpha sphere triplets to generate different poses of the ligands. The initial collection of poses was refined using the Amber99 forcefield [19]. The scoring was done with the placement method and rescored by the London dG scoring function to evaluate favorable contacts such as hydrogen bonds based on geometric accuracy as well as hydrophobic and ionic interactions. This scoring function calculates an estimate of the free energy of binding of each pose by taking the rotational and translational entropy into account. Also, the flexibility loss of the ligand is estimated based on the topology. Furthermore, the desolvation energy is considered to quantify the preference of the ligand for the pocket (MOE 2011.10) [18]. After docking, the list of various poses was reduced by only considering the ones with a scoring value lower than −30 kJ mol−1. The remaining ones were visually inspected and the one which fittest best in the pocket chosen for the initial structure.

Since the homology model may not reflect the optimal functional conformation of the protein the resulting poses are only considered to be initial approximations of the complexes. Therefore MD simulations were started from the docking pose (pose 1) as this allows for both protein and ligand flexibility and an ensemble of conformations reflecting the behavior of a system can be generated. Additionally rutinose was positioned in the pocket such that the rhamnose moiety of rutinose points to the conserved catalytic residues proposed by Cui et al. [8] (pose 2). These conserved catalytic residues were determined by performing a sequence alignment of several rhamnosidases of the GH78 family from different bacteria (see Supplementary data 1). p-NPR was positioned in the pocket of the protein according to the two poses of rutinose by overlaying the rhamnose moieties of the two substrates.

Atomic interactions for the protein–substrate complex were calculated using the GROMOS 54a7 force field [20]. Molecular topology building blocks for the substrates are available in Supplementary data 2. All simulations were performed using the GROMOS11 suite of simulation programs [21].

The complex of the protein with the ligand was placed inside a rectangular box of water with a size of approximately 8 nm × 8 nm × 8 nm which contains approximately 16 000 SPC water molecules [22]. Fourteen Na+ counter ions were added to obtain a neutral system. The ligands were also simulated free in water solution in order to obtain the interaction energies of the ligands free in water, necessary to calculate the free energy of binding. Here, the box size was about 3 nm × 3 nm × 3 nm with approximately 1000 water molecules. Before simulation, all computational boxes were energy minimized to relax unfavorable contacts between the protein–ligand complex and the solvent. For this a steepest descent energy minimization with an energy threshold of 0.1 kJ mol−1 was performed. Subsequently thermalization and equilibration steps were performed, gradually increasing the temperature from 60 K to the simulation temperature of 300 K and decreasing the force constant of an initial position restraint on all solute atoms from 2.5 × 104 kJ mol−1 nm−2 to zero.

The system was set to a constant temperature of 300 K with a relaxation time of 0.1 ps and to a constant pressure of 1 atm with a relaxation time of 0.5 ps [23]. The isothermal compressibility was set to 4.575 × 10−4 (kJ mol−1. nm−3)−1. The SHAKE algorithm [24] was used with a relative geometric accuracy of 10−4 constraining all bonds to their optimal bond lengths. For the calculation of the non-bonded interactions a triple-range cutoff was used. For the long ranged cutoff 1.4 nm and for the short ranged cutoff 0.8 nm were chosen to calculate interactions from a grid-based pairlist. A reaction-field contribution [25] with a dielectric permittivity of 61 [26] was added to the electrostatic interactions and forces to account for the contributions of a homogenous medium outside the long-range cutoff distance.

The leap frog algorithm [27] with a time step of 2 fs was used. For both poses of both substrates, four simulations of 5 ns each were performed, starting from different randomly assigned initial velocities. Simulations of this length are not expected to sample the complete relevant phase space, but in order to capture the structure of protein–ligand complexes multiple distinct simulations often prove more useful than performing fewer long simulations.

After the simulation several analyses were performed. First the stability of the protein during the simulation was checked with the GROMOS++ [28] program dssp which shows the secondary structure over time according to the Kabsch and Sander rules [29]. Also several frames of the simulations were visualized to observe the behavior of the protein with the ligand during the simulation. To visualize interactions, hydrogen bonds between the protein and the ligand were determined using a geometric criterion. A hydrogen bond is defined as such when the distance between the acceptor and the hydrogen atom is less than 0.25 nm and the donor–hydrogen–acceptor angle is at least 135° [30]. Hydrophobic interactions were monitored by calculating the distance between the center of mass of aromatic moieties and the center of mass of both sugar moieties of the ligand or the nitrophenol ring during the whole simulation.

The free energy of binding was calculated using the Linear Interaction Energy (LIE) approach [31]. In LIE, the binding free energy ΔG is estimated using

In words, the free energy of binding is calculated from ensemble the averages (〈 〉) of the Van der Waals and electrostatic 〈V〉 interaction energies between the ligand and its surroundings (l-s) from the simulations with the ligand bound to the protein (bound) and free in water (free), scaled by the parameters α and β, respectively. The LIE Method only uses the endpoints of binding to calculate “solvation” free energies according to the linear response theory. For the parameters α and β, the values from the original parameterization were used: α was set to 0.16 and β to 0.5 [31].

Additionally, the relative weights of individual simulations and the binding poses were calculated based on their separate free energies. For the assignment of the relative weights first the solvation free energies of the bound and the free ligand ΔG(free) for a given pose or simulation i needs to be calculat ed using Eqs. (2) and (3) [32].

can subsequently be used in Eq. (4) to calculate the relative weight of each simulation where k is Boltzmann's constant and T is the absolute temperature.

Eqs. (5) and (6) are subsequently used to calculate the averages of the Van der Waals and electrostatic energies, weighted over the different poses. These averages were then used to calculate the overall free energies of binding according to the LIE method by using Eq. (1).

The calculated binding free energies are compared to experimental estimates obtained from the substrate concentrations at which the enzyme activity is 50% of its maximum value [6].

Results

The stability of the protein in complex with the substrates is depicted by the atom-positional root-mean-square deviations from the initial structures in Supplementary data 3 and the secondary structure over time for the first simulations of poses 1 and 2 of rutinose and p-NPR in Fig. 2. These are representative for the other simulations. The helices in all simulations are very stable. There are four short β strands which are formed mostly in all simulations but do show some reversible instability. In the first simulation of the two poses of both ligands the helices are very stable. The β strand shown in panel B between amino acids Lys321 and Val331 is not observed in the initial structure but forms during the simulation. This β strand is also formed in some of the other simulations. In the simulation of rutinose in pose 2 it is observed very stably after 1.6 ns till the end of the simulation. The two β strands between Asp222 and Ile227 and between Ile264 and Phe269 are only stable during the first simulation with rutinose in pose 1.

Fig. 2

Stability of the secondary structure over time of Ram2 simulated with rutinose and p-NPR respectively. The first simulation of all poses of both ligands was taken as representative. (A) Pose 1 of rutinose, (B) pose 2 of rutinose, (C) pose 1 of p-NPR, (D) pose 2 of p-NPR. Helices including 310, α and π are shown in red and β strands in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

A sequence analysis of several rhamnosidases of different bacteria was performed to identify catalytic residues conserved within the GH78 family (Supplementary data 1). Some of these conserved catalytic residues are involved in interactions such as hydrogen bonds or hydrophobic stacking. But also other amino acids are involved in hydrogen bond interactions. See Fig. 3 for a general overview of the protein structure and a close up of the active site.

Fig. 3

Structure of Ram2 with important residues highlighted. Amino acids forming the catalytic site are colored in red. Residues which are involved in hydrogen bond formation or hydrophobic stacking (Trp226) are shown in yellow. (A) Top view and (B) catalytic site. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The calculation of the free energy according to the LIE Method (Table 1) shows that the free energy of binding is different for the different simulations. Using Eq. (4), the binding free energies are transformed into weights for the individual simulations reported in Table 2. It can be seen that pose 1 of rutinose has a lower free energy and higher weights than pose 2. The free energy values of p-NPR in pose 1 are more constant within the different simulations while the largest fluctuations are observed for p-NPR in pose 2, with values ranging from −0.6 kJ mol−1 to −31 kJ mol−1. This correlates with the fact that the ligand positions itself very differently in the protein during these simulations. In particular during simulations with the favorable values of ΔG the ligands are located similarly in the protein while in the other two simulations the ligand has partly moved out of the active site. Fig. 4 shows a superposition of the four poses with the highest weights from which it can be seen that all ligands are positioned in the same area.

Table 1

Average non bonded interaction energies between the ligand and its surroundings, separated into van der Waals , and electrostatic contributions, together with the binding free energies , calculated using the LIE approach (Eq. (1)). Sim indicates the simulation and pose one of the two poses. All values in kJ mol−1.

	Water	Pose 1				Pose 2
		Sim1	Sim2	Sim3	Sim4	Sim1	Sim2	Sim3	Sim4
Rutinose
(Vl−svdw)	−33.3	−73.0	−68.1	−62.8	−73.8	−80.6	−61.3	−92.0	−77.4
(Vl−sel)	−391.7	−407.1	−430.0	−446.8	−449.6	−407.5	−429.5	−423.1	−406.5
(ΔGbindi)		−14.1	−24.7	−32.3	−35.4	−15.5	−23.4	−25.1	−14.5
p-NPR
(Vl−svdw)	−72.9	−131.2	−120.7	−121.7	−121.8	−133.0	−122.0	−109.6	−118.6
(Vl−sel)	−197.0	−203.6	−220.1	−226.4	−215.5	−178.9	−205.3	−244.8	−245.2
(ΔGbindi)		−12.6	−19.2	−22.5	−17.1	−0.6	−12.0	−29.8	−31.4

Table 2

Relative weights of the simulations, calculated according to Eq. (4). Pose indicates the different positions of the ligands from which the simulations were started. Sim denotes the different simulation runs.

	Sim1	Sim2	Sim3	Sim4
Rutinose
Pose 1	0.0001	0.0104	0.2168	0.7539
Pose 2	0.0003	0.0061	0.0122	0.0002
p-NPR
Pose 1	0.0003	0.0048	0.0181	0.0020
Pose 2	0.0000	0.0003	0.3298	0.6447

Fig. 4

The positions of rutinose and p-NPR showing the highest relative weights in the protein. The third simulation with pose 1 of rutinose is colored in green. The fourth simulation with pose 1 of rutinose is colored in yellow. The third simulation with pose 2 of p-NPR is colored in violet. The fourth simulation with pose 2 of p-NPR is colored in pink. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The highest relative weights of the energies are observed for pose 1 of rutinose and pose 2 of p-NPR. The fact that for the rutinose simulations pose 1 gives the highest relative weights agrees with the observation that the glucose moiety points to the conserved catalytic residues and makes stronger interactions. Also, the lower free energies of pose 2 of p-NPR over pose 1 are not unexpected because in pose 2 rhamnose is close to these residues rather than the nitrophenol ring, which is not much involved in interactions. These observations agree with the fact that glucose is a stronger inhibitor of the enzyme than rhamnose [6].

For the binding free energy ΔG in Table 3 all simulations of both poses were considered for the calculation using Eqs. (1), (5) and (6). The experimental free energy of binding ΔG was not available but was estimated from the experimentally determined concentration at which the protein shows 50% activity [6]. The calculated relative free energies of binding of rutinose and p-NPR agree quite well with the experimental estimates. Note that these values were obtained using the standard values of α and β in Eq. (1) which are not optimized for this protein. Furthermore, the experimental estimate of ΔΔG should be considered a rough estimate only, as p-NPR seems to show cooperative binding and will not adhere to the steady-state assumption [33].

Table 3

Calculated ΔG of rutinose and p-NPR as well as the experimentally determined concentration at 50% enzyme activity (K50). The relative binding free energy ΔΔG can be compared between the two approaches.

	Exp. K50/ΔΔG	Calc. ΔGbind
Rutinose	2.54 mM	−34.4 kJ mol−1
p-NPR	22.3 mM	−30.6 kJ mol−1
ΔΔG	5.4a kJ mol−1	3.8 kJ mol−1

Calculated as .

The occurrence of hydrogen bonds between the protein and the substrates is summarized in Table 4. The occurrence of selected hydrogen bonds for the most relevant poses and simulations (according to the weights in Table 2) is presented as a function of time in Fig. 5. Table 4 shows that the total number of hydrogen bonds is much higher for rutinose than for p-NPR. This agrees with the binding free energies of the different poses in Table 1; for the poses with low free energy of binding more hydrogen bonds to the protein are observed.

Table 4

Hydrogen bonds between the protein and the substrates. All amino acids which are involved in the hydrogen bond formation are listed. The conserved catalytic residues are highlighted. For each simulation the interactions are indicated as the percentage of simulation time in which the hydrogen bond is observed. The first column shows the interacting amino acids of Ram2 and the first row the ligands. The second row shows the various poses and the third row the independent simulation numbers. The row denoted water gives the occurrence of hydrogen bonds formed between the substrates and the water in the protein simulation. In the last row all interactions of the ligand with the protein and the water are summed up.

	Rutinose								p-NPR
	Pose 1				Pose 2				Pose 1				Pose 2
	Sim1	Sim2	Sim3	Sim4	Sim1	Sim2	Sim3	Sim4	Sim1	Sim2	Sim3	Sim4	Sim1	Sim2	Sim3	Sim4
Asp217															227
Arg221			21	35
Asp222	10					48		13		59
Trp226											36	52
Asp229	262	218	199	226			237						73
Asp277								161
Thr279								37
Tyr284				22			42									33
Val338						34
Phe339		41	40		41	35	24	12	38		25	15
Asp341		30			170	106					13	21	10
Tyr449		15
Glu488					125			82						83
Asn493														115
Ser496					64
His510					71		11	99					87			38
Ala511									13	14
Trp512							40			34
Water	1074	1068	1227	1057	906	1178	893	924	637	633	629	623	408	460	519	625
Total	1346	1372	1487	1340	1377	1401	1247	1328	688	740	703	711	578	658	746	696

Fig. 5

Hydrogen bonds formed for at least 20% of the simulation time. The graphs show rutinose from the third (A) and fourth simulation (B) in pose 1 and p-NPR from the third (C) and fourth simulation (D) in pose 2.

The interactions between the protein Ram2 and the substrates rutinose and p-NPR are shown in Fig. 6 for the most relevant poses and simulations (according to the weights in Table 2). In all four panels, the last snapshot of the 5 ns simulations is shown. Panels (A) and (B) represent pose 1 of rutinose. In (A) rutinose accepts hydrogen bonds from Arg221, and donates hydrogen bonds to Asp229 and Phe339. The hydrogen bond to the backbone of Phe339 is formed through the rhamnose moiety, while the other hydrogen bonds involve the glucose moiety. In (B) the same interactions with Arg221 and Asp229 are more persistent (Table 4 and Fig. 5). Additionally, there are hydrogen bonds formed with Tyr284 rather than with Phe339. Here, only the glucose moiety forms interactions as only the glucose is in the protein while the rhamnose is pointing outwards. Pose 2 of p-NPR is shown in (C) and (D). In (C) only one strong hydrogen bond is observed between Asp217 and rhamnose but with a very high strength of 227% (Table 4). The interaction is made between both Asp oxygens and rhamnose is acting as hydrogen donor. Since p-NPR is only kept in the protein by its rhamnose moiety the nitrophenol ring points out of the active site. In (B) and (D) hydrogen bonds are formed with Tyr284 although for (B) the interaction is with the side chain while for (D) the interaction is with the backbone. p-NPR also forms a hydrogen bond with His510.

Fig. 6

Interactions between Ram2 and rutinose or p-NPR. The interacting amino acids are shown in yellow and the conserved catalytic ones are shown in red. The snapshots are taken from the simulations with the lowest free energy. (A) Simulation 3 with pose 1 of rutinose, (B) simulation 4 with pose 1 of rutinose, (C) simulation 3 with pose 2 of p-NPR and (D) simulation 4 with pose 2 of p-NPR. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

It can be seen that in the four simulations in Fig. 6 only two of the conserved catalytic residues, Asp229 and Asp217, are involved in hydrogen bond formation, which are colored in red. However, considering all simulations all of the conserved amino acids are at some point involved in interactions as shown in Table 4 where these residues are highlighted. Several amino acids are involved in hydrogen bonds and some of them occur for more than 70% of the time. Asp217 is involved only once in an interaction but with a very high value of about 227%. In experiments in which this amino acid was replaced K remained the same while k changed [8]. Therefore, this amino acid is believed to play an important role in the hydrolysis but should be less important for the binding of the substrate. Other amino acids which are involved in many interactions during the simulations are Asp229 and Phe339. Interactions with Phe339 occur in more simulations but Asp229 makes more persistent interactions (200%), mostly with rutinose. These simulations suggest that these two amino acids could be relevant for substrate binding. The amino acids Asp341 and His510 are involved in interactions with rutinose as well as with p-NPR.

The calculation of the number of hydrogen bonds formed between the substrate and the water in the protein simulation shows that the ligands are also interacting with the surrounding water. Rutinose is involved in more interactions with water than p-NPR because it contains more hydroxyl groups. Pose 2 of p-NPR shows even less hydrogen bonds with water than p-NPR of pose 1 because the sugar ring is pointing into the protein and is therefore not available for interactions with the water but rather involved in stronger interactions with the protein. Rutinose seems to be more stable in the protein than p-NPR because it is involved in more interactions with the protein as well, as can be seen in the calculated totals (Table 4).

Several hydrophobic and aromatic amino acids are in the vicinity of the ligands in the active site for possible hydrophobic interactions such as Trp226, Tyr284, Phe339, Tyr449 and Trp512 shown in Fig. 7. The average distances are also summarized in Supplementary data 4. During the simulations of rutinose the glucose moiety was always close to Trp226. Simulations 3 and 4 of rutinose pose 1 show very stable distances with an average value of 0.49 nm and 0.46 nm, respectively. The simulations with rutinose pose 2 show less stable distances, especially seen in simulation 2 where the value increases from 0.6 nm to more than 1.0 nm. During the simulations 3 and 4 of p-NPR pose 2, Trp226 was in a stable average distance of approximately 0.56 nm.

Fig. 7

Distances between hydrophobic amino acids and the substrates (rutinose, p-NPR) are shown. The different colors show various simulations of both ligands. Only curves are drawn for which the distance is smaller than 0.6 nm at some point during the simulation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Tyr284 occurs less often. In simulation 2 of p-NPR in pose 1 it is close to the nitrophenol ring but the distance immediately increases from 0.5 nm to over 1.0 nm, as occurs in simulation 2 of rutinose pose 2 with the glucose moiety. Stable distance values are shown in simulation 4 of p-NPR and rutinose both in pose 2. The nitrophenol ring of p-NPR is in an average distance of 0.47 nm and the glucose moiety of 0.58 nm to Tyr284.

Another hydrophobic amino acid in the active site is Phe339, which shows stable average distances of 0.47 nm (p-NPR pose 2 sim 4; nitrophenol ring) and 0.56 nm (rutinose pose 2 sim 4; glucose). In other simulations Phe339 is first in an interaction distance but the distance increases to 1.0 nm or more.

Tyr449 is only close by in two simulations where it is 1.0 nm away at the start and comes closer during the simulation shown in simulation 1 of rutinose pose 1 (rhamnose) and simulation 3 of rutinose pose 2 (glucose).

Cui et al. [8] mentioned Trp512 to be in close vicinity of rhamnose, which is also shown in our simulations to be close by at a very stable distance during the different simulations. In the simulations of p-NPR pose 2 the rhamnose moiety is in an average distance of approximately 0.64 nm to Trp512. During simulation 1 and 2 of rutinose pose 1 the glucose moiety is in a stable distance of approximately 0.55 nm during the whole simulation. The rhamnose moiety of rutinose in pose 2 further shows stable distances in simulation 1 (0.6 nm), simulation 3 (0.47 nm) and simulation 4 (0.62 nm).

Discussion

p-Nitrophenyl-glycosides are widely applied to characterize novel glycosidases, as these synthetic substrates allow rapid photometric enzyme assays and are convenient tools to determine overall kinetic parameters such as the response to pH or temperature. However, performing kinetic assays exclusively with such chromogenic analogs of natural substrates do not offer the best insight, as the physiological functionality of glycosidases cannot be accessed this way. In order to achieve that, cumbersome experiments have to be conducted to uncover the specificities and affinities for natural substrates of interest. In most cases, the results have then to be accessed analytically after the assays were stopped. Consequently, analyses of stopped assays provide estimates rather than the physiological constants satisfying steady state kinetics. However, the comparison of apparent specificities and affinities obtained under identical conditions provides valuable insights into the chemistry of a given glycosidase, especially in comparison to the commonly used p-NP-substrates. These studies helped to distinguish the two rhamnosidases Ram and Ram2 from P. acidilactici as mentioned before. The reported results suggest that Ram and Ram2 may represent two subclasses within GH78 [6] which is also indicated by a sequence alignment of bacterial GH78 rhamnosidases (see Supplementary data 1). The identity of the amino acid residues serving as catalytic acid/base of GH78 rhamnosidases is not yet resolved but the identification of several key amino acid residues which may be involved in substrate recognition or processing [8] offers a better insight into this class of proteins. While these amino acids are highly conserved within GH78, the motifs surrounding these residues are apparently conserved within the proposed subclasses (see Supplementary data 1). It can be assumed that such differences determine enzyme functionalities through distinct stereochemistry at the catalytic site or by influencing substrate binding affinities. Consistent with these considerations, our experimental results indicated that Ram may be classified as a broad aryl/alkyl-rhamnosidase with high activity for p-NPR but low capacity to hydrolyze the disaccharide rutinose. Contrarily, Ram2 was almost inactive toward p-NPR but displayed remarkably high catalytic efficiency with rutinose. It is not clear yet whether this distinctiveness is a general trait in the two proposed GH78 subclasses or rather a peculiarity of the two P. acidilactici rhamnosidases, as this was the first reported case where the substrate rutinose was included in detailed kinetic studies of GH78 glycosidases. However, this unusually sharp distinction between putatively structurally similar enzymes merits a closer look into the mechanisms of GH78 rhamnosidases. Especially the behavior of Ram2 is of high interest as this is the first reported case of a true rutinase. Since rutinose is not a naturally occurring sugar (i.e., serving as potential carbon source), the metabolic function of Ram2 remains so far unclear. However, it was previously discussed that this subclass of GH78 glycosidases might be involved in the turnover of bacterial films containing rhamnose [34].

The best studied glycosidase class regarding substrate binding is that of GH1 β-glucosidases, and considering the analogous situation (aryl-glucoside/aryl-rhamnoside and glucooligosaccharide/rutinose), it is reasonable to discuss our results in this context. Although the glycone-specificity of glycosidases is mainly determined by the fold/structure of the catalytic domain (the GH family), the specific details are still not completely understood. GH1 glycosidases recognize the non-reducing glycosyl residue through a hydrogen bond network formed with several highly conserved residues (Gln, His, Glu, Asn, Trp) at the glycone binding site (subsite −1) [35,36]. Additionally, the glycone is stabilized by hydrophobic stacking interactions with the indole ring of a Trp residue [37]. Due to the high conservation of the residues interacting with the glycone, it is still not understood how a GH1 enzyme can primarily be a β-glucosidase or a β-mannosidase [1]. Further, the orientation of the aglycone in the active cleft was discussed to be of significance for glycone binding/positioning as well [35]. This adds to the fact that studying aglycone specificities is highly important to understand the nature and function of individual glycosidases. In the case of GH1 β-glucosidases, aglycone binding is apparently mostly determined by hydrophobic interactions. In contrast to the highly conserved glycone binding residues, all of the so far investigated GH1 glucosidases showed different (though sometimes “overlapping”) sets of active site residues interacting with the aglycone [1,38]. Analysis of maize GH1 β-glucosidase (ZmGlu1) showed that the (aromatic) aglycone was sandwiched between Trp378 and three Phe residues. The Trp378 residue which is highly conserved among GH1 glucosidases, is responsible for binding the aglycone (of p-nitrophenyl-β-d-thioglucoside) moiety by hydrophobic stacking [39]. It was concluded that variations of the other, less conserved sites (the Phe residues of ZmGlu1) and the resulting stereochemical changes determine the high variability in aglycone binding modes among GH1 β-glucosidases [40]. Experiments with DIMBOAGlc (2-O-β-d-glucopyranosyl-4-hydroxy-7-methoxy-1,4-benzoxazin-3-one) and cellotriose indicated that the aglycone (DIMBOA) is bound by the three Phe residues of Glu1, while as determined by docking experiments, the reducing end of cellotriose was oriented toward a distinct site. Interestingly, docking of cellotriose to the ZmGlu1 active site indicated surprisingly few (citing the authors) possibilities for likely hydrogen bonds with several amino acid side chains (Glu, Arg, Gln, Asn, Thr, His) at subsites +1 and +2. Apparently, non-polar interactions are important even when reducing end sugar(s) form the aglycone moiety of the substrate. Interestingly, Trp378 of ZmGlu1 was found to stabilize the DIMBOAGlc aglycone as well as both reducing end glucose residues of cellotriose [37].

Cui et al. [8] showed experimentally that the rhamnose of p-NPR is stabilized by Trp226 through hydrophobic stacking. This also occurs during the simulations of p-NPR with Ram2 and leads to a stronger total binding and also to a lower free energy seen in simulation 3 and 4 of pose 2 because during these simulations the distance between both centers of mass is stable during the entire simulation. In the simulations of rutinose Trp226 is parallel to the glucose but not to the rhamnose moiety. Furthermore, simulations with low free energy values and accordingly high weights display Trp226 in a stable distance during the simulations. Trp226 seems to be an important amino acid for hydrophobic interactions since it is close by in all simulations. Trp512 could also play an important role because it also is close by in most simulations and the distances fluctuate less than for Trp226. Also Cui et al. [8] suggested Trp512 to be in close approximation of the rhamnose.

Quercetin, a widely distributed flavonoid, forms several glycosides including rhamnosides (quercitrin, rutin). Tribolo et al. [38] positioned quercetin-4′-glucoside into the active site of human cytosolic GH1 glucosidase (hCBG) through homology modeling and confirmed that binding of the aglycone is mainly determined by hydrophobic interactions. No polar or charged residues close enough to form hydrogen bonds with the substrate could be identified. A similar observation was made with a flavonoid glycosyltransferase [41]. Therefore, it is important to take hydrophobic interactions into account when analyzing the substrate specificity as discussed with GH1 β-glucosidases as well as other enzymes involving glycosylated or aromatic substrates [42-44].

The aim of the present study was to investigate whether the unique characteristics of Ram2 could be rationalized through computer simulations. This should also demonstrate that MD simulations are a convenient and interesting approach to predict the substrate binding affinities of glycosidases based on a homology model. Therefore, we analyzed the binding affinities of Ram2 to p-NPR and rutinose and also evaluated which amino acid residues of Ram2 are most likely to interact with these substrates. To find the proper binding site in the protein for the ligands several approaches were used for approval. The site finder tool of MOE was used to get a first idea, this showed a site also identified by Cui et al. [8]. This finding was confirmed by a sequence alignment of several rhamnosidases which shows catalytic residues as part of the binding site, to which the substrates were subsequently docked. Finally, to verify the binding site MD simulations with the ligand in this binding site were performed to show the affinity of the ligand to the binding site. During the simulations the homology model of the protein and its complex with the substrates can adjust to a more appropriate structure, while the dynamic character of the simulation also reflects the presence of thermal fluctuations more appropriately. The calculated relative binding affinity agrees well with a rough experimental estimate. However, the experimental setup is not necessarily directly comparable to the computational one. First, the experiments are measuring the kinetics of substrates undergoing a chemical reaction, while the simulations only consider the binding process. Second, the experimental values do not necessary reflect the physiological constants in terms of steady state kinetics but possibly involve cooperativity effects. Further, in the experiments additional processes may take place such as inhibition by products and matrix, or unknown experimental errors. Such a complexity cannot be reproduced in silico, as computer power still restricts simulations to truly microscopic systems. Keeping these considerations in mind, we demonstrate that computer simulation can reproduce binding affinities of glycosidases in a comparative manner. Using Ram2 as a case study, the free energy of binding was calculated for the substrates p-NPR and rutinose, reproducing the experimental trends and allowing us to identify the favorable orientations of the ligands in the pocket. Rutinose presents its glucose moiety toward the catalytic center, in its most favored orientation. Interestingly, p-NPR prefers orientations in which the rhamnose moiety is buried in the active site, but in doing so binds with a significantly lower affinity. These observations are supported by the experimental finding that glucose is a stronger inhibitor of Ram2 than rhamnose. Several analyses were subsequently performed on the simulations to understand which amino acids might be involved in the recognition of rutinose.

As such, our simulations may be considered as a tool to suggest putative mutations to further characterize the role of these amino acids experimentally.

Conclusion

The molecular dynamics simulations indicate that rutinose binds stronger and more favorably to Ram2 than p-NPR based on the hydrogen bond analysis and free energy calculations. It was further observed that rutinose preferably orients its glucose moiety toward the active site while p-NPR shows a low free energy when it positions its rhamnose moiety in the same area. According to the hydrogen bond analysis the protein seems to prefer interactions with the glucose moiety of rutinose which agrees with the observation that glucose itself is a much stronger inhibitor than rhamnose. p-NPR seems to be a less favorable binder because only the rhamnose moiety is involved in the binding. However, the sugar group is not positioned optimally and interactions are made with different amino acids compared to rutinose. The most frequently occurring amino acids making hydrogen bonds during the rutinose simulations are Asp229 and Phe339. The most important amino acid for binding seems to be Asp229 because it is involved in many hydrogen bonds and shows a high occurrence. Also when simulated with p-NPR it is close by for an interaction although it does not bind to the nitrophenol ring.

Our simulations offer insights into this rhamnosidase at an atomic resolution. They form a complementary tool to cumbersome experimental determinations and give direct suggestions for new mutagenic experiments to characterize this intriguing enzyme.