Understanding Thermostability Factors of Barley Limit Dextrinase by Molecular Dynamics Simulations

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.

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. 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. There are several works focused on improvement of the themostability of pullulanase derived from bacteria 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 (R g ), 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.

Systems Preparation
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  were also built based on the coordinate of the corresponding residues in the free form of HvLD (PDB ID: 4AIO). A mutant of LD D317A was constructed to evaluate the effect of salt bridge between Asp317 and Arg378.

MD Simulations
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 Ca 2+ . 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 Ca 2+ at 298, 318, and 343 K were also constructed. The temperature was set as 298 and 343 K for LD D317A . 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 LD D317A 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.
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 Ca 2+ , 318 K without Ca 2+ , and 343 K without Ca 2+ are shown in black, magenta, red, cyan, deep teal, and blue, respectively.
FIGURE 3 | Radius of gyration (R g ) plot. The simulation data obtained at 298 K, 318 K, 343 K, 298 K without Ca 2+ , 318 K without Ca 2+ , and 343 K without Ca 2+ is shown in black, magenta, red, cyan, deep teal, and blue, respectively.

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 R g 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.
A one-way ANOVA was conducted here to evaluate whether the differences are significant for systems containing Ca 2+ , 318 K/298 K and 343 K/298 K, and systems without Ca 2+ (318 K/298 K and 343 K/298 K) for the RMSD, R g , and SASA. The difference is considered significant in the case of P < 0.05.

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 The distribution of total SASA. The simulation data obtained at 298, 318, 343, 298 K without Ca 2+ , 318 K without Ca 2+ , and 343 K without Ca 2+ are shown in black, magenta, red, cyan, deep teal, and blue, respectively. investigate the effect of calcium ions for the structural stability, MD simulations were also conducted for the three systems without Ca 2+ (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 Ca 2+ 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. The radius of gyration (R g ) reflects the compactness of protein structure. To detect the compactness of the overall structure, the residues. e Hydrogen bonds between hydrophilic residues. f Hydrogen bonds between hydrophobic residues. g Hydrogen bonds between hydrophilic residues and hydrophobic residues.  radius of gyration (R g ) 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 R g ranging from 25 to 100 ns in each system. The averaged R g value indicated that the compactness of HvLD increases when the temperature rises. HvLD has the most compact structure at the lowest temperature (298 K, R g : 29.95 ± 0.09 Å). Moreover, HvLD exhibits similar compactness at 318 K (R g : 29.97 ± 0.08 Å) and 343 K (R g : 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 R g values at three systems without Ca 2+ , 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 Ca 2+ was less compact than those systems with Ca 2+ . In addition, with the temperature rising, the structure of HvLD becomes less compact.
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 nm 2 , 348.90 ± 3.95 nm 2 , and 349.05 ± 3.78 nm 2 at 298 K, 318 K, and 343 K, respectively, while in the systems without Ca 2+ , the total SASA values exhibit significant increase at both temperatures. The averaged SASA is 365.80 ± 6.52 nm 2 , 357.43 ± 4.75 nm 2 , and 368.26 ± 6.79 nm 2 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 nm 2 at 298 K, with the major peak at 349 nm 2 (20.73%). The same major peak appears at 318 K (18.94%) and 343 K (22.00%), respectively. In the systems without Ca 2+ , the total SASA increases dramatically. The range of SASA is from 337 to 383 nm 2 , with the major peak at 367 nm 2 (14.91%) at 298 K. The range of SASA is from 337 to 371 nm 2 , with the major peak at 359 nm 2 (18.95%) at 318 K. The SASA at 343 K distributes at a range of 339-391 nm 2 , with the major peak at 365 nm 2 (13.31%). This profile is consistent with the trend of R g values, which indicates that HvLD become less compact with more solvent penetration into the core of the enzyme at high temperatures.
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 Ca 2+ . 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, hydrophobichydrophobic 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 hydrophobichydrophobic residues reduces at 343 K compared with that at 298 K. In the systems without Ca 2+ , the amount of the total HBs, MM HBs, MS HBs, and HBs for hydrophobichydrophobic 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 Ca 2+ , 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 Ca 2+ at 343 K. The unstable interaction would be unfavorable for the stability for the catalytic triad of Asp473-Glu510-Asp642.

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 Ca 2+ , which suggests that ions Ca 2+ contribute to the stability of salt bridges. We also analyzed the salt bridges with occupancy higher than 50%. The length of each salt bridge was averaged over 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 Ca 2+ , which are colored by white, wheat, and light blue, respectively. The high flexible region was colored by red.
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 Ca 2+ 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.
To evaluate the function of salt bridge, the salt bridge between Asp317 and Arg378 was selected. The MD simulations on LD D317A 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 LD D317A exhibited higher RMSD, R g , and SASA than LD WT at 298 K and 343 K (Supplementary Figure S1).

Hydrophobic Contacts
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 Ca 2+ .

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 Ca 2+ , 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. 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 Ca 2+ , 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 Ca 2+ 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 Ca 2+ 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 Ca 2+ , which is located between α1 and α2, had a weaker effect on the conformation of HvLD at different temperatures. This Ca 2+ contributes less to the stability of HvLD.

CONCLUSION
In the present study, thermostability factors of barley LD were investigated by MD simulations. The higher value of RMSD, R g , 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.

AUTHOR CONTRIBUTIONS
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.
FIGURE S2 | The conformation of Ca1 site. Initial structure and the extracted structure form system without Ca 2+ at 343 K are colored by cyan and purple, respectively.
TABLE S1 | Statistical analysis results of RMSD, R g and SASA for HvLD at different simulation temperatures.