Ion–ion interactions in the denatured state contribute to the stabilization of CutA1 proteins

In order to elucidate features of the denatured state ensembles that exist in equilibrium with the native state under physiological conditions, we performed 1.4-μs molecular dynamics (MD) simulations at 400 K and 450 K using the monomer subunits of three CutA1 mutants from Escherichia coli: an SH-free mutant (Ec0SH) with denaturation temperature (Td) = 85.6 °C, a hydrophobic mutant (Ec0VV) with Td = 113.3 °C, and an ionic mutant (Ec0VV_6) with Td = 136.8 °C. The occupancy of salt bridges by the six substituted charged residues in Ec0VV_6 was 140.1% at 300 K and 89.5% at 450 K, indicating that even in the denatured state, salt bridge occupancy was high, approximately 60% of that at 300 K. From these results, we can infer that proteins from hyperthermophiles with a high ratio of charged residues are stabilized by a decrease in conformational entropy due to ion–ion interactions in the denatured state. The mechanism must be comparable to the stabilization conferred by disulfide bonds within a protein. This suggests that introduction of charged residues, to promote formation of salt bridges in the denatured state, would be a simple way to rationally design stability-enhanced mutants.

The unusually high ratio of charged residues in hyperthermophile proteins might contribute to protein stability through an effect on the denatured state ensemble. We performed molecular dynamics (MD) simulation [25][26][27][28][29][30][31] to elucidate the characteristics of charged residues in the denatured state, using three E. coli CutA1 mutants: an SH-free mutant with T d = 85.6 °C (Ec0SH), a hydrophobic mutant also lacking SH groups with T d = 113.3 °C (Ec0VV), and an ionic mutant with T d = 136.8 °C (Ec0VV_6), as previously reported 23 . The hydrophobic mutant Ec0VV differs from Ec0SH by only two residues and has a T d 27.7 °C higher than that of Ec0SH. The ionic mutant Ec0VV_6 differs from Ec0VV by the addition of six charged residues and has a T d 23.5 °C higher than that of Ec0VV.
The tertiary structure of EcCutA1 clearly resembles that of PhCutA1. The monomer structure consists of three α-helices and five β-strands, and three monomers are assembled into a trimer through interactions between the edges of three β-strands (Fig. 1A). We performed 1.4-μs MD simulations at 400 K and 450 K for three monomers (112 residues) of the described EcCutA1 mutants, and the structure ensembles in the heat-denatured state were evaluated by examining changes in features of the structure due to substitution of residues. The average occupancy of salt bridges by the six substituted residues of Ec0VV_6 was 89.5% at 450 K, indicating that about 90% of each charged residue forms salt bridges even in the heat-denatured ensemble. It appears that stabilization due to salt bridges in the native state is largely cancelled when salt bridges form in the denatured state. However, this interaction might suppress the flexibility of the protein in the denatured states, contributing to protein stability by decreasing entropy in the denatured state, similar to the mechanism responsible for stabilization of protein disulfide bonds 18 .

