Original Research ARTICLE
Understanding Thermostability Factors of Barley Limit Dextrinase by Molecular Dynamics Simulations
- 1State Key Laboratory of Biological Fermentation Engineering of Beer, Tsingtao Brewery, Qingdao, China
- 2Shandong Province Key Laboratory of Applied Mycology, College of Life Sciences, Qingdao Agricultural University, Qingdao, China
Limit dextrinase (LD) is the only endogenous starch-debranching enzyme in barley (Hordeum vulgare, Hv), which is the key factor affecting the production of a high degree of fermentation. Free LD will lose its activity in the mashing process at high temperature in beer production. However, there remains a lack of understanding on the factor affecting the themostability of HvLD at the atomic level. In this work, the molecular dynamics simulations were carried out for HvLD to explore the key factors affecting the thermal stability of LD. The higher value of root mean square deviation (RMSD), radius of gyration (Rg), and surface accessibility (SASA) suggests the instability of HvLD at high temperatures. Intra-protein hydrogen bonds and hydrogen bonds between protein and water decrease at high temperature. Long-lived hydrogen bonds, salt bridges, and hydrophobic contacts are lost at high temperature. The salt bridge interaction analysis suggests that these salt bridges are important for the thermostability of HvLD, including E568–R875, D317–R378, D803–R884, D457–R214, D468–R395, D456–R452, D399–R471, and D541–R542. Root mean square fluctuation (RMSF) analysis identified the thermal-sensitive regions of HvLD, which will facilitate enzyme engineering of HvLD for enhanced themostability.
Limit dextrinase (LD), also termed R-enzyme, pullulanase, isoamylase, or amylopectin 6-glucanohydrolase, is the only endogenous starch-debranching enzyme in barley (Hordeum vulgare, Hv) that digests amylopectin and dextrins (Manners et al., 1970; Yang et al., 2008). HvLD belongs to the glycoside hydrolase family 13 subfamily 13 (GH13_13) and can cleave α-1,6-glucosidic bonds in limit dextrins derived from amylopectin (Stam et al., 2006).
Barley is a major raw material in beer production. The major biochemical process in brewing is to degrade barley starch into fermentable sugars, which are further converted into alcohol by yeast metabolism. HvLD as a specific enzyme to digest amylopectin and dextrins is the key factor affecting the production of a high degree of fermentation (Wang et al., 2015). There are three different forms of LD existing in barley: insoluble bound, soluble inactive, and active free. Only the free form is capable for degrading amylopectin (Sissons et al., 1994). The essential industrial process of brewing includes three steps: malting, mashing, and fermentation. Mashing is usually performed at 60–70°C and at a pH of below 4.5 (Moshi et al., 2015). The heat resistance of LD in free form is poorer than the bound or latent form (Sissons et al., 1995). Free LD will lose its activity in the mashing process at a temperature higher than 63°C (Sissons et al., 1995). High themostability and activity of LD is desirable for the beer production.
Several crystal structures of HvLD have been reported (Vester-Christensen et al., 2010; Møller et al., 2012a, Møller et al., 2015a, b). The HvLD structure contains four domains (Vester-Christensen et al., 2010) (shown in Figure 1): the N-terminal domain, a carbohydrate binding module 48 (CBM48), a catalytic (β/α) 8 domain containing the two catalytic residues (Asp473, nucleophile; Glu510, general acid/base) and the transition-state stabilizer (Asp642), and a C-terminal domain. The N-terminal domain includes residues 2–124 resembling carbohydrate binding module 21. CBM48 includes residues 125–230. The catalytic domain contains residues 231–774 and the C-terminal domain contains residues 775–885.
Figure 1. The overall structure of Barley limit dextrinase. N-domain, Magenta; CBM48, cyan; catalytic domain, wheat; C-domain, light green; Ca2+ ions, red spheres. The catalytic residues (Asp473, Glu510, and Asp642) are represented in sticks.
There are several works focused on improvement of the themostability of pullulanase derived from bacteria (Chen et al., 2015; Li et al., 2015; Chang et al., 2016; Wang et al., 2016). However, it is still lacking the understanding on the factor affecting the themostability of HvLD at the atomic level. It is suggested that enzymes keep their structural stability by various kinds of non-covalent interactions, such as hydrogen bonds, salt bridges, disulfide bonds, and hydrophobic interaction (Nick Pace et al., 2014; Nilofer et al., 2017). Recently, molecular dynamics (MD) simulation, as a useful tool, has been widely applied to find important characteristics of protein stability (Alizadeh-Rahrovi et al., 2015; Sharma and Sastry, 2015; Jiang et al., 2016; Idrees et al., 2017; Gu et al., 2019).
In this work, MD simulations were carried out for barley limit dextrinase (HvLD) to explore the key factors affecting the thermal stability of LD. The root mean square deviation (RMSD), radius of gyration (Rg), and surface accessibility (SASA) were calculated to explore the dynamics of HvLD. Intra-protein hydrogen bonds, protein–water hydrogen bonds, salt bridges, and hydrophobic interaction were analyzed to find the factors about the thermal stability of HvLD. Finally, root mean square fluctuation (RMSF) analysis was performed to identify the thermal-sensitive regions of LD. The structural and dynamic details will help to understand the driving forces that lead to the stability of HvLD at different temperatures, which will facilitate enzyme engineering of HvLD.
Materials and Methods
The X-ray structure of barley LD (PDB ID: 4CVW) (Møller et al., 2015a) was obtained from the RCSB Protein Data Bank. The LD inhibitor was removed from this structure. The structure of barley LD (PDB ID: 4CVW) was superimposed on the free form of HvLD (PDB ID: 4AIO) (Møller et al., 2012b). The missing residues (43-PSN-45, 102-FGADGK-107) were also built based on the coordinate of the corresponding residues in the free form of HvLD (PDB ID: 4AIO). A mutant of LDD317A was constructed to evaluate the effect of salt bridge between Asp317 and Arg378.
Standard AMBER ff03 force field (Wang et al., 2004; Hornak et al., 2006) was assigned to the protein. The force field parameter developed by Bradbrook et al. (1998) was assigned for the Ca2+. The protonation state of ionizable residues was set under pH 5.5 based on the pKa values calculated by the H++ server (Anandakrishnan et al., 2012). Na+ ions were added to neutralize the overall system. Each system was embedded in a rectangular box of the TIP3P water molecule (Jorgensen et al., 1983), maintaining a distance of 10 Å from any solute atom to the boundary.
The MD simulations were performed using AMBER12. Energy minimization was carried out with a decreasing harmonic force constraint on the protein. The minimized system was gradually heated from 0 K to the desired temperature within 200 ps under the NVT ensemble condition. The temperature was set as 298 K, 318 K (optimum temperature), and 343 K (the highest mashing temperature), respectively. To investigate the effect of calcium ions for the structural stability, three systems without Ca2+ at 298, 318, and 343 K were also constructed. The temperature was set as 298 and 343 K for LDD317A. Then, the system was relaxed within 1.55 ns under the NPT ensemble condition. Finally, a total of 100 ns was simulated to produce trajectories under the NPT ensemble condition for each system. A 50-ns MD simulation was conducted for LDD317A at both temperatures. The covalent bonds to hydrogen atoms were constrained using the SHAKE algorithm (Coleman et al., 1977) and the Particle Mesh Ewald (PME) method (Darden et al., 1993) was employed to calculate long-range electrostatic interactions. The real space cutoff was set at 10.0 Å, the same as that for van der Waals interactions. The grid-spacing and convergence criteria of PME calculation was set to 1 Å and 1.0E-05, respectively. The time step used for the simulations was set to 2 fs. The atom coordinates were saved every 10 ps for subsequent analysis.
All of the analyses were performed using the analytical tools cpptraj (Roe and Cheatham, 2013) module in AMBER tools and VMD (Humphrey et al., 1996). The RMSD was calculated as the deviation between backbone atoms of the protein with respect to the initial structure’s backbone atoms, averaged over the backbone atoms. For the RMSF and SASA, the average value was calculated over time. SASA was calculated for all of the residues, hydrophobic residues, hydrophilic residues, and catalytic residues (Asp473, Glu510, and Asp642), respectively. The Rg was calculated by VMD.
The hydrogen bond was calculated based on a maximum cutoff distance between the donor and the acceptor at 3.5 Å and the angle of donor–hydrogen–acceptor larger than 120°. The average value of the number of hydrogen bond was calculated as the ratio of the sum of the total number of HBs in each frame to the total number of frames. The redundant hydrogen bonds between the same donor and acceptor but with different hydrogen atoms were removed, saving the one with the highest occupancy. The HBs were analyzed by considering the chemical properties of different residues, including charged residues (Arg, Lys, Asp, and Glu), polar residues (Gln, Asn, Ser, Thr, Tyr, and Cys), and hydrophobic residues (Ala, Ile, Leu, Phe, Val, Pro, Gly, Met, and Trp).
The salt bridges (SBs) were considered to be formed if the distance between an oxygen atom of an acidic residue (Oδ1 and Oδ2 of ASP and Oδ1 and Oδ2 of Glu) and the nitrogen atom of a basic residue (Nε, Nη1 and Nη2 of Arg and Nζ of Lys) was less than 4 Å. The average value of the number of SBs was calculated as the ratio of the sum of the total number of SBs in each frame to the total number of frames. The SBs between the same two residues but different atoms were regarded as unique, keeping the one with the highest occupancy.
It is suggested that all C atoms within 3.9 Å interacts through hydrophobic contacts (Stojanovic and Zaric, 2009). We calculated the hydrophobic contacts between all the hydrophobic atoms (C, Cα, Cβ, Cδ, Cδ1, Cδ2, Cε, Cε1, Cε2, Cε3, Cγ, Cγ1, Cγ2, Cζ, Cζ2, Cζ3, and Cη2) with a cutoff of 4 Å, without redundancy. The structures were visualized by VMD and PyMOL (Schrödinger, 2010).
A one-way ANOVA was conducted here to evaluate whether the differences are significant for systems containing Ca2+, 318 K/298 K and 343 K/298 K, and systems without Ca2+ (318 K/298 K and 343 K/298 K) for the RMSD, Rg, and SASA. The difference is considered significant in the case of P < 0.05.
Results and Discussion
Dynamics of Barley LD
In order to identify the key factors responsible for instability of HvLD at high temperature, MD simulations were performed at different temperature conditions (298, 318, and 343 K) to predict the molecular behavior over the period of time using AMBER. To investigate the effect of calcium ions for the structural stability, MD simulations were also conducted for the three systems without Ca2+ (298, 318, and 343 K). The structural stability of six systems was examined by calculation of the RMSD of the backbone atoms relative to the initial structure. Figure 2 shows the RMSD variations of the six systems with respect to simulation time. It is observed that the systems are equilibrated and thus suitable for exploring the dynamics of HvLD. The first three systems achieved equilibrium at 20, 10, and 25 ns, respectively. The RMSD values for the backbone atoms of HvLD converge at 1.40 ± 0.11 and 1.39 ± 0.10 Å at 298 and 318 K and 1.59 ± 0.16 Å at 343 K. The statistical analysis suggests that the differences are significant (P < 0.05) (Supplementary Table S1). HvLD in the systems without Ca2+ exhibits large variation compared with the initial structure at 298 and 343 K. The backbone RMSD increases rapidly and major structural distortion occurs at 343 K.
Figure 2. Time evolution of the backbone RMSD versus the starting structure of HvLD. The simulation data obtained at 298, 318, 343, 298 K, without Ca2+, 318 K without Ca2+, and 343 K without Ca2+ are shown in black, magenta, red, cyan, deep teal, and blue, respectively.
The radius of gyration (Rg) reflects the compactness of protein structure. To detect the compactness of the overall structure, the radius of gyration (Rg) of the protein in six systems was also calculated, as shown in Figure 3. According to RMSD plot, all the systems converged after 25 ns. We calculated the average Rg ranging from 25 to 100 ns in each system. The averaged Rg value indicated that the compactness of HvLD increases when the temperature rises. HvLD has the most compact structure at the lowest temperature (298 K, Rg: 29.95 ± 0.09 Å). Moreover, HvLD exhibits similar compactness at 318 K (Rg: 29.97 ± 0.08 Å) and 343 K (Rg: 30.14 ± 0.09 Å). The structure of HvLD is less compact at the higher temperature than it at the lower temperature, indicating expansion of protein structure at higher temperature. The HvLD exhibited higher Rg values at three systems without Ca2+, with the value of 30.69 ± 0.15 Å (298 K), 30.47 ± 0.10 (318 K), and 31.13 ± 0.34 Å (343 K), respectively. This result indicated that the structure of protein in the system without Ca2+ was less compact than those systems with Ca2+. In addition, with the temperature rising, the structure of HvLD becomes less compact.
Figure 3. Radius of gyration (Rg) plot. The simulation data obtained at 298 K, 318 K, 343 K, 298 K without Ca2+, 318 K without Ca2+, and 343 K without Ca2+ is shown in black, magenta, red, cyan, deep teal, and blue, respectively.
To evaluate the exposure of protein atoms to solvent, SASA was also obtained by calculating the surface area of atom in contact with solvent molecules. From Figure 4A, it is found that the total SASA values show a slight increase with a rise of temperature. The averaged SASA is 347.63 ± 4.29 nm2, 348.90 ± 3.95 nm2, and 349.05 ± 3.78 nm2 at 298 K, 318 K, and 343 K, respectively, while in the systems without Ca2+, the total SASA values exhibit significant increase at both temperatures. The averaged SASA is 365.80 ± 6.52 nm2, 357.43 ± 4.75 nm2, and 368.26 ± 6.79 nm2 at 298 K, 318 K, and 343 K, respectively. Total SASA, SASA of hydrophilic residues and catalytic residues, increases relatively from 298 K to 343 K (Table 1). The distribution of SASA is displayed in Figure 4B. It can be observed that the total SASA is from 331 to 359 nm2 at 298 K, with the major peak at 349 nm2 (20.73%). The same major peak appears at 318 K (18.94%) and 343 K (22.00%), respectively. In the systems without Ca2+, the total SASA increases dramatically. The range of SASA is from 337 to 383 nm2, with the major peak at 367 nm2 (14.91%) at 298 K. The range of SASA is from 337 to 371 nm2, with the major peak at 359 nm2 (18.95%) at 318 K. The SASA at 343 K distributes at a range of 339–391 nm2, with the major peak at 365 nm2 (13.31%). This profile is consistent with the trend of Rg values, which indicates that HvLD become less compact with more solvent penetration into the core of the enzyme at high temperatures.
Figure 4. Total solvent-accessible surface area (SASA). (A) Time evolution of total SASA. (B) The distribution of total SASA. The simulation data obtained at 298, 318, 343, 298 K without Ca2+, 318 K without Ca2+, and 343 K without Ca2+ are shown in black, magenta, red, cyan, deep teal, and blue, respectively.
The hydrogen bond interaction is considered important in protein folding, stability, and function. It can be seen that HvLD lost 55 hydrogen bonds at 343 K with respect to 298 K (Table 2). The number of hydrogen bonds also decreases in the systems without Ca2+. Based on the percentage existence time of HBs (Figure 5), short-lived HBs (0 < X ≤ 10%) increase significantly at 343 K, indicating that most of them appear transiently at high temperature. In contrast, the number of long-lived HBs decreases at 343 K, indicating that these interactions are unable to maintain at high temperature. In addition, the number of substantially live HBs (10 < X ≤ 90%) at 343 K is more than the corresponding value at 298 K, indicating that these HBs contribute to the stability of HvLD. We analyzed the HBs with occupancy higher than 50% at 298, 318, and 343 K. There are 56 HBs becoming weak and their occupancy decreases with the increase in temperature (Supplementary Table S2). These HBs are very sensitive to temperature changes and they affect the thermostability of HvLD.
Besides, we also analyzed the effect of temperature on main chain–main chain HBs, main chain–side chain HBs, and side chain–side chain HBs. There are more MM HBs compared with two other types. There are 18 MM HBs broken at 343 K compared with those at 298 K. MM HBs are important in secondary structure formation. The decrease of the number of MM HBs indicates that the stability of the secondary structure of the enzyme would be impaired. MS HBs and SS HBs exhibit slight decrease at 343 K with respect to 298 K. Moreover, HBs among the residues having similar/different chemical properties (charged–charged residues, hydrophobic–hydrophobic residues, hydrophilic–hydrophilic residues, and hydrophobic–hydrophilic residues) were analyzed. Among these, the number of HBs between hydrophilic–hydrophilic residues and hydrophobic–hydrophilic residues does not exhibit a significant difference at different temperatures, while the amount of HBs between charged–charged residues and hydrophobic–hydrophobic residues reduces at 343 K compared with that at 298 K. In the systems without Ca2+, the amount of the total HBs, MM HBs, MS HBs, and HBs for hydrophobic–hydrophobic residues, hydrophilic–hydrophilic residues, and hydrophobic–hydrophilic residues also decreases at 343 K compared with 298 K. The amount of SS HBs and HBs for charged–charged residues increases at 343 K compared with 298 K, indicating that there are new hydrogen bonds formed in the distorted structure.
The changes of HBs between protein and water were also explored in the temperature range of 298–343 K, which is displayed in Table 3. There is a decrease in protein–water HBs from 298 to 343 K, due to the loss of both M-water HBs and S-water HBs from 298 to 343 K. Besides, the number of HBs for hydrophobic residues to water, hydrophilic residues to water, and charged residues to water reduces at higher temperatures, indicating that the network of HBs between HvLD and water molecules is broken. In the systems without Ca2+, the amount of total protein–water HBs and other types of protein–water HBs decreases at 343 K compared with 298 K, indicating the large change of the network of HBs between HvLD and water molecules in these two systems.
The Stability of the Hydrogen Bond Between D403 and H404
Structural stability may also affect the catalytic activity of this enzyme. The hydrogen bond between D403 and H404 is favorable for the stability of the catalytic triad, which is suggested in a previous work (Møller et al., 2015b). The distance between these two residues was monitored during the whole trajectory. It can be seen that the interaction is more stable at 298 and 318 K (Figure 6). When the temperature rises, the interaction becomes unstable (at 343 K). In addition, most of this interaction disappeared in the system without Ca2+ at 343 K. The unstable interaction would be unfavorable for the stability for the catalytic triad of Asp473–Glu510–Asp642.
Figure 6. Monitoring the hydrogen bonds between the OD1 atom of Asp403 and the ND1 atom of His404 during MD simulation. (A) 298 K; (B) 318 K; (C) 343 K; (D) 298 K without Ca2+, (E) 318 K without Ca2+, (F) 343 K without Ca2+.
Salt Bridge Interaction Analysis
Besides hydrogen bonds, salt bridges are also very important for the stability of protein (Horovitz et al., 1990; Strop and Mayo, 2000; Jelesarov and Karshikoff, 2009). Salt bridges in enzyme may contribute to its stability at high temperature (Vieille and Zeikus, 2001; Kundu and Roy, 2010). Recently, Guo et al. (2018) summarized factors may contribute to the thermostability for pullulan-hydrolyzing enzymes. They found that there are more salt bridges in thermophilic pullulan-hydrolyzing enzymes than mesophilic ones, suggesting the importance of salt bridge for the enzyme pullulan-hydrolyzing thermostability. In this work, the salt bridges were identified using a 4-Å distance cutoff. We can observe that the number of short-lived salt bridges increases from 298 to 343 K (Figure 7), while the number of long-lived salt bridges decreases from 298 to 343 K. It is indicated that the short-lived salt bridges form transiently at high temperature. Some long-lived salt bridges are broken at high temperature. The number of long-lived salt bridges also decreases in the systems without Ca2+, which suggests that ions Ca2+ contribute to the stability of salt bridges.
Figure 7. The number of unique salt bridges based on the percentage existence time at cutoff of 4 Å.
We also analyzed the salt bridges with occupancy higher than 50%. The length of each salt bridge was averaged over the whole trajectory. There are 33, 30, and 29 salt bridges with occupancy higher than 50% during the simulations at 298, 318, and 343 K and 19 and 16 salt bridges in two systems without Ca2+ at 298 and 343 K. The salt bridges are noted to decrease in number with an increase in temperature. There are 21 common salt bridges identified at 298, 318, and 343 K. Among these salt bridges, five could still be maintained well even at high temperatures, suggesting an essential role in stabilizing this enzyme. There are eight salt bridges exhibiting comparable occupancy at three temperatures. These salt bridges were also not affected by increases in temperature. These 13 salt bridges are correlated with the partial structural stability of HvLD, while 8 salt bridges become weak and their occupancy decreases with an increase in temperature (Table 4), including E568–R875, D317–R378, D803–R884, D457–R214, D468–R395, D456–R452, D399–R471, and D541–R542. These salt bridges are very sensitive to temperature changes and they affect the thermostability of HvLD. All of the salt bridges are located on the surface of HvLD.
Table 4. Average length and occupancy of important salt bridge interactions of HvLD at 298, 318, and 343 K.
To evaluate the function of salt bridge, the salt bridge between Asp317 and Arg378 was selected. The MD simulations on LDD317A at 298 and 343 K were conducted to evaluate the effect of salt bridge D317–R378 on the thermal stability of HvLD. It is indicated that the break of this salt bridge destabilizes the protein. The LDD317A exhibited higher RMSD, Rg, and SASA than LDWT at 298 K and 343 K (Supplementary Figure S1).
Protein stabilization is maintained by various interactions. Hydrophobic interaction is one of the important part of them (Van Dan Burg et al., 1994). Unique hydrophobic contacts were calculated as described in the “Materials and Methods” section. The percentage of time is displayed in Figure 8. It can be seen that the short-lived unique hydrophobic contacts (0 < X ≤ 10%) increase from 298 to 343 K. Long-lived hydrophobic contacts (90 < X ≤ 100%) decrease from 298 to 343 K. It exhibits the same trend with hydrogen bonds and salt bridges. A similar situation also occurs in systems without Ca2+.
Figure 8. The number of unique hydrophobic contacts based on the percentage existence time at cutoff of 4 Å.
Identification of Thermal-Sensitive Regions
To explore the structural and dynamic changes, we identified the thermal-sensitive regions of HvLD by analyzing the structural mobility based on the RMSF of the backbone atoms with respect to the initial structure. Figure 9 shows that the RMSF values of most regions of HvLD fluctuate slightly at a temperature of up to 343 K, suggesting that these regions are relatively thermostable. Some regions showed steep RMSF fluctuations at high temperature, such as 318 and 343 K, indicating that those are thermal-sensitive regions. It can be observed in Figure 9 that the highest fluctuations occur at the N-terminal because it was not restrained. Also, some loops that always exhibit high fluctuation at different temperatures may be due to their intrinsic flexibility, including residues 102–107, 135–138, and 575–582. Besides, regions that fluctuated higher than 0.5 and 1.0 Å are highlighted by magenta and red in Figure 9. Larger fluctuations are observed at 318 K or 343 K, such as residues 21–29, 42–46, 341–349, and 410–415. These regions are predicted to be thermal-sensitive regions. In a previous work, it was suggested that residues 23–27, 42–48, 102–109, and 806–810 exhibited high flexibility, with low level electron density (Vester-Christensen et al., 2010), which is consistent with our results. To compare with those systems without Ca2+, it can be observed that the flexibility of LD is all higher in these systems than that at 298 K. In particular, these loops exhibit high fluctuations, including residues 21–29, 42–46, 322–350, 410–415, 428–442, 550–557, 575–582, 720–740, and 805–812.
Figure 9. RMSF plots of the backbone atoms of six systems. (A) 298 and 318 K. (B) 298 and 343 K. (C) Simulation at 298 and 318 K without Ca2+. (D) Simulation at 298 and 343 K without Ca2+.
To exhibit these regions in HvLD, the higher fluctuation regions were mapped onto the tertiary structure. As shown in Figure 10, it is observed that some of them are located near the catalytic crevice, including residues 550–557 (Figure 10A). The high flexibility of these regions would decrease the stability of the catalytic triad. Residues 21–29, 42–46, 341–349, and 410–415 are located in surface loops (Figure 10A). In the systems without Ca2+, there are large conformational change of residues 322–350, 428–442, 720–740, and 805–812 at high temperature (343 K) (Figure 10B). The loss of Ca2+ made the significant conformational change of a long loop in catalytic (β/α) eight domain, previously described as loop 2 (Jespersen et al., 1991; Mikami et al., 2006; Turkenburg et al., 2009) (the site for Ca1) (Supplementary Figure S2), indicating that this Ca2+ are important for maintaining the stable conformation of HvLD at high temperature. The conformational change of loop 2 affected the conformation of residue 322–350 and 720–740. The removal of another Ca2+, which is located between α1 and α2, had a weaker effect on the conformation of HvLD at different temperatures. This Ca2+ contributes less to the stability of HvLD.
Figure 10. Snapshots of HvLD structures from MD trajectories. (A) Comparison of the structure extracted from the system at 298, 318, and 343 K, which are colored by white, wheat, and light blue, respectively. (B) Comparison of the structure extracted from the system at 298, 298, and 343 K without Ca2+, which are colored by white, wheat, and light blue, respectively. The high flexible region was colored by red.
In the present study, thermostability factors of barley LD were investigated by MD simulations. The higher value of RMSD, Rg, and SASA suggests the instability of HvLD at high temperatures. Intra-protein hydrogen bonds and hydrogen bonds between protein and water decrease at high temperature. Long-lived hydrogen bonds, salt bridges, and hydrophobic contacts are lost at high temperature. The salt bridge interaction analysis suggests that these salt bridges are important for the thermostability of HvLD, including E568–R875, D317–R378, D803–R884, D457–R214, D468–R395, D456–R452, D399–R471, and D541–R542, which are located on the surface of HvLD. Based on RMSF calculations for HvLD at various temperatures, several thermally sensitive regions of HvLD were identified, such as residues 21–29, 42–46, 341–349, and 410–415. The structural and dynamic details will help to understand the driving forces that lead to the stability of HvLD at different temperatures, which will facilitate enzyme engineering of HvLD for enhanced thermostability.
Data Availability Statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.
JDu conceived and supervised the experiments. JDo, SD, and KZ performed the MD simulations. JDo, JDu, JY, and SH analyzed the data. JDu and HY wrote the manuscript.
This work was supported by the Open Research Fund of State Key Laboratory of Biological Fermentation Engineering of Beer (grant no. K201804), the Natural Science Foundation of Shandong Province, China (grant no. ZR2018QB004) and the Off-site collaborative calculation of marine medicinal biological resources (grant no. 2018ASKJ01-05).
Conflict of Interest
JDu, JDo, JY, SH, and HY were employed by company Tsingtao Brewery.
The remaining 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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2020.00051/full#supplementary-material
FIGURE S1 | RMSD, Rg, and SASA for HvLDD317A. (A) RMSD plot. (B) Rg plot. (C) SASA plot.
FIGURE S2 | The conformation of Ca1 site. Initial structure and the extracted structure form system without Ca2+ at 343 K are colored by cyan and purple, respectively.
TABLE S1 | Statistical analysis results of RMSD, Rg and SASA for HvLD at different simulation temperatures.
TABLE S2 | Length and occupancy of important hydrogen bond interactions of HvLD at 298, 318, and 343 K.
Alizadeh-Rahrovi, J., Shayesteh, A., and Ebrahim-Habibi, A. (2015). Structural stability of myoglobin and glycomyoglobin: a comparative molecular dynamics simulation study. J. Biol. Phys. 41, 349–366. doi: 10.1007/s10867-015-9383-2
Anandakrishnan, R., Aguilar, B., and Onufriev, A. V. (2012). H++ 3.0: automating pK prediction and the preparation of biomolecular structures for atomistic molecular modeling and simulations. Nucleic Acids Res. 40, W537–W541. doi: 10.1093/nar/gks375
Bradbrook, G. M., Gleichmann, T., Harrop, S. J., Habash, J., Raftery, J., Kalb, J., et al. (1998). X-Ray and molecular dynamics studies of concanavalin-A glucoside and mannoside complexes Relating structure to thermodynamics of binding. J. Chem. Soc. Faraday Trans. 94, 1603–1611. doi: 10.1039/A800429C
Chang, M., Chu, X., Lv, J., Li, Q., Tian, J., and Wu, N. (2016). Improving the thermostability of acidic pullulanase from bacillus naganoensis by rational design. PLoS One 11:e0165006. doi: 10.1371/journal.pone.0165006
Chen, A., Li, Y., Nie, J., McNeil, B., Jeffrey, L., Yang, Y., et al. (2015). Protein engineering of Bacillus acidopullulyticus pullulanase for enhanced thermostability using in silico data driven rational design methods. Enzyme Microb. Technol. 78, 74–83. doi: 10.1016/j.enzmictec.2015.06.013
Coleman, T. G., Mesick, H. C., and Darby, R. L. (1977). Numerical integration: a method for improving solution stability in models of the circulation. Ann. Biomed. Eng. 5, 322–328. doi: 10.1007/BF02367312
Guo, J., Coker, A. R., Wood, S. P., Cooper, J. B., Keegan, R. M., Ahmad, N., et al. (2018). Structure and function of the type III pullulan hydrolase from Thermococcus kodakarensis. Acta Crystallogr. Sec. D 74, 305–314. doi: 10.1107/S2059798318001754
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
Horovitz, A., Serrano, L., Avron, B., Bycroft, M., and Fersht, A. R. (1990). Strength and co-operativity of contributions of surface salt bridges to protein stability. J. Mol. Biol. 216, 1031–1044. doi: 10.1016/S0022-2836(99)80018-7
Idrees, D., Rahman, S., Shahbaaz, M., Haque, M. A., Islam, A., Ahmad, F., et al. (2017). Estimation of thermodynamic stability of human carbonic anhydrase IX from urea-induced denaturation and MD simulation studies. Int. J. Biol. Macromol. 105, 183–189. doi: 10.1016/j.ijbiomac.2017.07.010
Jelesarov, I., and Karshikoff, A. (2009). “Defining the role of salt bridges in protein stability,” in Protein Structure, Stability, and Interactions, ed. J. W. Shriver (Totowa, NJ: Humana Press), doi: 10.1007/978-1-59745-367-7_10
Jespersen, H. M., MacGregor, E. A., Sierks, M. R., and Svensson, B. (1991). Comparison of the domain-level organization of starch hydrolases and related enzymes. Biochem. J. 280, 51–55. doi: 10.1042/bj2800051
Jiang, X., Chen, G., and Wang, L. (2016). Structural and dynamic evolution of the amphipathic N-terminus diversifies enzyme thermostability in the glycoside hydrolase family 12. Phys. Chem. Chem. Phys. 18, 21340–21350. doi: 10.1039/C6CP02998A
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
Kundu, S., and Roy, D. (2010). Structural study of carboxylesterase from hyperthermophilic bacteria Geobacillus stearothermophilus by molecular dynamics simulation. J. Mol. Graph. Model. 28, 820–827. doi: 10.1016/j.jmgm.2010.03.001
Li, S., Xu, J., Bao, Y., Zheng, H., and Song, H. (2015). Structure and sequence analysis-based engineering of pullulanase from Anoxybacillus sp. LM18-11 for improved thermostability. J. Biotechnol. 210, 8–14. doi: 10.1016/j.jbiotec.2015.06.406
Mikami, B., Iwamoto, H., Malle, D., Yoon, H.-J., Demirkan-Sarikaya, E., Mezaki, Y., et al. (2006). Crystal structure of pullulanase: evidence for parallel binding of oligosaccharides in the active site. J. Mol. Biol. 359, 690–707. doi: 10.1016/j.jmb.2006.03.058
Møller, M., Hachem, M. Abou, Svensson, B., and Henriksen, A. (2012a). Structure of the starch-debranching enzyme barley limit dextrinase reveals homology of the N-terminal domain to CBM21. Acta Crystallogr. Sect. F Struct. Biol. Cryst. Commun. 68(Pt 9), 1008–1012. doi: 10.1107/S1744309112031004
Møller, M. S., Abou Hachem, M., Svensson, B., and Henriksen, A. (2012b). Structure of the starch-debranching enzyme barley limit dextrinase reveals homology of the N-terminal domain to CBM21. Acta Crystallogr. 68(Pt 9), 1008–1012.
Møller, M. S., Vester-Christensen, M. B., Jensen, J. M., Hachem, M. A., Henriksen, A., and Svensson, B. (2015a). Crystal Structure of Barley Limit Dextrinase-Limit Dextrinase Inhibitor (LD-LDI) Complex Reveals Insights into Mechanism and Diversity of Cereal Type Inhibitors. J. Biol. Chem. 290, 12614–12629. doi: 10.1074/jbc.M115.642777
Møller, M. S., Windahl, M. S., Sim, L., Bøjstrup, M., Abou Hachem, M., Hindsgaul, O., et al. (2015b). Oligosaccharide and Substrate Binding in the Starch Debranching Enzyme Barley Limit Dextrinase. J. Mol. Biol. 427(6 Pt B), 1263–1277. doi: 10.1016/j.jmb.2014.12.019
Moshi, A. P., Hosea, K. M. M., Elisante, E., Mamo, G., and Mattiasson, B. (2015). High temperature simultaneous saccharification and fermentation of starch from inedible wild cassava (Manihot glaziovii) to bioethanol using Caloramator boliviensis. Bioresour. Technol. 180, 128–136. doi: 10.1016/j.biortech.2014.12.087
Nilofer, C., Sukhwal, A., Mohanapriya, A., and Kangueane, P. (2017). Protein-protein interfaces are vdW dominant with selective H-bonds and (or) electrostatics towards broad functional specificity. Bioinformation 13, 164–173. doi: 10.6026/97320630013164
Sharma, R., and Sastry, G. N. (2015). Deciphering the dynamics of non-covalent interactions affecting thermal stability of a protein: molecular dynamics study on point mutant of thermus thermophilus isopropylmalate dehydrogenase. PLoS One 10:e0144294. doi: 10.1371/journal.pone.0144294
Stam, M. R., Danchin, E. G. J., Rancurel, C., Coutinho, P. M., and Henrissat, B. (2006). Dividing the large glycoside hydrolase family 13 into subfamilies: towards improved functional annotations of α-amylase-related proteins. Protein Eng. Design Select. 19, 555–562. doi: 10.1093/protein/gzl044
Stojanovic, S. D., and Zaric, S. D. (2009). Hydrogen bonds and hydrophobic interactions of porphyrins in porphyrin-containing proteins. Open Struct. Biol. J. 3, 34–41. doi: 10.2174/1874199100903010034
Turkenburg, J. P., Brzozowski, A. M., Svendsen, A., Borchert, T. V., Davies, G. J., and Wilson, K. S. (2009). Structure of a pullulanase from Bacillus acidopullulyticus. Proteins 76, 516–519. doi: 10.1002/prot.22416
Van Dan Burg, B., Dijkstra, B. W., Vriend, G., Van Dar Vinne, B., Venema, G., and Eijsink, V. G. H. (1994). Protein stabilization by hydrophobic interactions at the surface. Eur. J. Biochem. 220, 981–985. doi: 10.1111/j.1432-1033.1994.tb18702.x
Vester-Christensen, M. B., Hachem, M. A., Svensson, B., and Henriksen, A. (2010). Crystal structure of an essential enzyme in seed starch degradation: barley limit dextrinase in complex with cyclodextrins. J. Mol. Biol. 403, 739–750. doi: 10.1016/j.jmb.2010.09.031
Wang, X., Nie, Y., Mu, X., Xu, Y., and Xiao, R. (2016). Disorder prediction-based construct optimization improves activity and catalytic efficiency of Bacillus naganoensis pullulanase. Sci. Rep. 6:24574. doi: 10.1038/srep24574
Wang, X., Zhang, X., Cai, S., Ye, L., Zhou, M., Chen, Z., et al. (2015). Genetic diversity and QTL mapping of thermostability of limit dextrinase in barley. J. Agric. Food Chem. 63, 3778–3783. doi: 10.1021/acs.jafc.5b00190
Keywords: barley limit dextrinase, thermostability, molecular dynamics simulation, hydrogen bond, salt-bridge
Citation: Du J, Dong J, Du S, Zhang K, Yu J, Hu S and Yin H (2020) Understanding Thermostability Factors of Barley Limit Dextrinase by Molecular Dynamics Simulations. Front. Mol. Biosci. 7:51. doi: 10.3389/fmolb.2020.00051
Received: 02 December 2019; Accepted: 16 March 2020;
Published: 16 April 2020.
Edited by:Guang Hu, Soochow University, China
Reviewed by:Jian Tian, Chinese Academy of Agricultural Sciences, China
Hao Dong, Nanjing University, China
Copyright © 2020 Du, Dong, Du, Zhang, Yu, Hu and Yin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hua Yin, email@example.com