The ß-mannanase gene was synthesized by Shanghai Sangon after codon optimization. pPICZ vector (preserved in our laboratory) was used to form a pPICZ-Anman recombinant vector with Anman. Escherichia coli JM109 and P. pastoris X33 were kept in our laboratory for cloning and heterologous expression of mannanase, respectively. Restriction endonuclease and DNA polymerase were purchased from Takara (Dalian, China). PCR purification kit was purchased from GENEWIZ (Suzhou, China). LBG was purchased from Sigma (St. Louis, MO, USA) as substrate. The rest of the chemical reagents not specified were analytically pure.

Because of the high protein sequence similarity between 3WH9 (PDB ID) and Anman (Huang et al., 2014), we selected 3WH9 as the template for Anman using the SWISS-MODEL server for homology modeling (Supplementary Figure 1). The ΔΔG of Anman calculated using the Discovery Studio 2019 software (Fu et al., 2017). Briefly, using the DeepView software to complement the missing side chains and residues, the Charmm36 force field was imparted to the protein, and ΔΔG was calculated under pH-dependent conditions. ΔΔG was calculated as shown in Equations (1) and (2):

ΔΔG_folding indicates the Gibbs unfolding free energy change of WT or mutant in the folding and unfolding states, and the mutant is predicted to be stable if it possesses negative ΔΔG_mut, and the opposite is predicted to be unstable.

The B-factor of wild-type mannanase Anman was calculated by ba2r,² and we subsequently obtained the secondary structure information of Anman using the PDBsum database.³ To reduce the loss of enzyme activity caused by the mutation effect, the average distance between catalytic residues and ß-turns in loops was measured using Pymol. Subsequently, we further narrowed down the range of mutation sites according to the amino acid preference of ß-turns and then combined it with ΔΔG screening to obtain hotspots. Meanwhile, Anman was compared with ß-mannanase with 13 excellent thermostability by multiple sequence alignment to explore semi-rational design’s role in finding hotspots.

The recombinant vectors of Anman and mutants were linearized by restriction endonuclease Sal al restriction endonuclease. The linear vector was transferred into P. pastoris X33 and incubated in yeast extract peptone dextrose (YPD) solid medium containing 100 mg/mL bleomycin (zeocin) at 30°C for 72 h. Colonies were picked and incubated in 30 Ml YPD liquid medium at 30°C and 220 rpm for 24 h, and cells were then collected by centrifugation at 8,000 rpm and transferred to 30 mL YP medium containing 0.5% methanol at 28°C and 220 rpm for 72 h. Methanol was added every 12 h to maintain its concentration at 0.5% for induction of enzyme production.

The induced fermentation broth was centrifuged at 4°C at 8,000 rpm, and the collected supernatant was dialyzed for 16 h in 0.1 M NaH₂PO₄-Na₂HPO₄ containing 20 mM NaCl (pH 5.5, buffer A) to remove salt and metal ions. Anion exchange chromatography was used for the purification of crude enzymes. Briefly, after pre-equilibrating the anion exchange column (Hi Trap™ 1 mL Q HP) with Buffer A, the dialyzed crude enzyme was loaded onto the column, followed by a linear gradient elution using Buffer A with Buffer B (Buffer A with 1 M NaCl), and the collection containing the enzyme activity was utilized for the subsequent study. The purification of the samples was analyzed by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) (Laemmli, 1970). The Bradford method was used to determine the protein concentration (Bradford et al., 1976).

The determination of ß-mannanase activity was performed by the 3,5-dinitrosalicylic acid (DNS) method. Briefly, the purified enzyme solution and LBG were diluted and dissolved with 0.1 M acetic acid-sodium acetate buffer (pH 3.0, Buffer C), and the reaction system consisted of 80 μL of pure enzyme and 160 μL of LBG (0.3 mg/mL), incubated at 37°C for 30 min and then terminated by adding 200 μL of DNS. After boiling for 5 min, 560 μL of water was added to the reaction system, and the absorbance of the reaction solution at 540 nm was measured using a microplate reader. Under the same conditions, the absorbance of mannan oligosaccharides was used as the standard curve. One unit of ß-mannanase activity was defined as the amount of enzyme that releases 1 μmol of reducing sugar per minute under the above conditions. All experiments were performed in three parallel experiments.

