Chitosanases, enzymes that catalyze the endo-hydrolysis of glycolytic links in chitosan, belong to glycoside hydrolase families. They are useful for producing chit oligosaccharides a GlcN using in the food and pharmacological industries [ 1 , 2 ]. Chitosanases EAG1 is a classical glycoside hydrolase from Bacillus ehimensis. It can hydrolyze different sizes of β 1,4-glycosidic bonds of chitosan and it can also keep high catalytic activity in organic solvents with metal ions. Nevertheless, EAG1 is relatively unstable, especially when the temperature above 50˚C its activity decreases abruptly [ 3 ]. In industrial enzyme applications, thermo stability is a key factor. Thermo stable chitosanases active between 60˚C and 100˚C is of special interest [ 4 , 5 ]. In order to improve the thermo stability of proteins, a number of protein engineering techniques have been developed, such as random mutagenesis, DNA shuffling, truncation, circularization [ 6 ] and so on [ 7 , 8 ]. However, because these methods are often required to create large enzyme mutant libraries for positive clones screening, it is difficult to acquire beneficial mutants. Thermo stability is presumably based on the protein structure. From protein structure to clarify its structural basis of thermo stable character, and then guide site directed mutagenesis will great facilitate to improve the thermo stability. In the previous research, sites in flexible regions were mutated to cysteine and these mutation were proved to be effective in improving thermo stability and catalytic efficiency of EAG1 [ 9 ]. Are there any other residues that can also play this role?
In the present study, we focused on the investigation the theoretical basis for thermo stability of EAG1 by using the structural modeling and molecular dynamical simulation. The modeled structure will give us a better understanding of the relationship between protein structure and its function. Meanwhile, the modeled structure of EAG1 will be of importance in the screening and designing of its inhibitor. This structure may also facilitate its industrial applications. Meanwhile, using the molecular dynamic simulation, we predicted the flexible resides for EAG1. Our results showed that residues Leu84, Gly113, Asp116, Ala207 and Leu286 are of important for the stability of EAG1. This result will give information on sites mutant to improve the thermostability of EAG1.
2. Material and method
2.1. Homology Modeling
Homology Modeling was used to construct the tertiary structure of EAG1. The target sequence of EAG1 (GenBank accession number AB008788) was downloaded from the National Center for Biotechnology Information (NCBI) protein database. Meanwhile, a sequence similarity search for this protease against sequences from the Protein Data Bank (PDB) database was performed on the BLAST online server (http://blast.ncbi.nlm.nih.gov). The tertiary structure of EAG1 was built using Rosetta modeling [ 10 , 11 ], which was used for homology modeling of protein three-dimensional structures. The predicted structures with the highest Rosetta score was selected for further refinement. The refinement is taken by all-atom molecular dynamics (MD) simulations in Amber16 [ 12 ]. To assess the quality of the optimized models, PROCHECK (http://servicesn.mbi.ucla.edu/SAVES/) analyses were also undertaken.
2.2. Molecular Dynamics Simulations
All MD simulations were performed by the software package Amber16 using the TIP3P for water molecules [ 13 ]. Before MD simulation, the systems were equilibrated in four steps: first of all, fixed the protein while allowed water to equilibrate for 2 ns. Secondly, the protein alpha carbons were then restrained by a harmonic potential with force constant reducing gradually from 10 kcal/mol to 0.1 kcal/mol in 4 steps and each step runs 10 ns. All bonds involving hydrogen atoms were constrained with the SHAKE20 algorithm. The formal simulations were carried out for a total of 260 ns at constant temperature (300 K) and pressure (1 bar) conditions maintained using a Langevin Piston and Langevin dynamics respectively. Six Cl-ions were added to obtain charge neutrality. To avoid artifacts, the MD simulations were run twice with different starting atomic velocities. The resulting trajectories were analyzed using the cpptraj module of the AMBER16 package.
2.3. Site-Direct Mutation
All chemicals were of analytical grade for biochemical use and were used as described by the manufacturer. Restriction endonucleases were purchased from New England Biolabs (MA, USA). All of the mutagenic oligonucleotides and PCR primers were synthesized by Shenergy (Shanghai, China). DNA sequences were determined by Bioasia (Shanghai, China). The detail procedure of protein assay and the thermal inactivation assays are all similar with the previous research [ 9 ]. One unit of chitosanase was defined as the amount of enzyme required to liberate 1 mol reducing sugar per min under the conditions described above. Three replicates were performed per analysis. The L84A mutations were separately introduced into the EAG1 gene using the QuickChange Multi Site-Directed Mutagenesis Kit (Stratagene). All mutations in the EAG1 were identified by DNA sequencing.
2.4. Protein Expression and Enzyme Assay
The recombinant plasmids were linearized with BglII and transformed into P. pastoris for protein expression. Positive transformants were initially grown in BMGY medium (25 ml) in 100 ml flasks at 28˚C with shaking (200 rpm) until the OD600 reached 4.5. The culture was centrifuged at 30,009 g for 10 min, the resuspended pellets were added to 1 l BMMY medium. To maintain expression, methanol was added every 24 h to give 0.5% (v/v). After 120 h induction, about 1l supernatant was concentrated to 80 ml by ultrafiltration (10 kDa cut-off). One ml concentrated crude enzyme was diluted in 5 ml of buffer A (20 mM sodium phosphate containing 500 mM NaCl and 20 mM imidazole, pH 8) and applied to Ni2+-charged 1 ml HisTrap FF crude column (GE Healthcare) equilibrated with buffer A. The recombinant EAG1 were eluted with buffer B (buffer A plus 400 mM imidazole). Chitosanase activity was determined at 40˚C by estimating the amount of the reducing ends of sugars using a modified dinitrosalicyclic acid (DNS) method with glucosamine HCl as the calibration standard. One unit of chitosanase was defined as the amount of enzyme required to liberate 1 lmol reducing sugar per min under the conditions described above. Three replicates were performed per analysis.
3. Results and discussion
3.1. Homology Modeling
The blast result showed that EAG1 shared 91% sequence identity with a chitosanase from bacillus circulans (PDBid: 1QGI), Hence, it was selected as the template. The automated sequence alignment and analysis of the template and target was made using the online service of [ 14 ]. The alignment (Figure 1) shows that EAG1 shares ~91% sequence identity with the chitosanase from bacillus circulans (PDBid: 2D05). The sequence differs mainly in the N-terminal. Therefore, the modeled structure of EAG1 is reliable.
The ROSETTA homology modeling method was used to model the structure of EAG1. Firstly, 10,000 models were constructed and then all of the models were selected through clustering and scoring. The modeled structure with highest Rosetta energy and dDFIRE score was selected. Then, the selected models was refined and checked by using Ramachandran plots (Figure 2(b)). It is found that the whole structure contains 12 α-helices and 5 β-strands (Figure 2(a)) that fold into upper and lower domains of comparable size. Residues 1 - 153 form the upper domain, it spans residues 1 - 22 and 46 - 176 and includes helices α1 - 7 and strands β1 - 5; the lower domain of 158 amino acids spans residues 23 - 45 and 177 - 311 and includes helices α7-12. The longest helix, α7, forms the backbone that connects the two globular domains. The Ramachandran plots for the selected model was shown in Figure 2(b) and the residue dihedral angles of all selected models are located in the allowed regions. Finally, the refined models was selected for further analysis.
3.2. Molecular Dynamical Simulation
The minimized structure was used as the initial structure for further MD simulation.300ns MD simulation were taken and the backbone RMSD of the conformation showed a rapidly increases to 3.5Å, then fluctuated around 3.5Å (Figure 3), indicating the model attained equilibration during MD simulation. This result demonstrated that during the whole dynamics simulation, the fluctuation of RMSD is mainly due to the flexibility of P-loops. The equilibrated structure is shown in Figure 2(a).
3.3. Prediction of the Key Residues for Stability
In this research, we identified the flexibility of EAG1 and its amino acid residues by using the molecular dynamics simulation and root mean square fluctuation (RMSF). The difference between RMSD and RMSF is that with the later the average is taken over time, giving a value for each particle. Highly flexible residues could trigger protein unfolding due to their large thermal fluctuations; thus, B-factors can be used as identifying characteristics to search for thermo labile residues [ 16 ]. By predicting the flexible motion of proteins, molecular dynamic simulation can provide more information about protein flexibility than the simple B-factors of static structures. RMSF calculates the mobility of atoms during the simulation. Thus higher RMSF values indicate higher mobility and lower RMSF values indicate restricted mobility. The amino acid residues have higher RMSF values may with higher flexibility. From Figure 4 we can easily find that residues Leu84, Gly113, Asp116, Ala207 and Leu286 show higher.
3.4. The Activity and Thermal Stability of L84A
The activity and thermalstability of the L84A were determined. The same as the wild-type EAG1, the optimal activity also showed at 30˚C. However, the activity of the L84A (8200 U/mg) was 1.4-folds higher than that of the WT enzyme (5860 U/mg). Compared with Leu, Ala has shorter side chains, which provides
Figure 1. Sequence alignment of EAG1 and 2D05. Secondary structure elements are indicated in blue. The relative accessibility of each residue is rendered as blue boxes from dark blue (accessible) to white (buried residues).
Figure 2. Representation and structure characters of EAG1 after refinement (a) the whole structure. This figure was generated by pymol [ 15 ], shown in rainbow cartoon. (b) Ramachandran plot.
Figure 3. Root-Mean-Square Deviation of the backbone atoms from EAG1. The horizontal axis is time (ns) and the longitudinal axis is RMSD (Å).
ample space around the catalytic triad and allows access of the substrate to the active site. So when Leu84 is changed into ALA84, at the same activation energy the substrate will become easily entry and exit the active site and the catalytic activity can be enhanced greatly.
Using the purified enzymes, the wild-type EAG1 and L84A were heated at 80˚C for 2 min to 10 min respectively. As heating time extended, EAG1 activity reduced gradually, and after being heated for 5 min, their residual activities were 11.42% and 34.57%, respectively (Figure 5). The subtraction mutant L84A showed 23.15% thermo stability enhancement compared with that of the wild-type. Moreover, the thermostability decreasement of L84Awas exhibited obviously just after being heated for 5 min. For the other four residues (Gly113, Asp116, Ala207 and Leu286).
Figure 4. The RMSF value of each residue of EAG1 after 200 ns MD simulatio207n, calculated by AmberTools2019. The horizontal axis is RMSF (Å) and the longitudinal axis is residue number of EAG1.
Figure 5. The thermo stability of wild-type EAG1 (Blue) and L84A (orange) after being heated at 80˚C for 2.5 min, 5 min, 7.5 min and 10 min. Enzyme activity without being heated was defined as 100%. All data presented here are the mean values derived from triplicate measurements.
Thermostability is a key factor in the applications of industrial enzymes. Generally, higher thermostability is often corresponding to a longer half-life. A 2˚C - 15˚C increase of thermostability often corresponds to a ten-fold longer life time. Here, we employed molecular dynamics simulation and root mean square fluctuation (RMSF) to identify the flexibility of EAG1 and its crucial amino acid residues. Highly flexible residues could trigger protein unfolding due to their large thermal fluctuations. RMSF calculates the mobility of an atom during the MD trajectory, thus higher RMSF values indicate higher mobility and lower RMSF values indicate restricted mobility. By molecular dynamical simulation, residues Gly113, Asp116, Ala207 and Leu286 have the higher RMSF values and they were thought to be the key residues for the instability of EAG1. Residues Asp116, Ala207 and Leu286 have been reported to the residues that can improvement in the thermostability of EAG1 by forming artificial disulfide bonds. By using the site-direct mutation, we demonstrated that L84A can also improve the thermostability of EAG1. These results will give information on the theoretical bases to improve the thermostability of EAG1.