Stability Mechanisms of a Thermophilic Laccase Probed by Molecular Dynamics

Laccases are highly stable, industrially important enzymes capable of oxidizing a large range of substrates. Causes for their stability are, as for other proteins, poorly understood. In this work, multiple-seed molecular dynamics (MD) was applied to a Trametes versicolor laccase in response to variable ionic strengths, temperatures, and glycosylation status. Near-physiological conditions provided excellent agreement with the crystal structure (average RMSD ∼0.92 Å) and residual agreement with experimental B-factors. The persistence of backbone hydrogen bonds was identified as a key descriptor of structural response to environment, whereas solvent-accessibility, radius of gyration, and fluctuations were only locally relevant. Backbone hydrogen bonds decreased systematically with temperature in all simulations (∼9 per 50 K), probing structural changes associated with enthalpy-entropy compensation. Approaching T opt (∼350 K) from 300 K, this change correlated with a beginning “unzipping” of critical β-sheets. 0 M ionic strength triggered partial denucleation of the C-terminal (known experimentally to be sensitive) at 400 K, suggesting a general salt stabilization effect. In contrast, F− (but not Cl−) specifically impaired secondary structure by formation of strong hydrogen bonds with backbone NH, providing a mechanism for experimentally observed small anion destabilization, potentially remedied by site-directed mutagenesis at critical intrusion sites. N-glycosylation was found to support structural integrity by increasing persistent backbone hydrogen bonds by ∼4 across simulations, mainly via prevention of F− intrusion. Hydrogen-bond loss in distinct loop regions and ends of critical β-sheets suggest potential strategies for laboratory optimization of these industrially important enzymes.


Extended Reference Simulations
As discussed in the main text, the 3 ns NVT reference simulation (glycosylated protein, 0.3 M NaCl, 300 K) was extended by 20 ns in four additional simulations with different seeds. Similarly, we extended the following simulations to 20 ns: TvL without glycosylation in 0.3 M NaCl at 300 K, TvL with glycosylation in 0.3 M NaCl at 400 K, and TvL without glycosylation in 0.3 M NaCl at 400 K. The backbone RMSD curves for these simulations are shown in Figure S3a, S3b, and S3c, respectively.
The RMSD curve for the non-glycosylated protein at 300 K in 0.3 M NaCl background ( Figure S3a) is well-behaved throughout the entire 20 ns simulation. The RMSD curves for the 400 K simulations of the glycosylated protein ( Figure S3c) and non-glycosylated protein ( Figure S3d) both shows substantial increases after ~3 -5 ns in agreement with the magnitude of the temperature perturbation.    Table S2. Glycosylated TvLαin NaCl background: Salt bridges found in 3 ns NVT simulations and their persistence (Cons. %) in percentage of 100 frames sampled from 2 last ns of simulations. The simulation average distance (r) in Å between the center of mass of interacting charged groups is also given. Color codes indicate stable salt bridges found in the majority of simulations.  Table S3. Non-glycosylated TvLαin NaCl background: Salt bridges found in 3 ns NVT simulations and their persistence (Cons. %) in percentage of 100 frames sampled from 2 last ns of simulations. The simulation average distance (r) in Å between the center of mass of interacting charged groups is also given.
acceptor donor location VAL27 VAL30 loop, exposed SER33 ARG121 β-sheet start ASP42 LYS39 turn, exposed GLN45 PRO4 first β-sheet, exposed THR51 VAL10 β-sheet end, exposed HIE66 TRP107 HIE66 involving T3 ASP118 THR114 helical segment, burried ALA134 ASP131 loop, exposed THR180 SER177 loop, exposed ALA241 SER202 turn, burried pointing to T3 site ASN304 ILE301 loop/α-helical segment, exposed ASP364 THR361 loop, exposed LEU365 THR361 loop, exposed SER370 PRO367 turn, exposed Table S8. Simulation averaged electrostatic interaction energy between the residue pairs listed in bold in Table S6 and S7 and persistence of backbone hydrogen bonds for the same residues. The analysis was carried out on the last 2 ns of the 300 K, 350K, and 400 K simulations of glycosylated TvL in 0.3M NaCl. Hydrogen bonds in secondary structure and loop regions are indicated in red and green, respectively.  Figure S51. Backbone hydrogen bond persistence (a) and electrostatic interaction (b) between the involved residues, averaged across the labile hydrogen bond pairs marked in Table S6 and Table S7. Standard deviations are indicated with error-bars. Figure S52. Correlation between backbone hydrogen bond persistence (%) and electrostatic energy (kcal/mol) evaluated between the entire residues for (a) residue pairs in structured parts of the protein (red in Table S8), (b) residue pairs in loosely structured parts of the protein (green in Table S8), and (c) residue pairs in both structured and unstructured parts of the protein. Figure S53. Electrostatic analysis of the C-terminal unfolding observed for glycosylated TvL in zero ionic background at 400 K (b, d, f), but not at 300 K (a, c, e). The time series show electrostatic interactions between three groups: "water" consisting of all TIP3P molecules, "helix" consisting of the last C-terminal residues 489 -499, and "protein_nohelix" consisting of the remainder of the protein.

Persistence of Backbone Hydrogen Bonds in Extended Simulations
As discussed in the article and above (under Extended Reference Simulations), extended simulations were made in 0.3 M NaCl background for TvL at 300 K with glycosylation (NAG_0P3M_NACL_300K), and without glycosylation K (noNAG_0P3M_NACL_300K). The corresponding simulations at 400 K with and without glycosylation (NAG_0P3M_NACL_400K and noNAG_0P3M_NACL_400K) were also extended.
The number of persistent backbone hydrogen bonds (HB) was calculated for the 500 MD snapshots between t = 10 ns and t = 20 ns in each extended molecular dynamics trajectory (Table S9). With reference to the RMSD plot for the simulations (Figure 1 in the article main text and Figure S3 in Supporting Information), it is seen that well-behaved RMSD curves are associated with a larger number of persistent HB. Thus it is only reasonable to compare HB numbers from simulations with equally well-behaved RMSD curves.
The extended noNAG_0P3M_NACL_300K simulation yields 164 HB. However, in three out of four simulations with new seeds, the extended NAG_0P3M_NACL_300K simulations have more HB (165 -167) than noNAG. This is in agreement with the numbers of persistent HB in the 3 ns NVT simulations for NAG (165) and noNAG (162).
The temperature effect is substantial also in the extended simulations: At 400 K, NAG and noNAG have 142 and 144 HB, respectively. In the original 3 ns NVT simulations NAG and noNAG had 148 and 147 HB, respectively. Thus the longer simulation time has decreased the number of persistent hydrogen bonds by ~5 in the high temperature case.
In conclusion, the analysis of hydrogen persistence for the longer simulations from new seeds agrees with the major conclusions from the analysis of 3 ns NVT simulations. In particular, the increase in temperature from 300 K to 400 K is associated with a marked reduction of persistence HB, whereas the presence or absence of NAG has a more subtle effect. Table S9. Persistent backbone hydrogen bonds (HB) calculated for the 500 MD snapshots between t = 10 ns and t = 20 ns in extended NVT simulations. NAG_0P3M_NACL_300K_Seed1  155  NAG_0P3M_NACL_300K_Seed2  167  NAG_0P3M_NACL_300K_Seed3  170  NAG_0P3M_NACL_300K_Seed4  165  noNAG_0P3M_NACL_300K  164  NAG_0P3M_NACL_400K  142  noNAG_0P3M_NACL_400K 144