The pPIC9K plasmid and Pichia pastoris GS115 (Invitrogen, Carlsbad, CA, United States) were used for subcloning and overexpression. Escherichia coli DH5α (Weidi Biotechnology Co., Ltd., Shanghai, China) was used for gene cloning and DNA sequencing. Plasmid pPIC9K-pepsinogen was previously constructed and stored in our laboratory. Primers for mutation construction were synthesized by Ruimian Biotech Co., Ltd. (Shanghai, China). Plasmid Extraction and Gel/PCR Purification Kits were purchased from TOROIVD (BVI, United Kingdom). The Seamless Cloning Kit was purchased from Beyotime Biotechnology (Shanghai, China). All enzymes used in the experiments were supplied by Thermo Scientific (Waltham, MA, United States). The bovine hemoglobin substrate was purchased from Maokang Biotech Co., Ltd. (Shanghai, China). All other chemicals were of analytical grade and are commercially available.

Site-directed mutagenesis was performed using the Seamless Cloning Kit according to the manufacturer’s instructions, with the plasmid pPIC9K-pepsinogen (with a C-terminal His-tag) as the template. The primers used for mutagenesis are listed in Supplementary Table S1. Subsequently, the seamless cloning products were transformed into E. coli DH5α cells. Recombinant plasmid construction was confirmed by DNA sequencing. They were then linearized with SalI, and transformed into P. pastoris GS115 competent cells through electroporation. Positive clones were screened using minimal dextrose (MD) plates, their activity was demonstrated using skim milk plates, and confirmed by DNA sequencing.

Positive clones were fermented as previously described [24]. Cultures were centrifuged at 8000 ×g for 10 min at 4°C and crudely purified by the slow addition of ammonium sulfate until 80% saturation was achieved. After overnight dialysis against TE buffer (20 mM Tris-HCl, 500 mM NaCl, pH 7.5) at 4°C, the crude samples were purified using a HisTrap HP affinity column (General Electric Company, Boston, MA, United States). The fractions containing pepsinogen were pooled and dialyzed against the TE buffer, then activated by adding 0.9 M HCl to adjust the pH to ∼2 and subsequently incubated at 25°C for 10 min. The prosegment was removed using Centricon filtration units with a 10 kDa molecular mass cutoff (Millipore, Boston, MA, United States). All samples were stored in 20 mM sodium acetate (pH 5.3) and analyzed using 12% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE, Coomassie Brilliant Blue R250-stained). Finally, the concentration of recombinant pepsin was measured using a BCA protein quantificationkit (Yeasen Biotech Co., Ltd., Shanghai, China).

The proteolytic activities of pepsin and its variants were determined by measuring the amount of exposed tyrosine residues during hydrolysis in 100 mM sodium citrate buffer (pH 2.0), with 1% bovine hemoglobin as the substrate. A mixture containing 0.5 ml of 1% bovine hemoglobin and 0.5 ml of appropriately diluted recombinant proteins was incubated at 60°C for 10 min. The reactions were terminated by adding 1.0 ml of 0.4 M trichloroacetic acid (TCA). Exposed tyrosine residues were quantified by incubating reactions containing 0.5 ml of the filtered supernatant, 2.5 ml of 0.4 M sodium carbonate solution, and 0.5 ml of Folin’s phenol reagent at 40°C for 20 min and then measuring absorbance at 680 nm. Control samples were treated in a similar manner, but TCA was added before adding the enzyme.

The initial reaction rates of recombinant proteins were also determined at 60°C and a pH of 2.0, using 0.5–10.0 mg ml⁻¹ bovine hemoglobin as the substrate. Kinetic parameters (K _m, k _cat , and V _max ) were obtained by fitting the initial data to the Michaelis–Menten nonlinear regression equation in the Origin 9.0 software (Origin Lab Corporation).

To evaluate thermostability, the half-lives (t _1/2 ) of pepsin and its mutants were measured at 65°C and 70°C. Recombinant proteins were appropriately diluted with 100 mM sodium citrate buffer (pH 2.0) and pre-incubated at 65°C and 70°C. Samples were taken at regular intervals, and the residual activities measured at 60°C and pH 2.0, as described in Enzymatic Activity and Kinetic Assays. Initial activity, before pre-incubation, was set to 100%.