The specific activity of mannanase at different temperatures (50–80°C) was determined using the method mentioned above, and the temperature with the highest specific activity was determined as its optimum temperature. Thermostability was determined by incubating mannanase at 70°C for different times (10–60 min) and then reacting with LBG at 37°C to determine the residual activity. Tris-HCl (pH 2.0) and acetic acid-sodium acetate (pH 3.0–7.0) were used to determine the optimum pH of mannanase, the specific activity of mannanase was measured at different pH conditions (2.0–7.0), and the pH with the highest specific activity was determined as the optimum pH.

To determine the t_1/2 of mannanase, the enzyme was incubated at 70°C for different periods (5–90 min) and then reacted with LBG to determine its remaining specific activity. The kinetic parameters were determined using 1–10 mg/mL LBG as the substrate (concentration gradient of 1 mg/mL) with the pure enzyme at 37°C and pH 3.0 to determine specific activity. The K_m and K_cat parameters were obtained by non-linear fitting using the Origin 2017 software according to the specific activity and LBG concentration. All the above experiments were performed in three parallel experiments.

Molecular dynamics simulations of wild-type Anman and its mutants were performed using Gromacs 2018. Briefly, the missing atoms were added to the PDB file of Anman using the DeepView software before adding the Amberff14sb force field, followed by adding the TIP3P water model to the ortho-dodecagonal box and neutralizing the charge in the system with a sufficient amount of sodium ions. Next, the entire system was minimized using the fastest descent method (1,000 steps), followed by NVT and NPT simulations for 200 ns, respectively, using the LINCS algorithm to limit the heavy atoms and all bonds. The above processes were performed under periodic boundary conditions (PBC). Finally, the system was simulated at 370 K for 50 ns, with a trajectory saving interval of 2 fs. The parameters of root mean square deviation (RMSD), root mean square fluctuation (RMSF), and hydrogen bonds during the simulation were analyzed by the saved trajectory file, and VMD carried out the analysis of salt bridges.

The protein melting temperature (T_m) represents the thermodynamic stability and reflects the conformational changes of the enzyme at high temperatures. In this study, T_m was determined using the Nano-Differential SCNAning Calorimeter (Nano-DSC, Waters, USA) instrument. Briefly, the pure enzyme was diluted to a concentration of 0.2 mg/mL with Buffer C and loaded into a capillary cell. The scanning temperature range was 20–100°C with a ramp rate of 1°C/min, and Buffer C was used as the control. The data were analyzed using the Nano Analyze software.

Constraint network analysis (CNA) simulates the enzymes as body-and-bar networks, where the forces between atoms are replaced by bars. The number of bars represents the strength, and E_hb measures the magnitude of the energy contained in the force. CNA simulates a thermal unfolding process with increasing temperature, and the E_cut–hb decreases as the temperature increases. During the CNA simulation, there is a sudden decrease in the rigidity order parameter (P_∞), which means that the network changes from rigid to flexible, just like the protein folding-unfolding transition, and the E_cut–hb is the phase transition point. The temperature of conformational transition (T) can be calculated by using Equation (3):

Constraint network analysis was accomplished on a web service.⁴ In this study, the range of E_hb was set from –0.1 to –0.6 kcal/mol, and the decreasing gradient of E_cut–hb was –0.1 kcal/mol.

To identify the more flexible regions in Anman, the average B-factor of α-helices, β-sheets, and loops was calculated, which were 20.5, 17.7, and 21.4, respectively (Supplementary Figure 2). The regions with a higher B-factor level tended to contain more flexible residues; therefore, the present study focused on the mutation range in loop regions, and it has been also confirmed that the modification in loop regions could improve thermostability (Cheng et al., 2012). Due to the different flexibilities of the enzyme surface and internal loops, the B-factors of residues exposed on the surface vs. those buried in the interior were calculated as 24.1 and 18.2, respectively. Meanwhile, the calculated average B-factor of residues within 5 Å of the catalytic residue was only 14.8, suggesting that the protein interior was more stable than the surface. Thus, we narrowed the selection of mutation sites to loop regions located on the surface of Anman.

