How Modification of Accessible Lysines to Phenylalanine Modulates the Structural and Functional Properties of Horseradish Peroxidase: A Simulation Study

Horseradish Peroxidase (HRP) is one of the most studied peroxidases and a great number of chemical modifications and genetic manipulations have been carried out on its surface accessible residues to improve its stability and catalytic efficiency necessary for biotechnological applications. Most of the stabilized derivatives of HRP reported to date have involved chemical or genetic modifications of three surface-exposed lysines (K174, K232 and K241). In this computational study, we altered these lysines to phenylalanine residues to model those chemical modifications or genetic manipulations in which these positively charged lysines are converted to aromatic hydrophobic residues. Simulation results implied that upon these substitutions, the protein structure becomes less flexible. Stability gains are likely to be achieved due to the increased number of stable hydrogen bonds, improved heme-protein interactions and more integrated proximal Ca2+ binding pocket. We also found a new persistent hydrogen bond between the protein moiety (F174) and the heme prosthetic group as well as two stitching hydrogen bonds between the connecting loops GH and F′F″ in mutated HRP. However, detailed analysis of functionally related structural properties and dynamical features suggests reduced reactivity of the enzyme toward its substrates. Molecular dynamics simulations showed that substitutions narrow the bottle neck entry of peroxide substrate access channel and reduce the surface accessibility of the distal histidine (H42) and heme prosthetic group to the peroxide and aromatic substrates, respectively. Results also demonstrated that the area and volume of the aromatic-substrate binding pocket are significantly decreased upon modifications. Moreover, the hydrophobic patch functioning as a binding site or trap for reducing aromatic substrates is shrunk in mutated enzyme. Together, the results of this simulation study could provide possible structural clues to explain those experimental observations in which the protein stability achieved concurrent with a decrease in enzyme activity, upon manipulation of charge/hydrophobicity balance at the protein surface.


Introduction
Horseradish peroxidase (HRP) is a Classical Secretory plant peroxidase isolated from horseradish (Armoracia rusticana) roots which catalyzes the oxidation of a wide variety of substrates, using hydrogen peroxide or organic peroxides [1,2]. HRP C, the most studied heme peroxidase, consists of a single polypeptide chain of 308 amino acid residues, including a heme prosthetic group, two calcium ions, four disulfide bridges, and eight carbohydrate chains [3,4]. The first crystal structure of glycan-free recombinant HRP has been reported in 1997 [5]. Since then, the three-dimensional (3D) structure of catalytic intermediates and several substrate complexes of HRP have also been reported and subjected to molecular dynamics (MD) simulations [6][7][8][9][10][11][12]. Some major advances in understanding the structural and functional aspects of HRP C have been achieved through these simulations, including the role of key residues in substrate binding, reactivity and stability [10][11][12].
HRP contains six lysine residues, among them K174, K232, and K241 are located on the protein surface [3]. It has been reported that thermostability of HRP increases by chemical modification of these three accessible lysines, but the stabilizing effect is lost after complete modification of all six lysine residues [18].
Most of the stabilized derivatives of HRP reported to date have involved chemical or genetic manipulations that neutralize or reverse the positive charges on lysine residues. Accordingly, these modifications can be classified into three distinct groups: The first and second are those which alter lysines to neutralized aromatic or aliphatic hydrophobic residues. The third group, are chargereversing modifications. In this study, we modeled a case from the first group by mutating the most accessible lysine residues to phenylalanine. In genetic manipulations, it is possible to target a specific residue, but there is no exact control over the chemical modifications and the accessibility of residues determines the rate and location of modifications. Moreover, it has been shown that increasing the number of modified lysine residues up to three increases the stability gain in HRP [18]. To mimic the case globally, all the three accessible lysine residues were modified simultaneously.
In earlier studies, we reported covalent attachment of an electron relay (anthraquinone 2-carboxylic acid) to the surfaceexposed lysine residues of HRP. Experimental results showed that the modification alters all three accessible charged lysines and improves electron transfer properties, catalytic efficiency, and stability of the enzyme [29,30]. MD simulations clarified some structural changes relating to stability and activity enhancement, including decreased flexibility of the protein backbone and redistribution of hydrophobic patches at the protein surface as well as peroxide and aromatic binding sites [30]. Here, we try to study a different case in which accessible charged lysines (K174, K232, and K241) are substituted by aromatic hydrophobic phenylalanine residues. We use MD simulation to investigate the structural and dynamical consequences of these substitutions.