The crystal structure of pepsin was obtained from the Protein Data Bank (PDB: 4PEP), and all the crystal waters were removed. The structures of the mutants were prepared using Discovery Studio 2019. The Leap module of AMBER18 was used to fill the missing atoms. The parameters were generated using the AMBER14SB force field. The protein was placed in a TIP3P water box with 12 Å buffer, and counter ions were added to neutralize the system. The steepest descent and conjugate gradient methods were used to optimize the structure. The entire system was subsequently heated to 310 K and the protein was constrained to 10 kcal/(mol Å²). The SHAKE algorithm was employed to constrain chemical bonds involving hydrogen atoms [25]. Finally, simulations were run for 100 ns.

All statistical analyses were performed using Origin 9.0 software, with values expressed as the mean ± standard error of the mean. A one-sample Student’s t-test was used to compare the mutant and wild-type means. Statistical significance was set at p < 0.05. For comparisons between the wild-type and mutants, ^*represents p < 0.05, **represents p < 0.01, and ***represents p < 0.001.

Members of the aspartic protease family have a similar catalytic mechanism. A nucleophile is formed from a water general-base activated by Asp215, which subsequently attacks the carbonyl oxygen of the scissile bond to form a gem diol of a tetrahedral intermediate with Asp32. The proton transfer from Asp215 to the leaving group is coupled with the regeneration of the proton of Asp32 [26]. Finally, the peptide bond is broken. Thus, the catalytic dyad composed of Asp32 and Asp215 is crucial for the catalysis of pepsin. Two key points need to be avoided in the thermostability design.

Many aspartic proteases have been reported to have similar three-dimensional structures. For simplicity, the pepsin structure can be divided into three regions: N-lobe, C-lobe, and connecting domain (Figure 1). The N-lobe is made up of a series of ß-strands that are folded into three layers of orthogonal packing: sheets IV, VI, andII.The C-lobe has two orthogonal packing sheets, VII and III, and a rather separate sheet V [27]. Several a-helices and loops are also present in both lobes. The connecting domain is comprised mainly of a six-stranded antiparallel interdomain ß-sheet connecting the N-lobe and the C-lobe to form the catalytic dyad.

Structural analysis revealed that the three ß-sheets in the N-lobe form a continuous hydrophobic region, and the a-helix near the connecting domain and adjacent loops cap the hydrophobic core, further stabilizing the structure. Similarly, a hydrophobic core is formed in the C-lobe, and the connecting domain simultaneously interacts internally with the N- and C- lobes [27]. The loops in pepsin are mainly ß-turns or ß-hairpins that connect ß-strands. We speculated that pepsin may possess two thermodynamic structural weaknesses: the relatively loose internal hydrophobic packing and the high flexibility of the surface a-helices and loops that are not part of the hydrophobic cores. Therefore, there are two strategies for improving the thermostability of pepsin: 1) designing ß-strands that directly improve the interaction between ß-sheets or ß-strands, and 2) designing a-helices or loops on the surface to reduce local disorder.

To investigate correlations between computational design software programs and structural characteristics, we selected several thermostability prediction software programs with successful cases and with varying principles. PoPMuSiC is a web server that predicts folding free energy changes ( Δ Δ G ) caused by single-site mutations in proteins. Δ Δ G is a linear combination of statistical potentials, considering the amino acid type, torsion angles defining the backbone conformation, and solvent accessibility. Δ Δ G < 0 indicates that the mutant is more stable than the wild-type. The thermostability of multiple proteins has been improved using the PoPMuSiC algorithm [9, 28, 29]. HoTMuSiC is a computational tool that uses an artificial neural network to predict the change in melting temperature Δ T m , upon the introduction of single-point mutations, based on the imperfect correlation of thermodynamics and thermal stabilities [30]. Δ T m > 0 indicates that the mutants are more stable than the wild-type. Enzyme thermal stability system (ETSS) is an enzyme redesign protocol that improves stability based on the TK-SA model calculation and surface charge–charge interaction analysis protocol [31, 32]. Charged amino acids are computed using the ETSS and residues with Δ G q q > 0 are deemed unstable for the given protein and considered to be uncharged or oppositely charged amino acid mutations. ETSS is effective at improving the thermostability of many enzymes [33, 34]. DeepDDG is a neural network-based method used to predict changes in protein stability resulting from single-point mutations. Δ Δ G > 0 indicates that the mutant is more stable than the wild-type. The neural network was trained on the largest dataset of experimentally determined stability data from various sources, and it predicted stability changes based on the three-dimensional structure of a protein. This algorithm gave Pearson correlation coefficients of 0.48–0.56 for three independent test sets, outperforming 11 other methods [35]. Fireprot is a web server developed by Musil et al., which integrates 16 computing tools, including Clustal ω, Rosetta, and FoldX [36]. The algorithm combines energy- and evolution-based approaches to predict stable mutations. Prior to virtual mutation, correlation and conserved analyses are performed to exclude important residues. Surprisingly, the phylogenetic analysis enabled the identification of mutations stabilized by entropy, which could not be predicted by force field calculations [36].