There are many β-turns on loop regions, thus, we can design mutations based on the amino acid preferences on different residues (position i, i + 1, i + 2, i + 3) of β-turns (Supplementary Figure 3). In addition, proline has a higher frequency at position i + 1, it is the most rigid of all amino acids, and Xaa→Pro is beneficial for the improvement of thermostability (Wang et al., 2014; Land et al., 2019). Thus, the mutation range was further narrowed to β-turns of Anman (types as furth, and we prefer mutation to proline. The Pymol software was used to calculate the average distance of C_a for each β-turn from catalytic residues Glu168 and Glu276, and their correlation with the B-factor is shown in Figure 1. Except for β-turn 6, all other β-turns had a higher correlation in B-factor and distance from catalytic residues, and β-turns with a higher B-factor tended to be further away from catalytic residues. To obtain mutants that do not affect the catalytic efficiency, several regions (β-turn 3, 7, 9, 11) that were closer to the catalytic residues were excluded, and the remaining β-turns were used for the next round of screening.

The remaining β-turns were subjected to virtual saturation mutation using the Discovery Studio 2019 software (Fu et al., 2017), and their ΔΔG was calculated. The heatmap in Figure 2 shows the ΔΔG for 380 mutants, and negative values represent mutations that are predicted to be stable. Since the calculation error for ΔΔG was about 0.5 kJ/mol, a total of 315 mutants were predicted to be neutral or deleterious (their ΔΔG > –0.5 kJ/mol). Some residues with high ΔΔG values (including I12, Y154, R197, V300, and L312) were excluded. Finally, six single-point mutations (E15C, S65P, G83T, A84P, A195P, and T298P) were selected using a β-turn preference for flexible modifications. D114A and E222T (predicted to be negative mutations) were selected for use as the controls due to their higher ΔΔG values (0.52 and 2.22 kJ/mol, respectively) and closer distances from catalytic residues (11.9 and 10.5 Å, respectively).

Using the wild-type Anman gene as the template, all mutants were generated using polymerase chain reaction point mutation and verified via DNA sequencing. The primers used for all mutants are shown in Supplementary Table 1. Anman and positive mutants were overexpressed using P. pastoris X33 and then purified and verified by SDS-PAGE. The crude enzyme in lane 2 had an inconspicuous protein band at 80 kDa, while all of the remaining lanes showed a distinct single band at 47 kDa, which was inconsistent with the calculated molecular weights (37.1 kDa) due to the processing of glycosylation modification of Anman before its extracellular secretion by yeast cells (Supplementary Figure 4).

Figure 3A shows the specific activity (black bars) and thermostability (red bars) of Anman and all mutants. Most of the mutants (excluding G83T and D114A) showed no significant decrease in a specific activity, which was due to the long distance between the mutation site and catalytic residue. As expected, the thermostability of mostly single-point mutants was improved, and A84P is the best-performing single-point mutant (from 46.2 to 68.7%). However, the relative activity of G83T was only 0.28-fold the level of Anman, which was an unexpected result. Finally, five single-site mutants with improved thermostability and no significant reduction in specific activity were combined to obtain the combinatorial mutant Mut5 (E15C/S65P/A84P/A195P/T298P). Mut5 had a 0.2-fold higher specific activity and a 0.7-fold higher thermostability than Anman (Table 1). It retained 82.4% of its activity after 10 min of incubation at 70°C, which meant that the obtained thermophilic β-mannanase mutant had a significantly improved thermostability without the loss of enzyme activity.

To further investigate the performance of the mutant and Anman at high temperatures, their residual activity was measured after incubation at 70°C for 0–60 min (Figure 3B). The single-point mutant showed a slight increase over Anman, and E15C became the single-point mutant with the highest residual activity after 60 min of incubation, still retaining 22.3% of its specific activity, while Anman was almost completely inactivated. The combined mutant Mut5 had a remaining activity of 55.73%, which was significantly higher than that of Anman throughout the incubation process, indicating that Mut5 was more suitable for long-term function and industrial production at high temperatures than Anman.

Anman and all mutants showed the highest activity at 72.5°C, with no significant change in the specific activity below 70°C. But, the activity of Anman and mutants decreased dramatically when the temperature was increased to 75°C (Figure 3C). The mutation effect did not increase the optimal temperature of β-mannanase, which implies that it is more challenging to engineer thermophilic enzymes than mesophilic enzymes. The optimum pH of both Anman and the mutant was 3.0, which is close to the acidic environment of the animal stomach, suggesting that Anman and Mut5 can be used as feed additives in the breeding industry (Figure 3D).

To further explore the reasons for the decrease in enzyme activity or thermostability of G83T, D114A, and E222T, the protein sequence of Anman was aligned with 13 GH5 family thermostable β-mannanases, and their derived strains are shown in Supplementary Table 2. The sequence conservation of E15, S65, A195, and T298 was less than 50%, G83 and E222 reached 92.3%, while D114 was up to 100% (Supplementary Figure 5). The comparison results indicated that G83, D114, and E222 were highly conserved residues, and their mutations occurred against the natural evolution, which was why they had negative mutation effects. The above results also demonstrated that relying on rational design alone could not completely exclude negative mutants, and introducing other methods in the screening process will be more beneficial in finding mutation sites. For example, Liu generated 44 mutants by modifying the surface-charged amino acids of ManAK, but only three mutants performed well (Liu et al., 2021).

To investigate the effect of mutation on the kinetics of Anman, the kinetic parameters of Anman and mutants were measured using 1–10 mg/ml LBG as the substrate (Table 1). Mostly, mutants had a similar K_m as Anman, but E15C reached 6.7 mg/mL. Mut5 showed a slight increase in K_cat, indicating that it had a higher affinity for LBG than Anman. There was no significant difference in catalytic efficiency K_cat/K_m for single-point mutants, but Mut5 was 0.25-fold higher than Anman, reaching 239 mL/mg s (Supplementary Figure 6). The above results revealed that the mutation did not impair the catalytic ability of Anman for LBG, which further verified that mutations farther away from the catalytic residues had a lower effect on enzyme activity and that it could be one of the conditions used for screening thermostable mutants.

The T_m for Anman was 75.7°C, which was close to its optimal temperature, and this was likely due to the protective effect of LBG on the catalytic pocket of β-mannanase. Except for A195P, the T_m for other mutants was higher than that for Anman, but the ΔT_m for all single-point mutants was less than 1°C, and the T_m value for Mut5 was only improved by 2°C compared to that for Anman (Supplementary Figure 7). Anman only had a 1.2% residual activity at 75°C, while Mut5 retained 36.0% activity (Supplementary Figure 8), and both Anman and Mut5 were completely inactivated when the temperature reached 80°C. This property made Mut5 more suitable for high-temperature conditions in industrial applications.

The t_1/2 was obtained by incubating at 70°C for different amounts of time. The t_1/2 for A84P was the highest of all single-point mutants and 1.2-fold higher than that for Anman. Although the t_1/2 for S65P and A84P was significantly higher compared to that for Anman, Mut5 was 7.8-fold higher than that for Anman, reaching 73.4 min, and this phenomenon was caused not only by epistatic effects but also by the modification of most of the flexible sites of Mut5. When there are multiple flexible sites on the surface, a single-point mutation often cannot prevent the process of thermal inactivation; therefore, if multiple or all flexible sites on the surface can be rigidified, enzyme thermostability will be dramatically enhanced, which is why Mut5 performed much better than the single-point mutants.

Constraint network analysis can be used to visually compare the overall rigidity of Anman and Mut5. The P_∞ of Anman and Mut5 decreased sharply when E_cut–hb was –1.84 and –2.76 kcal/mol, respectively (Figure 4A). The simulated temperatures were 336.8 and 355.2 K, according to Equation (3), which were consistent with the previous experimental results for T_m values. This also meant that Anman was inactivated earlier than Mut5. To further investigate the local and global stability of mannanase, a 50-ns kinetic simulation was carried out. It showed that the RMSD for Anman was usually higher than that for Mut5 (Supplementary Figure 9), and Mut5 had more hydrogen bonds for the first 45 ns (Figure 4B). These results indicated that the flexible engineered Mut5 was more stable than Anman in global conformation and more resistant to unfolding at high temperatures. MD simulation was also used to compare the local conformational stability of Anman and Mut5. The RMSF for all five mutation sites was reduced to various degrees, indicating that the mutation effect increased the rigidity of these sites (Figure 4C). The RMSF values in the region near the P65, P84, and P195 mutation sites were substantially reduced because the pyrrolidine proline ring did not only make the C_a and N atoms of the backbone more rigid in order to not rotate easily but also reduce the conformational flexibility of the nearby residues (Zhu et al., 2021). In addition, regions 24–35, 54–58, 106–117, and 307–321 also significantly reduced RMSF in the absence of mutations, probably because the mutation effect caused rearrangement of the internal force network through long-distance perturbation pathways (El Harrar et al., 2021). To investigate this phenomenon, the dynamic cross-correlation matrices (DCCMs) were calculated for Anman (Figure 5). There was a prominent cyan spot in the regions 54–58 and 106–117, which means that the atomic motions in these two regions were correlated and moved in the same direction, so their RMSF rise and fall trends should also be the same. In addition, region 54–58 was also spatially distributed adjacent to region 106–117, and there are many hydrogen bonds between the two regions; therefore, it is reasonable that the mutation causes a decrease in flexibility in both regions. The same process occurs in regions 24–35 and 307–321. The motion correlation of residues can be further used for the trade-off between thermostability and activity (Yu and Dalby, 2018).

The present study analyzed the chemical bonds within 5 Å of each mutation site (Figure 6). Several hydrophobic residues were located around E15, and the hydrophobicity of cysteine was much higher than that of glutamate. Thus, the E15C mutant site formed hydrophobic interactions with the surrounding residues to stabilize the region. A similar phenomenon was observed in lipolytic enzymes (Gallardo et al., 2010). The N atom of C15 also formed an additional hydrogen bond with the O atom of I12 (Figure 6A), which increased the thermostability of the E15C mutant. Even though A195P did not form more hydrogen bonds than Anman, D194 formed a salt bridge with R197 in A195P due to the structural perturbation caused by the mutation. To verify this phenomenon, the number of frames forming a salt bridge was divided by the total number of frames in the MD simulation. As a result, the probability of forming a D194-R197 salt bridge increased from 51.1 to 91.1%, and the distance between the OD2 atom of D194 and the NH2 atom of R197 decreased from 6.3 to 3.5 Å. Therefore, the improved thermostability of A195P was due to the formation of a salt bridge (Figure 6B). For the T298P mutant, no hydrogen bonds formed since the angle between the O atom of T298 and the N and H atoms of V300 was lower than 120°. In comparison, the mutated proline changed the local conformation and formed hydrogen bonds between P298 and V300, thereby reducing the flexibility of the local region (Figure 6C).

There were no hydrogen bonds near the S65 region of Anman, and there were four hydrogen bonds in S65P. The O atom of S64 formed two hydrogen bonds with the N and OG1 atoms of T67 (Figure 6D), which made the unfolding of the S65P region require more heat. Therefore, S65P did not easily deactivate at high temperatures. The same happens in the A84P region (Figures 6E, F), where the NE2 atom of Q88 approached the T82 residue to form an additional hydrogen bond, and the OD1 atom of D85 and the O atom of P84 formed two hydrogen bonds. A situation where one atom becomes an acceptor for two hydrogen bonds can strengthen the binding between residues. In addition, the S65P and A84P mutation sites were relatively close to each other in space, and the mutation effects could be transmitted between them (Dinmukhamed et al., 2020). Therefore, the force network between the S65P and A84P regions could be significantly changed or even redistributed, and the mutational effect resulting from structural rearrangement may also have a positive epistatic effect, which may account for the significantly reduced RMSF of residues near the two regions.