Materials and Methods
All Molecular dynamics simulations were performed using the GROMACS simulation package version 4.0.7 [40] and GRO-MOS96 force field [41]. The starting atomic coordinate of native HRP (n-HRP) was obtained from Protein Data Bank (PDB: 1ATJ) [5]. To generate the initial structure of the modified protein, p-HRP, the lysine residues 174, 232, and 241 of native HRP were modified to Phenylalanine residues. Each protein was centered in a cubic box and immersed in SPC water molecules so that the shortest distance between the protein and the box boundaries was 1.3 nm and periodic boundary conditions were applied. To achieve a neutral simulation box, the net charge of the protein was neutralized by replacing water molecules with necessary Cl 2 or Na + ions. Each solvated and neutralized system was energyminimized using the steepest descent algorithm until the maximum force was smaller than 500 kJ/mol.nm. After energy minimization, two separate position-restrained MD simulations were sequentially carried out to adjust temperature and velocities and to equilibrate the solvent and ions around the protein. First, to adjust the system temperature, an NVT MD simulation was performed for 200 ps at 300 K by imposing thermal energy in a constant volume condition using the velocity rescale algorithm (modified Berendsen thermostat) with t T = 0.1 ps [42,43]. After arrival at the correct temperature, the resulting atom velocities and coordinates were used to start an NPT MD simulation at 300 K and 1 bar for 200 ps by the Parrinello-Rahman algorithm with t P = 2.0 ps during which density of the system was stabilized at around 1000 kg/m 3 [44]. Finally, the production MD period of 100,000 ps at constant pressure (1 bar) and temperature (300 K) without position restraints was performed on native and p-HRP. Bond lengths were constrained using LINCS algorithm [45]. In Gromos96 force field, the Lennard-Jones and short-range electrostatic interactions are optimized with a cut-off radius of 1.4 nm. However, in order to reduce the computation time, we used a shorter cut-off of 1.0 nm based on the tested insensitivity of HRP simulation down to this value. The particle mesh Ewald algorithm was used for the long range electrostatic interactions [46]. The neighbor list was updated every 5 steps. Each component of the system was coupled separately to a thermal bath, and isotropic pressure coupling was used to keep the pressure at the desired value. A time step of 2 fs was used for the integration of equation of motion.
In order to evaluate the quality or sufficiency of conformational sampling of simulations, the production MD period was divided into four 25 ns parts and principal component analysis was performed on each sub-trajectory. Eigenvectors and eigenvalues were obtained from the diagonalization of the covariance matrices of the C a atoms, and the principal components were generated by projecting the trajectories on the respective eigenvectors [47]. The cosine content of the principal components was calculated to estimate whether the conformational fluctuations are connected with the potential (when the cosine content is close to 0) or with random diffusion (when the cosine content is close to 1) [48].