Results
We performed 1.4-μs MD simulations at 400 K and 450 K for each monomer (subunits A, B, and C) of three mutants (Ec0SH, Ec0VV, and Ec0VV_6). As a reference, 0.4-μs MD simulations at 300 K were also done for the native trimer structures of the three mutants. We found that all three mutant proteins are largely disordered and the average root mean square deviations (RMSDs) of all the Cα atoms reached a plateau after 200-400 ns simulations, at 400 K and 450 K (Fig. 2). The trajectories of the average radius of gyration for three mutant subunits are also shown in Fig. 3. The decrease in radius of gyration for all mutant monomers might be due to changes in structure, from ellipsoid in the native trimer (Fig. 1B) to globular in the disordered monomer (Fig. 1C). The average values of RMSD and radius of gyration are listed in Table S1. The average values of RMSD for each subunit at 450 K were greater than those at 400 K, suggesting that the structures at 450 K are more disordered.
At 450 K, the α-helical structures of all three mutants were nearly destroyed after 0.4 μs ( Figure S1A, Table 1), but about half of the β-sheets were sustained during 1.4 μs (Table 1). By contrast, at 400 K, all three mutants retained approximately 40% of their α-helical structures during 1.4-μs simulation ( Figure S1C, Table 1). Most of the α-helical residues sustained between 1.2-1.4 μs were seen in the α-1 helix in the trimer structure at 300 K (Fig. 4). Over the interval of 0.4-1.4 µs, average values of β-sheet for the three mutants were similar at 300 K and 400 K, but reduced by approximately half at 450 K ( Table 1).
The energy of hydrophobic interaction (ΔG HP ) for stabilization of a protein can be calculated from differences in the accessible surface area (ASA) of hydrophobic and hydrophilic atoms between the native and the denatured states, using the following equation 32,33 : HP non polar polar Here, ΔASA non-polar and ΔASA polar represent the differences in ASA of non-polar and polar atoms of all residues, respectively, between the native and the denatured states.  Table 3). The structure obtained from MD simulations at 300 K was assumed to represent the native state. The differences in hydrophobic energies    Table 3. Calculation of hydrophobic energy assuming that the ASA values at 400 K or 450 K represent those of the denatured states. Calculations were performed using the ΔASA values listed in Table 2 using equation (1). All values are given in kJ/mol. Subscript 1 indicates the calculated difference between the listed mutant and Ec0SH; subscript 2 indicates the difference between the listed mutant and Ec0VV. between Ec0VV and Ec0SH were estimated to be 20.8 kJ/mol and 18.2 kJ/mol at 400 K and 450 K, respectively. These results quantitatively indicate how the hydrophobic mutant, Ec0VV, is stabilized by hydrophobic interactions as compared with its template mutant, Ec0SH. On the other hand, the calculated hydrophobic energies of the ionic mutant substituted with six charged residues, Ec0VV_6, were approximately equal to those of its template mutant, Ec0VV (Table 3), indicating that the alkyl groups of charged residues do not contribute significantly to stabilization due to hydrophobic interaction. Two charged residues may form a salt bridge when the distance between Cε or Cδ of a positively charged residue (Lys or Arg) and Cδ or Cγ of a negatively charged residue (Glu or Asp) is less than 0.6 nm 27 . Figure 5 shows the typical trajectory of changes in distance between Cε of Lys87 and Cγ of Asp39 in Ec0VV_6B during 1.4-μs simulations at 450 K. Lys87 can probably make a salt bridge with Asp39 during this simulation, although the distance between the two residues is greater than 2.5 nm in the native structure. The interaction frequency of Lys87 and Arg88 with other favorable negatively charged residues is shown in Figs 6A and B, respectively. These individual interactions with Lys87 and Arg88 are also shown separately in Figures S2 and S3. As the figures show, each favorable salt bridge fluctuated drastically during simulations, but distance is less than 0.6 nm (suitable for formation of salt bridges) for a considerable portion of the time in both plots. As shown in Fig. 7, approximately one salt bridge with Arg88 (or Lys87) was preserved during MD simulations even at 450 K, even though favorable pair residues with Arg88 (or Lys87) frequently replace during simulations. The average number of salt bridges with Arg88 and Lys 87 was 1.12 ± 0.45 and 1.15 ± 0.53, respectively, during 1.4-μs MD simulations at 450 K, respectively. Furthermore, the occupancy of salt bridges by targeted charged residues during MD simulations could be calculated from the trajectory of the changes in distance between atoms of favorable ion pairs (Fig. 5). Table 4 shows the occupancies of salt bridges by the six substituted residues of Ec0VV_6 during MD simulation at  Figure S2 and Figure S3. 450 K, 400 K, and 300 K. An occupancy value of more than 100% indicates that a charged residue forms more than one salt bridge. The occupancies of salt bridges by the six substituted residues of Ec0VV_6 were 72.9-120.9% and 23.2-130.3% at 450 K and 400 K, respectively, compared to 61.7-285.4% at 300 K. The average percent occupancies were 89.5%, 84.3%, and 140.1% at 450 K, 400 K, and 300 K, respectively; thus, the average occupancy of salt bridges at 450 K was approximately 60% of the occupancy at 300 K. The data in Table 4 also shows how the pairs of native-state salt bridges at 300 K were replaced at 400 K and 450 K during simulations. These results suggest that many charged residues in proteins from hyperthermophiles also form a high percentage of salt bridges in the denatured state.