We chose the above-mentioned representative and unique software programs. By analyzing the successful cases of these software programs, mutations were found in a-helices, ß-strands, and loops distributed in the interior and surface of proteins, thus meeting the requirements of our analysis of thermodynamic weaknesses. Among these software programs, PoPMuSiC, HoTMuSiC, ETSS, and DeepDDG provided mutation rankings. Based on the highest-ranking batch and the principle of scattered mutation points, we selected four mutations predicted by each software program for follow-up experimental verification. The top 20 mutations predicted by PoPMuSiC, HoTMuSiC, ETSS, and DeepDDG are shown in Supplementary Tables S2–S5 to assist with explaining the screening principles. The number of mutations given by Fireprot was limited; thus all mutations were selected for further experiments.

Based on the screening described above, 34 single-site mutants were constructed and transformed into P. pastoris for expression, using the pPIC9K-pepsinogen plasmid as the template. The mutants predicted based on energy (from PoPMuSiC, HoTMuSiC, ETSS, DeepDDG, and Fireprot-Energy) were tested as a set. To rapidly screen the mutants, we used crude enzyme to measure thermostability. One unit of enzyme activity was defined as the amount of enzyme required to hydrolyze bovine hemoglobin toproduce 1 μg tyrosine per minute at 60°C and pH 2.0. HCl (0.9 M) was added to the fermentation broth supernatants and the pH was adjusted to 2.0, to convert pepsinogen into pepsin. The acidified fermentation broths were incubated at 65°C for 3 min and then cooled on ice. Changes in enzyme activity before and after incubation were measured to reflect the thermostability of recombinant proteins. Enzyme activity before incubation was plotted at 100%. The results are shown in Figure 2A. Apart from the mutations selected by energy-based Fireprot, some mutants predicted by PoPMuSiC, HoTMuSiC, ETSS, and DeepDDG resulted in complete loss of enzyme activity before incubation. These occurred mainly at Gly85, Asp87, Gly102, and Arg315. The remaining mutants did not improve pepsin thermostability. These results indicate that the energy-based mutants were not associated with the thermodynamic weakness of pepsin, and some were located on critical pepsin residues, resulting in the collapse of the hydrophobic cores.

To guide the follow-up design, we analyzed the positions of the above mutations to confirm the structural weakness of pepsin. As shown in Figure 2B, Ala66 and Asp159 are on loops, Tyr114 and Gln232 are on a-helices, and the remaining mutation sites are on ß-strands. Thr28, Gly85, Asp87, Asp96, and Gly102 are located on the three ß-sheets of the N-lobe, which are stacked orthogonally to form hydrophobic cores; Asp3, Glu4, Arg315, and Lys319 are located on the closely packed ß-sheet in the connecting domain; Thr222, Thr261, Thr283, and Tyr309 are located in the orthogonally packed ß-sheets in the C-lobe. The a-helices where Tyr114 and Gln232 are located are also close to hydrophobic cores; the Gln232 side chain is oriented towards the protein, and Tyr114 is an aromatic amino acid. Thus, most of the mutations are located on or near hydrophobic cores, and changes in hydrophobic cores can easily cause enzyme inactivation [37]. The results showed that pepsin’s hydrophobic cores were tightly packed, and were not the major thermodynamic structural weakness of the enzyme. Mutations in the compact region can have long-range effects. Therefore, when further strengthening the internal packing, mutations in the ß-strands and a-helices involving hydrophobic cores lead to enzyme inactivation. Although the computational design results showed that the overall energy of these mutants was greatly reduced, the probability of protein structure destruction was very high. Thus, mutations on ß-strands or inside a-helices are not suitable and mutations guided by the real structural weaknesses—flexible regions on the surface of pepsin—need to be selected.