Global Dynamics
Using the known x-ray crystallographic structure of HRP C (PDB: 1ATJ), two 3D molecular models of HRP C, n-HRP and p-HRP, were constructed differing in the residues 174, 232, and 241 ( Fig. 1). Keeping the amino acid pattern of the original X-ray structure, the control model (n-HRP) was built with cationic lysines at these positions. In the next model (p-HRP) charged lysine residues were substituted by hydrophobic aromatic phenylalanine residues. The overall stability and structural relaxation of the enzymes were monitored by computing time evolution of the root mean square deviation (RMSD) of the backbone atoms along the simulations. The RMSD of n-HRP and p-HRP with respect to their starting structure were 2.4660.47 and 3.0360.52 Å , respectively. This indicates that the structure of mutated protein undergoes more conformational changes along the simulation. Comparison of the RMSD per residue profiles of the MD trajectories revealed that these changes are concentrated in the N terminal part of the protein structure (including residues 1-12), C terminal part of the loop D9E (residues 140-144) and the helix E (residues 145-153) (data not shown). For n-HRP model, the backbone RMSD as a function of time reaches a relative plateau after about 50 nanoseconds (ns) of simulation, but for p-HRP such a plateau is achieved after about 70 ns ( Fig. 2A). Hence, to be statistically comparable, our analyses were focused on those trajectories obtained from the last 30 ns of simulations for both the native and mutated proteins. In order to check whether the sampling of conformational space is sufficient, the production MD period was divided into four 25 ns parts and the cosine content of the first 4 principal components was calculated for each subtrajectory. For both n-HRP and p-HRP models, the best results of cosine content were obtained for the last quarter part of the 100 ns simulations. The values of the first 4 principal components for this time window were calculated to be 0.146, 0.459, 0.036, 0.110 for n-HRP and 0.197, 0.517, 0.010, 0.001 for p-HRP.The stability of the fluctuation of the potential energy was also examined by calculating the ratio between the variance and average of potential energy. This ratio for n-HRP and p-HRP was respectively about 0.00099 and 0.00092, thus showing that energy was conserved during the simulations and providing additional evidences indicating that the simulations and models were stabilized.
In the last 30 ns of simulations, the backbone RMSD of n-HRP and p-HRP structures relative to their own 70th ns structures (i.e. the starting structures of the analyzed time frames) was calculated to be 1.6460.33 and 1.0860.21 Å , respectively (Fig. 2B). In agreement, the root mean square fluctuation (RMSF) values for this time window also show a clear decrease upon substitutions (Table 1), implying reduced flexibility of the protein backbone as a result of these mutations. To provide a more detailed description of the mobility of protein residues, the backbone RMSF per residue for the native and mutated HRP is shown in Fig. 3. According to the crystallographic structure (PDB: 1ATJ), HRP contains thirteen a-helices and two short antiparallel b-strands. It is obvious that the mobility of regular secondary structures, except for the helix F9 (including residues 181-183), is mainly reduced after substitutions. Reduced mobility is particularly significant in helices F0 and G, as well as the strands b1 and b2. Apart from regular structures, mobility of the connecting loop b2G is also greatly diminished. Furthermore, the mobility of protein backbone around the heme prosthetic group is declined from 1.10 Å in n-HRP to 0.69 Å in p-HRP. Taken together, it can be concluded that the overall backbone structure has become less flexible upon these mutations. Therefore, decreased mobility of the heme prosthetic group is to be expected (Table 1).
Additional information on the structural flexibility is offered by the analysis of time-dependent secondary structure fluctuations. The secondary structure content of the models was calculated as a function of time using the DSSP program [49] (Fig. 4 and Table 2). Analysis of Fig. 4 reveals that in both models, a-helices exist persistently throughout the simulations, although the values in Table 2 indicate a small increase in the a-helical content of p-HRP which is mainly due to the stability of two short helices F9 (residues 181-184) and G (residues 231-239). By comparing the occurrence of these helices along the simulations (Fig. 4), it becomes clear that both of them, particularly helix G, are more persistent in p-HRP. It is noteworthy to remind that the substituted residue 232 resides in this helix. In addition, the helix  patterns assigned by DSSP show subtle differences between the two models. The end of the helix A is extended by three residues (Arg-27, Ser-28, and Asp-29) in p-HRP compared with n-HRP. A similar trend was also observed in the helix D during the first 5 ns of analyzed time frames.
In order to gain a better insight into the effects of these substitutions on HRP structure, local changes around the substitution sites were also investigated. As mentioned earlier, HRP contains a b-sheet comprised of two short antiparallel bstrands, b1 and b2, which flank the helices F9 and G, respectively. K174 is located in strand b1, at the proximal side of the heme cavity. DSSP analysis of the MD trajectories shows that while in p-HRP, the b-sheet structure is maintained almost throughout the entire analyzed time frames, it is absent in 36.25% of the n-HRP trajectories and replaced by b-bridge approximately in one third of the simulation time. This is in accordance with the observed increase in existence percentage of those hydrogen bonds connecting backbone of the residues involved in b-sheet formation in p-HRP model (Table 3A). It seems that improved hydrogen bonding in this region accounts for lower mobility of the strands b1 and b2 in modified HRP compared with the native one. In addition to the changes mentioned above, K174F improves hydrogen bonding between protein moiety and the heme prosthetic group. Backbone amide of Phe-174 forms a strong hydrogen bond (occupancy 98.93%) with heme O2D propionate oxygen, whereas in n-HRP this hydrogen bond is almost absent. Furthermore, this substitution strengthens the hydrogen bonds of Ser-35 and Arg-31 with propionate oxygens of the heme prosthetic group. Together, these findings suggest that Phe at position 174 may have a local stabilizing effect on p-HRP structure.
Residues 232 and 241 are located in the helix G and the connecting loop GH, respectively. Helix G is longer by three residues in p-HRP than in n-HRP (Fig. 4). On the other hand, the data presented in Table 3B suggests that this helix is more stable in p-HRP due to the formation of more and stronger hydrogen bonds. Replacing Lys-232 with Phe also improves the hydrogen bonding strength between Phe-232 and Asp-230, Thr-171 and Asp-230, Asp-230 and Ile-228, Thr-225 and Asp-222 which are involved in proximal Ca 2+ coordination (Table 4). This substitution also introduces a new hydrogen bond between Asn-236 and Asp-222 that are bound to Phe-232 and proximal calcium ion, respectively. Such changes would enhance the stability of proximal Ca 2+ binding pocket. Moreover, the orientation of the connecting loop b2G between helix G and b-sheet has changed after modification. In p-HRP, formation of a new hydrogen bond between the side chain nitrogen of Asn-236 in helix G and the backbone oxygen of Asp-222 in loop b2G pulls the loop b2G towards the helix G. Also in p-HRP, Gln-240 and Gly-242 (from the connecting loop GH) which flank the substitution site 241 form two persistent hydrogen bonds with Asn-198 and Thr-196 (from the connecting loop F9F0), respectively. These new hydrogen bonds stitch end of loop F9F0 to the beginning of loop GH and enhance the integrity of this part of p-HRP structure.
Comparison of the structures also reveals that the connecting loop b1F9, helix F9, connecting loop F9F0 and helix F0 (residues 177-208) exhibit considerable conformational changes (Fig. 1B). The relative orientation of this part of p-HRP model, especially  To assess the effect of mutation of three positively charged lysine residues on the exposure pattern of surface charge/hydrophobic-ity, the water accessibility of hydrophobic and hydrophilic surface areas in each model was measured and averaged over the analyzed time frames. Surprisingly, both the hydrophobic and hydrophilic exposed areas are harmoniously declined in p-HRP (Table 1), so as the ratios between the areas remain almost identical after substitution (0.9696 in n-HRP and 0.9731 in p-HRP). Consistent with the total surface area, Table 1 shows that the average values for radius of gyration and total protein volume are also decreased upon substitution. Comparing the average distance between two structural calcium ions shows a decrease from 31.4262.03 Å in n-HRP to 28.2460.30 Å in p-HRP. This finding is particularly important because the substrate binding cleft in HRP is located  between the two calcium ions and any closing of this cleft could profoundly affect the heme availability and consequently the enzyme activity. Fig. 5 shows the existence map of hydrogen bonds for n-HRP and p-HRP models over the corresponding time frames of simulations. Analysis of the internal protein hydrogen bonds evinces that the number of stable hydrogen bonds is increased in p-HRP. For example, the total number of intramolecular hydrogen bonds with existence percentage greater than 80% has increased from 103 in n-HRP to 115 in p-HRP (Table 5). In agreement, the percentage occupancy of hydrogen bonds forming b-sheet and helix G in p-HRP is significantly higher than that of n-HRP (Table 3). Together, these findings suggest that consolidation of intramolecular hydrogen bonding network actively contributes to the stability enhancement of p-HRP.