Discussion
There are several reports 4-6,34,35 focusing on the denatured state ensembles that are in equilibrium with the native state under physiological conditions. Denatured state ensembles may be considerably compact, similar to the molten globule state, containing secondary structures 36 that are not suited to the random-coil model 37 . Denatured structures of pyrrolidone carboxyl peptidase under physiological conditions have been reported to form non-native hydrophobic clusters and native-like helices 5,6 . In the denatured state, charged residues are also involved in strong electrostatic interactions 3,34,38 .
In the current study, we performed 1.4-μs MD simulations of EcCutA1 stability-enhanced mutants at 400 K and 450 K to elucidate characteristics of denatured state ensembles. For the Ec0VV_6 mutant, at 450 K nearly all α-helices and β-sheets disappeared transiently at 500 ns (Fig. 1C), but on average 20% of the total residues maintained β-sheets during the period of 0.4-1.4 μs (Table 1), although the RMSD suggests that the conformation of Ec0VV_6 is drastically disordered (Table S1). At 400 K, one third of the α-helical content of Ec0VV_6 was maintained even after 1.0-1.4 μs simulation (Table 1).
Proteins from hyperthermophiles contain considerably more charged amino acid residues than those of mesophiles 8 , and ion-ion interactions are believed to stabilize the proteins to permit growth at high temperatures 39 . The CutA1 protein from Pyrococcus horikoshii (PhCutA1) has many charged residues and is very stable, with a denaturation temperature of nearly 150 °C 20 . Although the electrostatic energy due to charged residues seems to stabilize PhCutA1 in the native state 12 , it has not been shown how these charged residues contribute to stabilization of the protein in the denatured state. Our salt bridge occupancy results at high temperatures suggest that charged residues form considerable numbers of favorable ion-ion interactions in the denatured state; salt bridge occupancies were 89.5% and 84.3% at 450 K and 400 K, respectively. These values were approximately 60% of the salt bridge occupancy value in the native state (300 K). It is difficult to determine whether the simulated ensembles at 400 K or 450 K are more representative of the denatured state ensembles under physiological conditions. However, the occupancies of ion-ion interactions were similar at 400 K and 450 K, although the secondary structure content was different. These results suggest that when sufficient charged residues are present, a considerable number of charged residues can make salt bridges in the denatured state ensemble. Favorable ion-ion interactions in the denatured state would stabilize the denatured conformation 40 , resulting in destabilization of the native protein. On the other hand, favorable ion-ion interactions in the denatured state ensemble could also decrease the entropy of the conformation in the denatured state, resulting in stabilization of the native protein.
When one or two disulfide bonds are broken, a protein is substantially destabilized due to an increase in entropy in the denatured state 18 . In the current study, approximately 60% of the ion-ion interactions present at 300 K were maintained even at 450 K, suggesting that proteins from hyperthermophiles with a high ratio of charged residues are stabilized by decrease in conformational entropy due to ion-ion interactions in the denatured state.
Hyperthermophilic proteins differ greatly from mesophilic proteins in their change in heat capacity upon denaturation due to the presence of residual structure in the denatured state ensemble 16,41,42 . This feature results in increased protein stability at higher temperatures due to the broadening of temperature function curves for the Gibbs energy of unfolding. These results suggest that the denatured state conformation of hyperthermophilic proteins is more compact than that of mesophilic ones, indicating the role of entropy in determining stability. By analyzing thermodynamic data for 116 proteins, Sawle and Ghosh 43 have reported that entropic stabilization is responsible for the high melting temperature of proteins observed in hyperthermophiles because the gain in enthalpy upon folding is smaller in hyperthermophiles than in mesophiles, whereas the loss of entropy upon folding is greater in mesophiles than in hyperthermophiles.
To our knowledge, it has not been reported that ion-ion interactions in the denatured state ensemble contribute to protein stability by decreasing conformational entropy, similar to the mechanism of stabilization by disulfide bonds. Introduction of charged residues for the formation of the salt bridges in the denatured state may be a useful tool for engineering proteins with improved conformational stability.
MD simulations were performed using GROMACS software (ver. 4.5.5) 44,45 . The missing atoms in the coordinate file of Ec0SH (PDB ID, 4Y65), which are three N-terminal residues of the B subunit and eight N-terminal residues of the C subunit, were modeled in QUANTA2000 (Accelrys) using the coordinates of N-terminal residues of the A subunit as a reference. The structures of Ec0VV and Ec0VV_6 were modeled using FoldX, based on the structure of Ec0SH. Hydrogen atoms were added to each protein. The models were solvated in water boxes  with a minimum distance of 1.2 nm between the protein and the box. Counter-ions were added to the model to neutralize any net charge. The periodic boundary condition was adopted and the long-range electrostatic interactions were computed using the Particle-Mesh-Eward (PME) method 46 . The GROMOS 43A1 force field and SPC/E water model 47 were employed. The system was weakly coupled to a heat bath by velocity rescaling 48 with a relaxation time of 0.1 ps. A Parrinello-Rahman barostat 49 was used to maintain a pressure constant at 1.0 bar for 300 K and 6.0 bar for 400 K or 450 K with a relaxation time of 0.5 ps. Hydrogen atoms were constrained using LINCS 50 , and MD simulations at 300 K, 400 K, and 450 K were conducted with an integration time step of 1 femtosecond (fs). Energy minimizations were done to remove bad van der Waals contacts. Next, the temperature was raised from 50 K to 300 K in increments of 50 K, with 10,000 integration steps at each temperature and a harmonic constraint of C-alpha atoms. Thereafter, the ensemble was equilibrated through four 100-picosecond (ps) cycles with gradually released harmonic constraints: 1000, 100, 10, and 1 kJ mol −1 nm −2 . The subsequent MD stages for the EcCutA1 mutants were carried out without any restraint at 300 K. When the system temperature was increased to 400 K or 450 K from 300 K, pressure coupling was not set during 1000 ps at 400 K or 450 K. The obtained MD trajectories were analyzed using GROMACS software. The calculations for RMSD of Cα atoms and the radius of gyration were performed using the commands 'gmx rms' and 'gmx gyrate, ' respectively. For the salt bridges, 'gmx saltbr' was used with the option t = 0.4 nm, which means that groups that were never closer than this distance were not plotted. For the calculation of ASA values of each atom, the results of 'gmx sasa' (atomarea.xvg: average area per atom) were used. For the trajectory of secondary structures, the command 'do_dssp' was used. For the average α-helix at each residue, the command 'g_helix' was used.