Fireprot also identified 10 mutations based on evolution (Table 1). As mentioned above, mutations on tightly packed sheets and a-helices are less feasible, and surface a-helices and loops are the main thermodynamic weakness points of pepsin. Thus, we initially analyzed the positions of these 10 mutations. As shown in Figure 3A, E65Q is on the orthogonally packed ß-sheet of the N-lobe, L167F is on the tightly packed ß-sheet of the connecting domain, D200N is on the ß-turn of the C-lobe, L140M and W141M are on the a-helix, which is part of the hydrophobic core, and the other five mutations are on the loops or a-helices on the surface of pepsin. Thus, A24P, D52N, Q55R, S129A, and S270P mutants were set as the experimental group, and E65Q, L140M, W141M, L167F, and D200N were set as the control group. As described above, we measured the thermostability of the crude enzymes, the results of which are shown in Figure 3B. The wild-type pepsin retained 32.4% of its original activity after incubation for 3 min at 65°C. Among the mutants in the experimental group, the thermostabilities of A24P, D52N, E65Q, and S129A were significantly improved, with D52N and S129A retaining 75.4 and 63.5% of their initial activity, respectively. Only S270P showed a slight decrease in thermostability. The thermostability of mutants in the control group decreased, particularly L140M and W141M. The results indicate that the above strategy of being guided by the thermodynamic weakness of pepsin can rapidly and effortlessly guide mutant design.

Mutants A24P, D52N, E65Q, and S129A, which had improved thermostability, were purified (Supplementary Figure S1). The thick band around 35 kDa is unglycosylated pepsin, and the fuzzy fine band above it is O-glycosylated pepsin, both of which are consistent with previous reports. Glycosylation is used to enhance the thermostability of proteins, and Yoshimasu et al. confirmed that the thermostability of O-glycosylated pepsin was improved [24]. Therefore, it was retained in the experiment. The half-lives (t _1/2 ) at 65°C and 70°C were measured and 1% bovine hemoglobin was used as the substrate to determine changes in thermostability (Table 2). The t _1/2 values of A24P and Q55R at 65°C and 70°C were slightly improved. The t _1/2 of S129A at 65°C and 70°C increased by 55.6 and 66.3%, respectively, while that of D52N at 65°C and 70°C increased by 308.1 and 200.0%, respectively.

The activities of the mutants D52N and S129A were determined at pH 2.0 and 60°C, using 1% bovine hemoglobin as the substrate. The kinetic parameters were also measured under these conditions using gradient concentrations of bovine hemoglobin (Supplementary Figure S3). As shown in Table 3, compared with the wild-type, the activities of the D52N and S129A mutants decreased by 31.8 and 23.3%, respectively. The K _m values of D52N and S129A were 2.12 ± 0.05 and 1.40 ± 0.09 mg ml⁻¹, 130.4 and 52.2% lower than that of the wild-type, respectively. Thus, although k _cat values improved, the k _cat and k _cat /K _m values of D52N and S129A mutants decreased by 37.0 and 27.0%, respectively.

Dynamic simulation analysis was performed to explore the mechanism of the effect of mutations on the thermostability of pepsin and its mutants. Single-point mutations D52N and S129A were introduced into the crystal structure of pepsin using Discovery Studio 2019, and a stable structure was obtained through a 100 ns MD simulation. The root mean squared deviations (RMSDs) of the protein backbones are shown in Supplementary Figure S2. The RMSD of the wild-type and the D52N and S129A substituted structures were approximately 2.5 Å, 1.5 Å, and 1.5 Å, respectively, indicating that the three systemsare stable, and that D52N and S129A are more stable than the wild-type. A snapshot was extracted from each of the three stable structures for structural analysis (Figure 4). The D52N substitution was located on the a-helix on the protein surface, and formed hydrogen bonds with Leu48 and Ala49, making the a-helix more stable and improving structural rigidity. The S129A substitution was located in a loop on the protein surface. The replacement of Ser with Ala promoted the hydrophobic packing of Pro126, Gly132, Ala133, and Pro135, making this region more stable.

Furthermore, isotropic temperature factors (B-factors) of the three systems in the last 10 ns of the simulation were calculated to measure the residue fluctuations relative to their mean positions. As shown in Figure 4, we found that the D52N substitution reduced the flexibility of position 52, and the flexibility of Leu48 and Ala49 also decreased, but to varying degrees. After S129A substitution, the flexibility of position 129 decreased, and the flexibility of Pro126, Gly132, and Pro135 of the residues forming the hydrophobic packing decreased by different degrees, with the biggest decrease being noted in Pro126 and Pro135.