Peroxide-binding site
HRP utilizes hydrogen peroxide (or organic peroxides) to oxidize a broad range of organic and inorganic substrates. The main overall reaction catalyzed by HRP can be summarized as H 2 O 2 + 2AH R 2H 2 O + 2A*, where AH represents a reducing substrate and A* is a free radical product. Hydrogen peroxide oxidizes the native enzyme to an intermediate that called Compound I, which then can oxidize the reducing substrate [2]. To react with the heme prosthetic group, hydrogen peroxide permeates the protein at a fluctuating entry point located between Phe-68 and Phe-142. Conformational fluctuations of these Phe residues determine the accessibility of hydrogen peroxide to the interior [10]. In comparing the two simulated models, the average distance between the backbones of Phe-68 and Phe-142 decreases from 12.0860.95 Å in n-HRP to 10.6460.54 Å in p-HRP (Fig. 6). Accordingly, the average distance between their side chains decreases from 12.1961.87 Å to 10.2361.21 Å . Examination of the RMS fluctuations of these residues also reveals significant decreases in p-HRP. The RMSF of Phe-68 and Phe-142 has declined from 1.80 and 1.22 Å in n-HRP to 1.08 and 0.84 Å in p-HRP, respectively. These values indicate that the bottleneck entry has become tighter and less flexible upon substitution. These effects are expected to be less sensible for small peroxides and significant for bulky ones.
In the proposed mechanism for the reaction of H 2 O 2 with the active site of HRP two essential features are acid-base catalysis by the distal histidine (His-42) and charge stabilization of a precursor enzyme-substrate complex by the conserved distal arginine (Arg-38) [50,51]. The histidine is thought to facilitate formation of the initial iron-peroxide complex by deprotonating the peroxide and subsequently promoting cleavage of the oxygen-oxygen bond by protonating the distal oxygen [50]. Comparing the models shows that the surface accessibility of this histidine in p-HRP (8.0764.88 Å 2 ) is significantly less than n-HRP (23.9569.62 Å 2 ). This change could lead to decreased reactivity of the modified enzyme to peroxide substrates.   [2,13]. Reduction of compound II to the native enzyme is the rate limiting step in the catalytic cycle of HRP [52,53]. To react with compound II, the aromatic reducing substrate should diffuse from the bulk solution to the aromatic-binding pocket (i.e. the heme cavity). Substrate oxidation by HRP C occurs at the exposed heme edge, a region comprising the heme methyl C18 and heme meso C20 protons [2]. Spectroscopic and crystallographic studies have revealed a detailed picture of the site where aromatic substrates bind and react with the protein. A ring of three peripheral Phe residues 68, 142, and 179 guards the entrance to the exposed heme edge (Fig. 6). Amino acid residues Phe-68, Gly-69, Leu-138, Pro-139, Ala-140, Pro-141, Phe-142, and Phe-179 are flanking the substrate access channel and together with the heme C20-and C18-methyl groups form the aromatic-binding pocket of HRP [5]. Although in some cases hydrogen bonding between the reducing substrate and the active site residues of the distal heme pocket  contributes to the stability of the substrate-HRP complex, most HRP substrates do not possess the potential to make such interactions and will therefore depend more on the hydrophobic interactions which characterize the peripheral region of the substrate channel of HRP [54]. Accordingly, besides the size of reducing substrates, two determinant factors could affect the reaction rate; volume of the active site cavity which controls the penetration of reducing substrates into the active site and the extent of the hydrophobic patch functioning as a binding site or trap for reducing aromatic substrates which is related to the enzyme affinity for the aromatic substrates. Our results show that both factors are affected upon the substitutions. As mentioned before, comparison of the RMSD per residue profiles of the MD trajectories revealed that in p-HRP, the C terminal part of the loop D9E (residues 140-144) and the helix E (residues 145-153) undergo more conformational changes along the simulation (data not shown). Phe-142 of the loop D9E is located at the opening of the heme crevice and together with other residues of this loop (Leu-138, Pro-139, Ala-140 and Pro-141) covers one side of the inner wall of substrate access channel. Upon substitutions, conformational changes which occur in the loop D9E push these residues and the subsequent helix E toward the inside of the heme cavity, making the substrate access channel and its opening narrower in p-HRP. At the same time, the strand b1 (residues 174-176), which is located at the proximal side of the heme cavity and holds one of the substituted residues, also moves toward the inside of the heme crevice, so as to tighten it even more. Together, these movements make the active site cavity smaller in p-HRP. To evaluate this effect, average structure of each model was submitted to CASTp server (Computer Atlas of Surface Topology of protein) and the volume of the active site cavity was measured [55]. CASTp calculations showed that the area and volume of the active site cavity (i.e. the aromatic substrate binding pocket) are decreased from 675.2 Å 2 and 960.2 Å 3 in n-HRP to 441.4 Å 2 and 591.3 Å 3 in p-HRP, respectively. These changes are expected to hinder the access of bulky aromatic substrates to the heme active site and so affect the rate limiting step in the catalytic cycle of HRP. In agreement, the values in Table 6 show that hydrophobic residues forming the substrate-binding pocket in HRP are generally less exposed in p-HRP, implying that the hydrophobic patch functioning as a binding site or trap for reducing aromatic substrates is shrunk in the mutated enzyme. Values in Table 6 also indicate that the accessibility of the exposed heme edge has been significantly reduced after mutations. Summing up the above simulation results, the conclusion is reached that such conformational changes at the aromatic substrate-binding site of HRP may contribute in reducing the affinity and reactivity of the enzyme for aromatic substrates. Again, the effects are expected to be less sensible for small substrates and significant for bulky ones.

Discussion
Directed mutagenesis and chemical modification have been frequently used to characterize or stabilize HRP. HRP directed mutagenesis has been applied primarily to identify those amino acids which are critical for catalysis, structure maintenance, or stability of the enzyme. A detailed description of the introduced point mutations and their effects can be found in reviews [2,15,56,57]. However, genetic engineering tools have also been used to improve HRP stability. One of most extensive studies has been carried out by Ryan and Ó 'Fágáin, who described 22 HRP mutants to evaluate their stability against hydrogen peroxide. Most mutants displayed little or no alteration in H 2 O 2 stability, while three mutants (K232N, K241F and T110V) exhibited significantly increased H 2 O 2 tolerance, among them, the K241F single point mutation enhanced hydrogen peroxide tolerance about 12-fold [26]. In another study, they analyzed the effects of 13 single and 3 double point mutations on the stability of HRP. Three single mutants (K232N, K232F, K241N) demonstrated increased stability against heat (up to 2-fold) and solvents (up to 4-fold) [27]. These studies indicate that single substitution of either K232 or K241 by Phe residues is beneficial for the stability of HRP structure, although it cannot be inferred that simultaneous mutation of these two is additive or synergistic. There are also some reports indicating that modification of solvent-exposed lysines 174, 232 and 241 improves HRP stability [20,58], while exhaustive modification of all six lysines of HRP leads to decreased stability [18]. Accordingly, in this study, we replaced these surface exposed lysines of HRP by phenylalanine, in order to evaluate their simultaneous effects on HRP structure. In agreement to the aforementioned experiments, our simulations provide some structural clues implying that simultaneous substitution of these three lysines by phenylalanine residues can improve the enzyme stability.
Genetic engineering tools are not the only way to modify HRP properties. Actually, chemical modifications have been used more frequently to stabilize HRP. Clearly, chemical modification of solvent-exposed lysines 174, 232 and 241 has produced the best results. Chemical modification of these lysines, ranging from the use of a cross-linker through attachment of polyethylene glycol to simple acetylation, has succeeded in stabilizing HRP to varying degrees [18][19][20][21][22][23][24][25][28][29][30][31][32]. A comprehensive review has been published covering the structure and physicochemical properties of the modifiers and their effects on HRP stability [59], although more recent studies are also available [28,31,32]. Perhaps, one of the first and most extensive studies is reported by Ugarova et al., who studied thermostability of HRP after modification of its lysyl amino groups with a variety of modifiers including anhydrides of monocarboxylic and dicarboxylic acids as well as picryl sulfonic acid. Some of these compounds neutralized the positive charge on the lysine, whereas others reversed it. Authors concluded that modification of three out of six lysine residues produces the major effect on macromolecular conformation and thermostability of HRP [18].
In the case of this study, three accessible charged lysine residues were replaced by hydrophobic phenylalanine residues. So, one can expect the enzyme surface to be more hydrophobic after substitution, but our simulation results disagree with this speculation. Both the hydrophobic and hydrophilic exposed areas are harmoniously declined in p-HRP, so as the ratio between the areas remains almost identical after substitution. Consistent with the total surface area, the average values of radius of gyration and the total volume of protein also decreased upon substitution (Table 1). Although the magnitude of these changes is not statistically significant, they show a similar trend of decrease along with increasing the compactness of protein structure.
Further analysis of the trajectories seeking the local effects of substitutions led us to several clues in the intramolecular hydrogen bonding network indicating stabilization of p-HRP; (1) Analysis of the internal protein hydrogen bonds evinces that the number of stable hydrogen bonds has increased in p-HRP. For example, p-HRP benefits from 12 more intramolecular hydrogen bonds with existence percentage greater than 80% as compared to n-HRP. (2) Hydrogen bonds forming the b-sheet (i.e. between the strands b1 and b2) and helix G in p-HRP are more persistent than those of n-HRP. (3) We also found a persistent new hydrogen bond between the protein moiety (Phe-174) and the heme prosthetic group (O2D propionate oxygen), as well as improvements in the probability of some pre-existing hydrogen bonds in p-HRP. This additional hydrogen bond could strengthen the heme-protein interaction. It is important because in plant and fungal peroxidases, the binding of the heme prosthetic group to the protein moiety greatly depends on non-covalent interactions. (4) The hydrogen bonding strength between some residues involved in proximal Ca 2+ coordination shows considerable improvements upon substitutions (Table 4). Such changes could enhance the stability of proximal Ca 2+ binding pocket which is assumed to be essential for the structural and functional integrity of the enzyme [27,60]. (5) Formation of two persistent stitching hydrogen bonds between the connecting loops GH and F9F0 could enhance the integrity of this part of p-HRP structure. Although the contribution of these changes to the protein stability may not be great individually, the additive contributions of their small favorable enthalpic changes can produce an overall considerable enthalpic gain in stabilization of modified HRP.
In addition to the hydrogen bonding network, the results revealed subtle differences between the secondary structure content of the models. DSSP analysis of the MD trajectories showed that the helices F9 and G, as well as the b-sheet structure are more persistent in p-HRP. The structural stability of b-sheet in the modified model is in accordance with higher propensity of Phe for b-sheet conformation as compared to lysine [61,62]. The end of the helix A is also extended by three residues in p-HRP. A similar trend was also observed in the helix D during the first 5 ns of analyzed time frames.
These substitutions also decrease the RMSD and RMSF values. Such observation has also been reported in experimentally stabilized derivatives of HRP. Ugarova et al. observed that in successfully stabilized derivatives, the intensity of the Soret band in CD spectra of HRP decreases [18]. Accordingly, they inferred that thermostability of the modified enzymes increases due to the decreased conformational mobility of the protein backbone around the heme. Similar changes in CD spectra of HRP after modification with maleic anhydride, citraconic anhydride and anthraquinone 2-carboxylic acid has also been reported [24,30]. Our simulation results provide additional support for this inference, but interestingly indicate that this phenomenon is not merely limited to the heme pocket and the overall protein structure experiences a similar reduction in flexibility upon substitutions.
Mutagenesis and chemical modification of amino acids can also influence catalytic properties. To our knowledge, surface-exposed lysines of HRP were always modified or substituted with the aim of improving stability and in most reports, the activity of the stabilized forms of HRP has been described to be similar-to-wild type or with minor changes (see [59] and the references within), sometimes without making references to detailed kinetic parameters. Surprisingly, in some cases, precise evaluation of the reported data led us to a confusing conclusion. Actually, it seems that in modified HRPs, activity has attracted less attention than stability. This may arise from the fact that HRP has a very high catalytic turnover and some loss in this activity does not significantly affect its practical usefulness. For example, Ugarova et al. studied the effect of 12 charge-neutralizing and chargereversing lysine modifiers on thermostability of HRP. Comparison of the temperature-activity profiles indicated that at any temperature between 20-55uC, the modified HRPs exhibit lower activity than the native enzyme [18]. Considering their experimental measurements, the following conclusions can be inferred: (1) According to the reported temperature-activity profiles, at about 50uC, which is the temperature of HRP maximum activity, modified HRPs show about 20% less activity in comparison to the native enzyme.
(2) These profiles also indicate that thermostabilized HRPs require more heat (up to 15uC) to achieve their maximum activity. (3) Their CD spectra imply that the conformational mobility has decreased and so they inferred that the protein moiety around the heme has become more rigid. Taken together, the conclusion is reached that upon stabilization, the enzyme becomes more rigid, loses some of its activity and needs more heat to acquire the flexibility required for similar to wild-type enzyme activity level. In agreement to these experimental observations, our simulation results showed that upon chargeneutralizing substitutions, the mobility of the protein structure as well as the volume of the reducing substrate binding cavity decreases and the corresponding access channel becomes narrower. Accordingly, lower affinity of the enzyme for reducing substrate is expected. Experimental measurements of the apparent Michaelis constant (i.e. the K m value) for both K232F and K241F single mutants are also in line with this conclusion, although the apparent k cat measurements indicate some increases [27]. Since the number of substitutions in these experimental studies and our simulations are not the same, some differences between their outcomes are not unexpected. In support of this conclusion, it is noteworthy to mention that Ugarova et al. has shown that the number of lysine modification produces a major effect on the macromolecular conformation [18].
Among the charge-neutralizing lysine modifications, it has been reported that modification of HRP with glucosamine hydrochloride [63] and anthraquinone 2-carboxylic acid [29,30] generates not only more stable but also catalytically more active enzymes. In the case of chemical modifications, it is difficult to compare the results of these studies with our Lys-to-Phe substitutions, because the structure and physicochemical properties of these modifiers differ significantly from that of Phe side chain. Glucosamine is a highly hydrophilic carbohydrate having several hydroxyl groups [63]. Attachment of anthraquinone 2-carboxylic acid (a tricyclic aromatic electron relay) to the side chain of lysine residue forms a bulky side chain which differs greatly from that of Phe in size and hydrophobicity [29,30]. Such differences in structure and physicochemical characteristics could result in substantial variations in outcomes.
In conclusion, the results of this simulation study indicate that hydrophobization of the surface exposed lysines of HRP does not necessarily lead to an unfavorable change in the charge/ hydrophobicity balance at the protein surface and may produce stabilized derivatives with acceptable changes in enzyme activity. We have no evidence to think that this phenomenon is merely limited to the surface exposed lysines. Therefore, it would be a good idea to check whether substitution of phenylalanine for the other surface-exposed charged residues (apart from lysines) is beneficial for HRP stability. More generally, the results could also provide possible structural clues for understanding how manipulation of charge/hydrophobicity balance at the protein surface could be translated into structural and functional changes within a protein.