Unfolding Transitions of Peripheral Subunit Binding Domains Show

a Department of Chemical Sciences, Indian Institute of Science Education and Research ... c Center of Computational Natural Sciences and Bioinformatic...
0 downloads 0 Views 2MB Size
Subscriber access provided by UNIV AUTONOMA DE COAHUILA UADEC

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Unfolding Transitions of Peripheral Subunit Binding Domains Show Cooperative Behavior Monika Sharma, Gopalakrishnan Bulusu, and Abhijit Mitra J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.9b01114 • Publication Date (Web): 08 Apr 2019 Downloaded from http://pubs.acs.org on April 8, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Unfolding Transitions of Peripheral Subunit Binding Domains Show Cooperative Behavior Monika Sharmaa*, Gopalakrishnan Bulusub,c, and Abhijit Mitrac. Affiliations: a Department of Chemical Sciences, Indian Institute of Science Education and Research (IISER), Sector 81, Knowledge City, SAS Nagar, Punjab, India. b TCS Innovation Labs – Hyderabad (Life Sciences Division), Tata Consultancy Services Limited, Hyderabad 500081, India c Center of Computational Natural Sciences and Bioinformatics (CCNSB), International Institute of Information Technology (IIIT), Hyderabad, India. *Corresponding author: Monika Sharma, Department of Chemical Sciences, Indian Institute of Science and Education Research Mohali, Sector 81, Knowledge City, SAS Nagar, Punjab, India. Email: [email protected] Tel: +91 7087543458 Fax: +91-172-2240266

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 25

ABSTRACT Characterization of native, intermediate and denatured states is crucial for understanding the factors influencing the stability of proteins. We have carried out molecular dynamics simulations to study the unfolding of three peripheral subunit binding domains (PSBDs): E. coli BBL, Bacillus stearothermophilus E3BD, and human hbSBD, at three different temperatures: 300 K, 330 K, and 400 K; and in the presence of two solvents: water and 5 M guanidinium hydrochloride (GndCl) solution. These proteins share similar folds, with two parallel helices, maintained via a hydrophobic core comprising of residues from their interconnecting loop. BBL is more sensitive to thermal and chemical denaturation in comparison to hbSBD, and E3BD is the most stable of all the three proteins. Effect of temperature on stability of these proteins is more pronounced in ‘water only’ simulations compared to that in the presence of guanidium hydrochloride in high concentrations. Our results show cooperative unfolding transitions of these proteins, which are triggered by an initial melting of the C-terminal helix H2. The consequent loss of interhelical interactions or native contacts, as observed, leads to the subsequent melting of the N-terminal helix H1.

2

ACS Paragon Plus Environment

Page 3 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

INTRODUCTION The factors associated with the stability of compactly folded globular proteins can be identified by investigating their unfolding process. Capturing the events, occurring along their respective denaturation pathways, is crucial for understanding various aspects such as, loss of their enzymatic activities at higher temperatures or in presence of chemicals, which are of utmost importance for their biotechnological applications. In this study, we focus on three peripheral subunit binding domains, which form integral parts of the ubiquitous 2-oxoacid dehydrogenases (ODHc): Escherichia coli peripheral subunit binding domain of 2-oxoglutarate dehydrogenase (OGDH)1 (hereinafter referred to as BBL), Bacillus stearothermophilus peripheral subunit binding domain of dihydrolipoamide acetyltransferase from pyruvate dehydrogenase multienzyme complex (PDH)2 (hereinafter referred to as E3BD), and human subunit binding domain of human branched-chain -ketoacid dehydrogenase complex (BCKDC)3 (hereinafter referred to as hbSBD). These ODHc dehydrogenases constitute a family of very large multienzyme complexes consisting of multiple copies of at least three enzymes: E1 (termed as E1o, E1p, and E1b in OGDH, PDH, and BCKDC, respectively), dihydrolipoamide acetyl, succinyl, and branched-chain transferase E2 (E2o, E2p and E2b, respectively) and dihydrolipoamide dehydrogenase E3. The third component E3 catalyzes the same reaction and is similar in all the three. All ODHc complexes have a structural core of multiple copies of E2 subunits, with E1 and E3 subunits bound on the periphery. The E2 component is a tripartite domain comprised of: an N-terminal region consisting of 1-3 tandem repeated lipolyl domains (LD), a C-terminal catalytic domain, and the interconnecting peripheral subunit-binding domain (PSBD) responsible for binding E1/E3 chains. These three domains are separated by long, flexible linker regions allowing movements which enable active coupling. The central subunit binding domains binds either E1 or E3 at any given point of time, but not simultaneously, and the linker allows the subunit binding domain to move freely with respect to the core. These peripheral subunit binding domains (PSBDs) are one of the smallest folding domains. The canonical native structure of these ~35-50 residue long domains comprises of two parallel -helices, connected by an intervening loop, consisting of either five-residue helix turn or a short 310 helix in BBL and E3BD, respectively. The compact structure of the domain is stabilized mainly by hydrophobic interactions, and the interactions between these subunit binding domains and E3 are observed to be mediated by charged side chains, forming an electrostatic zipper. Structural characterization of E3BD with E34 and E15 subunits in B. stearothermophilus revealed that the same set of residues in E3BD are used in forming the two different complexes, suggesting E1 and E3 competing for association with E2 in the assembled multi-enzyme complex. E1 and E3 bind with the subunit binding domain with nearly same high affinity, and specifically via the N-terminal helix 3

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 25

H1. Figure 1 shows the comparative sequence and structures of the protein domains studied here. Though their sequences are highly similar, the thermal stabilities of these proteins are different. Previous studies on the peripheral subunit binding domains have shed light on the folding kinetics of BBL protein and its homologs. E3BD6 and hbSBD7,8 have been proposed to be two-state folders with melting temperatures of 330 K and 318 K, respectively. The high TM for E3BD is a result of the fact that B. stearothermophilus is adapted to grow at temperatures of 50º-80 ºC. However, BBL is observed to be on the cross-roads of non-cooperative downhill-type and classical two-state folder9. Different biophysical techniques, such as DSC, CD and FRET have been reported to suggest a wide range of melting temperatures, ranging from ~295 K to 335 K, and have indicated BBL as a global downhill folder10. Kinetic folding models obtained by comparing temperature jump (T-jump) relaxation methods and varied initial potential of BBL showed it as a folder with marginal barrier11. This study contradicted the conclusions reported in another study, of the two-state cooperative folding kinetics for BBL protein and its homologs12, which predicted melting temperature as 324-329 K. Replica exchange MD simulations13 inferred weakly cooperative folding mechanism for BBL protein where loss of native tertiary structure at low temperature, and complete loss of secondary structure at higher temperature was observed. Overall, these studies highlight the complexity in the folding behavior of BBL9. Temperature jump and continuous-flow fluorescence methods14 indicated folding half-times of E3BD and BBL as 25 s and 3-5 s, respectively. On the other hand, hbSBD folded cooperatively on 100 s scale8.

Figure 1: Sequence and structure comparison of the subunit binding proteins studied. The sequence alignment is predicted by T-coffee multiple sequence alignment software. Residues are colored with respect to alignment score, with red being very good, yellow being average and green depicting less than average. Proteins are shown as helical ribbons with polar residues highlighted as yellow sticks, and hydrophobic residues as grey sticks.

4

ACS Paragon Plus Environment

Page 5 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Besides thermal denaturation, chemical denaturation of these proteins has been studied in presence of guanidinium chloride and urea. Single molecule-FRET data on chemical denaturation of BBL at low temperature showed BBL to be a one-state downhill folder15. Equilibrium chemical denaturation experiments14 carried out in 5 M guanidinium hydrochloride (5 M GndCl) solution showed E3BD to be more stable than BBL. Circular dichroism studies of tryptophan mutants of hbSBD and E3BD in GndCl solutions indicated cooperative folding of these proteins8,14. There have been computational efforts to observe underlying conformational ensembles formed during such studies. For instance, Gomodel simulations of BBL16,17 showed barrier-less folding for pdb 1BBL, in line with the DSC and FRET experiments18. However, similar experiments on different sized constructs of BBL pointed towards a two-state folding mechanism14,19. When atomistic simulations of thermophilic E3BD were carried out at elevated temperatures and the trajectories mapped with heteronuclear NMR experiments to characterize the denatured ensemble, the overall folding pathway appeared to be via nucleationcondensation20. Langevin dynamics simulations of hbSBD7 showed it to be a two-state folder. These studies do capture the folding dynamics of these proteins. However, the atomistic details underlying the denaturation of these proteins with respect to temperature and chemicals are understood only to a limited extent, except for the case of one all-atom replica exchange MD (REMD) on BBL21. We have addressed this by carrying out extensive MD simulations of the three homologs: BBL, E3BD, and hbSBD. For BBL, we used an NMR structure (PDBid: 1BBL) for our study. This 37-residue construct is similar to Munoz construct with broad unfolding temperatures ranging from 295 K to 335 K18. Recently, similar unfolding experiments were carried out to study the effect of temperature and denaturants such as trifluoroethanol (TFE) on the unfolding and conformational stability, flexibility, and dynamics of N-terminal domain of TAR DNA binding protein-43 (TDP-43)22, Cu/Zn Superoxide dismutase (SOD1)23–25, and RNA recognition motifs (RRMs) of TDP-4326. These studies revealed the conformational heterogeneity within the structures related to the distinct structural stabilities and unfolding dynamics of these proteins and shed light on the misfolding and aggregation of these proteins. We have carried out long duration atomistic simulations of the three proteins, in water (often referred to here as ‘water only’ to highlight the absence of denaturant) and in 5 M GndCl solutions, respectively at three different temperatures: 300 K, 330 K and 400 K, with three objectives in mind. They are: to determine their respective stabilities; to investigate the dissimilarities in their unfolding dynamics notwithstanding their structural similarities; and to map the similarities and dissimilarities in the trajectories with their respective sequences. The simulations carried out in 5 M GndCl solution were compared with simulations performed in water to understand the net effect of the denaturant. We used higher concentration of GndCl to overcome the issue of low sensitivity of subunit binding domain to chemical denaturation14. 5

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 25

METHODS Simulations setup Simulations were performed in explicit solvent using Gromacs-2016.527 with Amber ff99SB-ILDN forcefield28. Coordinates for BBL, E3BD, and hbSBD systems were extracted from PDBid: 1BBL1, 1W3D2, and 1ZWV3. For water only simulations, systems were solvated in rectangular TIP3P solvent box and neutralized with Na+ or Cl- ions. Long range interactions were treated using Particle Mesh Ewald (PME) with tolerance of 1e-05, and grid spacing of 1.2 Å. A non-bonded cut-off of 12 Å was applied to the Lennard-Jones potential and a cut-off for the direct-space part of the Coulomb forces with a switching function starting at distance 10 Å and reaching zero at 12 Å. A 2 fs time step was used for all simulations while constraining bonds involving hydrogens with Lincs algorithm. Periodic boundary conditions were employed to eliminate surface effects. To reach the desired temperature, minimized systems were simulated for 10 ns while coupled to v-rescale thermostat implemented in Gromacs. For denaturant solution, Kirkwood-Buff derived forcefield (KBFF) parameters were employed for guanidinium chloride29. BBL, E3BD, and hbSBD systems were solvated in systems containing water/Guanidinium ions/Cl- ions/Na+ ions as 6083/553/620/66, 8570/772/779/2, 7403/667/720/57, respectively. The equilibrated structures were simulated with multiple independent runs: two at 300 K and 330 K, and five at 400 K, with different initial velocity assignments using random seed method. Two independent 0.5 s long simulations were carried out at 300 K and 330 K, and five 100-ns simulations were carried out at 400 K. All the simulations were carried out in NVT ensemble. Trajectory and structural analyses Structural analyses were done using built-in programs of Gromacs-2016.527 and VMD-1.9.330. Root mean square deviations and radius of gyration values for all backbone atoms were calculated with reference to the experimental structure. Root mean square fluctuations (RMSF) of Cα atoms in all the trajectories were also calculated. The changes in secondary structure of protein during the course of simulation were analyzed using DSSP31 program as implemented in Bio3D32 package. Helical propensities of each residue were calculated averaged over all the conformations sampled during the multiple simulation runs as analyzed by DSSP. The radial distribution functions, 𝑔(𝑟) were calculated between solute and solvent atoms using Gromacs. Three measures were used for evaluating inter-residue contacts: hydrogen bonds; Q, the fraction of native contacts and non-native contacts. Hydrogen bonds were calculated using hbond analysis module of Gromacs with cutoff values of 3.5Å and 30° for Donor-Acceptor distances and for Hydrogen-Donor-Acceptor angles, respectively. A contact (native or non-native) is defined when a 6

ACS Paragon Plus Environment

Page 7 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

heavy atom of one residue interacts with another heavy atom of another residue (excluding the adjacent three neighboring residues) within a cutoff distance of 6.0Å. A native contact is the contact that is present in the starting experimental structure, and non-native contacts are new contacts that are formed during simulation. Native contacts and non-native contacts between subunit binding domains were calculated using cpptraj33 after pooling in all the conformations sampled at every 2 ps during different trajectory runs. Dynamic cross correlation analysis The extent of correlated motion between two residues was calculated as the magnitude of the correlation coefficient between Cα atoms of two residues34. The cross-correlation coefficient 𝐶𝑖𝑗 for each pair of Cα atoms of residues 𝑖 and 𝑗 is calculated as, Cij  ri  rj

ri 2

1/2

rj2

1/2

where ∆𝑟𝑖

is the displacement from mean position of atom 𝑖 and the < > symbol represents the time-average. Dynamic cross correlation (DCC) analyses was carried out using Bio3D32 package for Cα atoms by pooling in all the conformations sampled at every 2 ps during multiple simulation runs. Cartesian Principal Component Analysis Essential dynamics or principal component analysis (PCA)35,36 was done using Gromacs-2016.5. A variance-covariance matrix was constructed using the coordinates of all heavy atoms from the multiple trajectories. The matrix was diagonalized to generate a diagonal matrix of eigenvalues and a transformation matrix comprising eigenmodes using covar module. The eigenmodes were analyzed using anaeig module, where the time-dependent projections of the principal components, onto the trajectories, were studied. Calculation of solvent accessibilities Solvent accessibilities values were calculated using sasa algorithm implemented within VMD-1.9.3 for heavy atoms of hydrophobic residues (designated as hydrophobic sasa value), and heavy atoms of hydrophilic residues (designated as hydrophilic sasa value). For hydrophobic sasa, amino acids: alanine, leucine, valine, isoleucine, proline, phenylalanine, methionine and tryptophan were considered, and rest of the amino acids were considered for hydrophilic sasa calculation. Community network analysis Protein structure networks were generated using functionality of Bio3D package implemented in R. The functionality builds network obtained from atomic correlation data derived from normal mode analysis (NMA) and molecular dynamics (MD) simulation. Normal mode analysis of experimental structures was done for coordinates obtained from PDB, and networks were analyzed using dynamic 7

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 25

cross-correlation data generated from normal modes. For simulated trajectories, cross-correlation data calculated using Pearson correlation method was used. Different communities were detected using Girvan-Newman betweenness method37, with cross-correlation cutoff value of 0.35; and all-residue networks were visualized and colored using igraph package installed in R. Node centralities assessing the density of connections per node were calculated as described in 38. Free energy landscapes Free energy landscapes (FEL) were constructed using three sets of reaction coordinates: rmsd versus radius of gyration (Rg); rmsd versus fraction of native contacts; and rmsd versus hydrophobic sasa. Radius of gyration values were calculated using gyrate utility of Gromacs. The values of reaction coordinates along the concatenated trajectories were binned into a two-dimensional histogram using a grid size of 64. For FEL corresponding to rmsd versus fraction of native contacts, grid size of 32 was used. The landscape was then generated by taking the negative natural logarithm of the bin counts at each position: 𝐹𝐸𝐿(𝑅𝑀𝑆𝐷,𝑅𝑔) = ― 𝑘𝐵𝑇𝑙𝑛[𝑃(𝑅𝑀𝑆𝐷,𝑅𝑔)] where 𝑘𝐵and 𝑇 are Boltzmann constant and temperature, respectively; and 𝑃 is the joint probability of structure having values of RMSD and Rg. RESULTS Structural stabilities of subunit binding domains The subunit domains studied contain two parallel helices stabilized by hydrophobic interactions, and an interconnecting loop enclosing a close-packed hydrophobic core (Fig. 1). The protein BBL (pdbid: 1BBL) is a small folded domain of 52 amino acids; with two parallel helices (H1: residues 14-23 and H2: 40-48), two short extended strands (residues 12-13 and 24-26), a well-ordered five-residue helixlike turn (26-30), and an irregular and more disordered loop (residues 31-39). E3BD contains two short parallel helices (H1: residues 133-142 and H2: 160-167), connected by a 310 helix (residues 145149) and a structured loop (residues 150-159). Human homolog hbSBD (pdbid: 1ZWV) contains two parallel helices (H1: residues 12-20, H2: 39-47) with a long irregular loop connecting them. We analyzed changes in root mean square deviations (RMSD), in response to temperature variations and that of presence or absence of denaturant in the aqueous solvent, for the protein backbone atoms of the respective hydrophobic cores of BBL (residues 12-48), E3BD (residues 130-170) and hbSBD (residues 10-48), with respect to their corresponding crystal structure coordinates. Figure 2 shows the time evolution of the backbone atoms RMSD and the radius of gyration values for each simulation run. The general trend, inferred from the observed increase in RMSD values (Fig. 2A), reveals that the increase in simulation temperature, in both water and denaturant solution, results in greater structural deviations. Of all the three, BBL protein showed larger RMSD values at higher 8

ACS Paragon Plus Environment

Page 9 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

temperatures. Thus, assuming that increased structural deviation correlates with reduced thermal stability, we may conclude that the E. coli BBL protein is the least thermostable of the three proteins. In comparison, based on similar considerations, E3BD should be more thermostable, followed by hbSBD. The same trend in RMSD values is also observed in the corresponding chemical denaturation studies, thus supporting the conclusion that E3BD is more stable than hbSBD, which is more stable than BBL. Interestingly, the overall backbone deviations are noticeably less for E3BD in presence of denaturant environment even at elevated temperatures, except for in one simulation run at 330 K. As expected, similar trends are also observed for time dependent variation of Rg values (Fig. 2B) as well as for root mean square fluctuations (RMSF) for each residue averaged over all the conformations (Fig. 3A).

Figure 2: Variation of root mean square deviation (RMSD) values and radius of gyration values of proteins at different temperatures and in presence of different solvents with respect to time. All the simulated trajectories shown in different colors depict different independent runs. 9

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 25

Detailed analysis of the RMSF values showed higher fluctuations at increased temperatures for residues of BBL, both in the presence and absence of denaturant, in the aqueous medium. The next in stability order is human homolog hbSBD, which showed higher fluctuations at elevated temperatures of 400 K in water only simulations. Interestingly, the residues of E3BD show higher RMSF values in water only at 330 K, whereas in 5 M GndCl, the RMSF values appear to be insensitive to temperature variations. For all these proteins, it is observed that residues constituting helices show lesser fluctuations than the residues in connecting loops.

Figure 3 A. RMSF for C atoms of subunit proteins averaged over all the simulation runs. B. Timeaveraged helical propensities for PSBDs. Blue colored lines and bars correspond to 300 K, black colored to 330 K and red colored to 400 K. C. Two-dimensional projections along the first two principal components for simulation trajectories sampled at different temperatures as projected onto the eigenvectors calculated for concatenated trajectory for each subunit binding domain. The red colored stars indicate the projections calculated for coordinates obtained from initial experimental structures. Conformational sampling of protein in response to temperature and denaturant We used the DSSP algorithm to analyze the variation in secondary structures during all the simulation runs. (Figs. S1,S2). Time evolution plots of the secondary structure elements for BBL protein at 300 K, 330 K, and 400 K suggest that unfolding of BBL begins via melting of the C-terminal helix H2 (as indicated by color inclusions in blue colored helical region), which is then followed by melting of the 10

ACS Paragon Plus Environment

Page 11 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

N-terminal helix H1. The loss of helical structure is more pronounced in water than in presence of GndCl solutions. For hbSBD, the helical structure for both helices is maintained at 300 K, whereas at 330 K, we observed that the unfolding of C-terminal helix H2 begins at ~200 ns during one of the simulations, and that the helix H1 remains intact during all the simulations. The elevated temperature unfolding of hbSBD at 400 K, on the other hand, also begins with the initial melting of the helix H2, but is followed by the melting of N-terminal helix H1. For E3BD, the helical propensities are maintained for both the helices at 300 K. Interestingly, the two independent simulation runs of E3BD at 330 K show mutually different behavior. During one run, the helical structure for both the helices is maintained, whereas during another, we observed that unfolding of E3BD commences with partial melting of the C-terminal helix H2 at around ~380ns. This induces structural changes in the interconnecting loop, which finally leads to the complete melting of the N-terminal helix H1 at around ~450 ns, while the helical structure of H2 is still partially retained. It may be mentioned that at 400 K, for E3BD we observed only minor secondary structure transitions for both the helices, but with no major loss of helical structure. We therefore carried out quantitative analyses of retained secondary structure by measuring the helical propensities of each residue averaged over all the conformations. As can be clearly seen from the bar plots (Fig. 3B), the helical propensities of residues in the Cterminal helix H2 decreased rapidly with increase in temperature, compared to those observed for residues in the N-terminal helix H1 in both water only and in GndCl solution. A simulation trajectory consists of fast local motions and slow global motions. However, most features of the trajectory can be described by a surprisingly small number of collective degrees of freedom. We carried out principal component analysis to resolve the conformational space of our simulated systems into essential subspaces by filtering global slow motions from fast local motions. To view the temperature and solvent dependent sampling of the protein conformations along the PC space, we concatenated the trajectories, for each PSBD, generated at different temperatures and in presence of different solvents and calculated the eigenvectors for this combined trajectory. Then we computed two-dimensional projections (Fig. 3C) of the respective trajectories projected onto the first two eigenvectors (PC) computed for the concatenated trajectory. It is quite evident that the conformations sampled different regions in PC space at different temperatures, and more diverse sampling is observed in water only than in GndCl solution. Correlation dynamics in response to environment We calculated the correlation matrix for C position fluctuations of the three proteins for different temperatures in the presence and absence of the denaturant, respectively. Dynamic cross-correlation maps for the two states, representing correlation coefficients calculated as averaged over all the conformations sampled during the simulations, are shown in Figure 4. The whole range of correlation 11

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 25

from -1 to +1 is represented in three ranges: blue color corresponding to positive correlation values ranging from 0.25 to 1; red color corresponding to negative correlation values ranging from -0.25 to -1; and white color corresponding to weak or no-correlation values ranging from -0.25 to +0.25. The extent of correlation or anti-correlation is indicated by variation in the intensity of respective blue or red color. A quick glance at the maps indicates that the motions during the unfolding dynamics are anti-correlated as shown by increase in red colored islands in DCC maps at higher temperatures. Comparative analyses of DCC maps with contact maps (Fig. S3) suggest that the residues forming contacts in the experimental structures retain positively correlated motions as indicated by blue colored islands, especially at 300 K. At elevated temperatures in both water and GndCl solutions for BBL and hbSBD proteins, we observed increase in anti-correlated motions. Mapping these interactions with the trajectories suggests that this effect may be related to protein unfolding. For E3BD, this trend is observed only in water only simulations at 330 K. This can be attributed to one of the two simulation runs where E3BD unfolds with consistent increase in RMSD and radius of gyration values. Free energy landscapes depicting dynamic behavior of proteins We assessed conformational landscapes of the subunit binding proteins in response to temperature and solvent environment. The effective two-dimensional free energy landscapes, as obtained from sampled RMSD and radius of gyration values (Rg), are shown in Fig. 5. Comparison of the FELs showed that proteins in water medium (Fig. 5A) explored larger conformational space at 330 K and 400 K. The large RMSD and radius of gyration values of extracted structures, corresponding to red colored islands representing FEL minima, reveal that all the three proteins are in their respective unfolded states with significant loss of helical structure. At the lower temperature of 300 K, representative structures for minima for BBL, E3BD and hbSBD retained the helical content, with lower rmsd and gyration values. BBL protein at 330 K in water and GndCl solution showed multiple minima with varying helical content, and correspondingly varying rmsd and gyration values. For E3BD at 330 K in water only medium, we observed two minima: one with retained helical content, and the other with an unfolded helix. No unfolding is observed during rest of the simulations of E3BD, as is also suggested by the representative structures. For human homolog hbSBD, minima are observed at lower rmsd and gyration values, at 300 K and 330 K, in both water only and GndCl solution. The representative structures corresponding to these minima are observed to retain helical content. At 400 K, we observed multiple minima for hbSBD in water only simulations with varying rmsd and gyration values. Representative structures for these minima showed conformations with varying helical content.

12

ACS Paragon Plus Environment

Page 13 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 4 Dynamic cross correlation maps calculated as time-averaged for Cα atoms for. Secondary structure elements: alpha helices, H1 and H2 are annotated with boxes of red color and green color. Interactions of amino acids during simulations Visualization of experimental structures of these subunit binding domains showed that their short parallel  helices are stabilized by the hydrophobic interactions which form a close-packed hydrophobic core (Fig. 1). In E. coli BBL protein, the two helices H1 (residues 14-23) and H2 (residues 40-48) are stabilized by I16, L19, L20 in H1 and V44, L48 in H2. These are further stabilized by hydrophobic interactions with L25, A27, I30, and L39 from the loop region. For E3BD domain, the parallel helices H1 (residues 133-142) and H2 (residues 160-167) enclose a close-packed hydrophobic core consisting of interactions between residues V135, A139 from helix H1; residues I163, L167 and F166 from helix H2, and residues A139, V144, I146, V149, and V156 from long loop connecting the two helices. Another interaction that keeps the two helices in E3BD intact is Y138 of helix H1 forming hydrogen bonding interaction with D164 of helix H2. Interestingly, this interaction 13

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 25

remains intact in all but one simulation (Fig. S4). The loss of interaction within that simulation is consistent with the increase in RMSD and radius of gyration values for E3BD (Fig. 2). In human SBD, it is observed that the helices H1 and H2 are stabilized by hydrophobic interactions between residues A13, V14, L17, A18, N21 from helix H1; K39, I42, L44, L46 from helix H2, and residues I23, L25, V28, I37, and L38 from the long loop connecting these two helices.

Figure

5

Two-dimensional free energy profile (FEP) in (kcal/mol) plotted as function of rmsd and

radius of gyration for proteins at different temperatures and in presence of different solvents. The representative structures are shown for each minima observed by red-colored islands. Helix H1 is shown in red color, and helix H2 in green color. We analyzed the hydrogen bonding interactions between residues of proteins and solvent molecules (Fig. S5). In both water only and 5 M GndCl simulations, we observed decrease in hydrogen bonding interactions within proteins as well as between protein and water molecules as the simulation temperature is increased. The loss of hydrogen bonds in solvent medium is suggestive of exposure of hydrophobic core to the solvent, and thus towards the unfolding of tertiary structure. To further assess the impact of temperature and solvent on the contacts between residues within the subunit binding domains, we calculated two type of contacts for conformations sampled during each simulation run: native contacts and non-native contacts. Whenever a heavy atom of one residue comes within 6Å of a heavy atom of another residue (offset by 3 in sequence), we define it as a contact. We wanted to understand the evolution of the native contacts present in the experimental structure into native and non-native contacts, and to further evaluate the structural differences in proteins with respect to the experimental crystal structure, during the course of the simulations. For this, we calculated the fraction of native contacts (Q) for each conformation sampled during all the simulation runs. Plotting the 14

ACS Paragon Plus Environment

Page 15 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

probability distribution functions for Q (Fig. 6) for these proteins at different temperatures in different solvent environments, we observed loss of native contacts at higher temperature. In water only simulations, there is significant loss of native contacts with rise in temperature. For BBL protein, there is considerable loss of native contacts at 330 K and 400 K in both water only and GndCl simulations. E3BD is thermostable with significant loss of native contacts observed only at 330 K, and ~45% are retained even at 400 K in water only simulations. This can be attributed to the unfolding of protein observed during one of the two longer simulation runs at 330 K. For human hbSBD, the fraction of native contacts observed during simulations is inversely related to the simulation temperature. Higher the temperature, greater is the loss of native contacts in water only simulations. However, in GndCl solutions, native contacts are minimally affected by temperature variations, for both E3BD and hbSBD. We computed the free energy landscape with rmsd and native contacts as reaction coordinates (Fig. S6). In water only simulations, as the simulation temperature increases, we observe that the rmsd values increase and the fraction of native contacts decrease, suggesting unfolding of the proteins. For instance, BBL conformations at 330 K and 400 K; E3BD at 330 K; and hbSBD at 400 K. For GndCl simulations, except for BBL simulation at 330K and 400K, the crowding effect of GndCl seems to keep the protein conformations intact, as is indicated by lower rmsd values and higher fraction of native contacts. Comparative analysis of loss of native contacts with formation of new non-native contacts (Fig. S7) showed that as the fraction of native contacts decreases with increase in temperature, an increase in non-native contacts is observed for BBL, in both water only and GndCl solution, as simulation temperature increases. For E3BD and hbSBD also, there is increase in non-native contacts with increase in temperature for water only simulations. The effect of temperature on formation of nonnative contacts in E3BD is less noticeable in presence of denaturant. Interestingly, we did observe an increase in new non-native contacts formed as temperature is increased for hbBSD protein in presence of denaturant.

15

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 25

Figure 6. Probability distribution functions for fraction of native contacts, Q observed during simulation at different temperatures and solvents.

Intra-molecular interaction dynamics represented as protein networks We represented the intra-molecular interactions of each protein as network with each node representing C atom of amino acid and their interactions as edges. We then weighted the edges based on the strengths of correlated motion and clustered the full network into communities of highly intraconnected but weakly inter-connected nodes using Girvan-Newman clustering method. Normal mode analysis of the experimental structures (Figs. 7, S8), revealed that all three subunit binding domains can be largely represented by three highly intra-correlated structural regions: two for helical residues (one for each helix), and one for the loop connecting the two helices. As the temperature increases, the inter-residue correlation patterns change, and it is observed that for simulations where unfolded conformations are present (For instance, BBL conformations at 330 K and 400 K, E3BD conformations at 330 K and hbSBD at 400 K in water only simulations., residues from either of the helices form communities with residues from another helix or with the residues from the interconnecting loop. In presence of GndCl solutions, changes observed in community patterns are less perceptible for all the three proteins at different temperatures, especially for E3BD, suggesting its tolerance for GndCl as denaturant. In addition, we assessed the distribution of the edges in the protein-networks by measuring the centrality values for each node (Fig. S9). The node centrality was measured by its betweenness value, that is, the number of unique-shortest paths crossing that node. As a general trend, it was observed that for all three subunit domains, the centrality values of helical residues decreased with rise in temperature, indicating loss of network connections or decreased level of coupling upon thermal unfolding.

16

ACS Paragon Plus Environment

Page 17 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 7. Protein Networks with backbone colored as different communities. Native corresponds to experimental structure. Centroids representative of communities is represented as spheres colored according to community they belong to, and the size of sphere is proportional to the number of residues constituting the communities. Correspondingly, full all-residue network is shown in Figure S6. Interactions of subunit binding domains with solvent We analyzed the distribution of solvent molecules around the heavy atoms of proteins within the first shell and second shell. We selected the two layers of solvent molecules in contact with the nonhydrogen atoms of proteins within a distance of 3.5 Å for first shell, and 5 Å for second shell (Fig. S10). We observed that similar number of solvent molecules are present around the shells irrespective of temperature, as inferred from the similar density peaks for distribution of solvent molecules. For simulation in GndCl solution, we observed that the guanidinium molecules promote the exclusion of water molecules from both first and second solvation shell. This distribution of guanidinium and water molecules around the protein also remains same, irrespective of the simulation temperature. We then analyzed the probability of finding solvent atoms in the neighborhood of the solute or protein. For this we selected the backbone oxygen and nitrogen atoms; and looked at the distribution of the distances of the hydrogen and oxygen atoms of solvent molecules from solute (protein) atoms, respectively (Fig. 8A). For all the proteins, at short distances, the radial distribution function g(r) is zero, and then it rises with first solvation shell at 2.5-3.0 Å, and second shell at 5Å. As a general trend, it is observed that the interaction pattern of water molecules with the protein backbone is not significantly affected by the presence of denaturant. However, the backbone oxygen atoms formed more favorable interactions with hydrogen atoms of guanidinium as observed from higher rdf peaks for GndCl solutions (comparing panel 2 with panel 6 of Fig. 8A). We decomposed the rdf calculations to charged, polar, and hydrophobic residues (Fig. S11). For charged residues, we considered aspartate, glutamate, 17

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 25

lysine and arginine atoms. For polar residues, we considered cysteine, methionine, tyrosine, serine and threonine atoms, and for hydrophobic residues, we considered isoleucine, leucine, valine, phenylalanine and tryptophan atoms. The most favorable interactions are observed for negatively charged residues: aspartate and glutamate interacting with hydrogens of guanidinium moieties. We further analyzed the solvent accessibilities of these proteins at different temperatures (Fig. 8B-C). As a general trend, we observed that solvent accessibilities of all three proteins increase considerably at higher temperatures in water only simulations. As the temperature increases, solvent accessibilities of their hydrophilic residues show slight decrease in both water only and GndCl solutions (Fig. 8B). On the other hand, solvent accessibilities of hydrophobic residues increase with temperature for water only simulations. For BBL and hbSBD proteins, we observe broadening of sasa peaks at 330 K and 400 K with respect to peaks observed at 300 K. For E3BD, however, the broadening is not that prominent, and this seems to be consistent with the fact that we observed unfolding of BBL and hbSBD at these temperatures, with E3BD showing thermostability. In presence of GndCl, solvent accessibilities of these proteins remain similar at all temperatures, suggesting crowding effect of guanidinium moieties around the protein. Based on the variation of rmsd values with sasa values, we constructed two-dimensional free energy landscapes for both water only and GndCl simulations (Fig. S12). As a general trend, we observed that in water only simulations the increase in hydrophobic sasa values for BBL and hbSBD correlates with increase in rmsd values, suggesting unfolding of the proteins. For GndCl simulations, the crowding effect of GndCl seems to keep the helical conformation intact for these proteins, and thus, we do not see commensurate increase in sasa values corresponding to the increase in rmsd values.

18

Figure 8 A. Radial distribution plots (rdf) between protein backbone atoms and solvent atoms as computed from simulation data from different simulation runs at different temperatures and in different solvent environment. Distribution peaks depicting solvent accessibilities values of B. hydrophilic and C.ACS hydrophobic SASA values are in Å2. Paragonresidues. Plus Environment

Page 19 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

DISCUSSION In our study, we investigated how these proteins with similar structures and similar functions unfold at different temperatures. We compared the simulations carried out in TIP3P water only model with those carried out in an ionic denaturant solution of 5 M GndCl. From the combined structural deviations (rmsd, radius of gyration) and secondary structure analysis, we observed that effect of elevated temperature is more pronounced in water only simulations compared to that in GndCl solution. PCA showed that these simulations sample different conformations along the first two principal components at different temperatures in presence of water. The time evolution plots showed partial unfolding of BBL at 300 K, and complete unfolding at 330 K and 400 K. For E3BD, partial unfolding is observed at 330 K and 400 K, with complete melting of helix H1 during one of the simulations run at 330 K. For hbSBD, partial unfolding is observed at 330 K, with complete unfolding observed at 400 K. The variation in secondary structure indicated that unfolding of these subunit domains commences with unfolding of one of the helices, preferably Cterminal. Comparative analysis of chemical denaturation studies at different temperatures suggested that crowding of Gnd+ ions around the protein may be screening the protein from temperature effect. Our radial distribution (rdf) calculations for protein atoms and Gnd+ ions suggested that more than the peptide backbone, the side chains of negatively charged residues formed favorable interactions with Gnd+. Presence of Gnd+ within the first and second hydration shell did not modify the rdf of water in the first two solvent shells appreciably. This effect appeared to be in agreement with the results of hydrogen exchange NMR experiments, which showed that GndCl cannot form hydrogen bonds with peptide39. Another 2D-IR study carried out to investigate changes in secondary structure of two proteins showed that Gnd+ disrupts -sheets much more effectively than -helices, and a denaturation mechanism was suggested involving a strong decrease of the hydrophobic effect via direct interactions of Gnd+ with the side chains of the residues, thus promoting the exposure of hydrophobic core residues40. Similarities and differences in the unfolding behavior of the three subunit binding domains The PSBDs studied shared similar structures with two parallel helices connected by a long interconnecting loop. The two helices are maintained by hydrophobic core including residues from both helices and the interconnecting loop. The overall native structures of these proteins are composed of networks between three major communities of residues, one for each helix and the third for the interconnecting loop. Based on the correlations observed during elevated temperatures and in the presence of different solvents, residues in these proteins form new communities. Most of the residues 19

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 25

get correlated in completely unfolded states, and they give rise to fewer, albeit larger, communities. For all three PSBDs, the computed free energy landscapes corresponding to parameters such as radius of gyration, fraction of native contacts, and hydrophobic sasa versus common reaction coordinate of rmsd show consistent patterns. With increase in rmsd values the structure unfolds, increasing the radius of gyration values; and within these partially or fully unfolded structures, the fraction of native contacts decreases and solvent accessibilities of hydrophobic residues increase. -value analyses and MD simulations of pseudo wild-type BBL-H142W41, QNND-BBL13 and replica exchange molecular dynamics (REMD) simulations13 of BBL constructs showed higher helical propensity of BBL (than E3BD, or its hyperthermophilic homolog POW), to form secondary structure near its N-terminal helix. In addition, REMD simulations of BBL constructs showed melting temperature of helix H1 to be higher than that of helix H2 and residual helical structure for helix H1 could be observed even when the helix H2 is completely lost. Our BBL simulations at 300 K showed partial melting of C-terminal helix H2 at 300 K with almost 50% loss of native contacts. At higher temperatures of 330 K and 400 K, BBL showed cooperative unfolding transitions, where unfolding of BBL initiates via melting of C-terminal helix H2. This resulted in the loss of inter-helical interactions or native contacts, which is followed by melting of N-terminal helix H1. Our results are in agreement with the reported REMD simulations13 of BBL constructs. E3BD is observed to be the most stable of all, secondary structures being retained during all (but one) simulations. Our results are consistent with the -value analyses carried out on pseudo wildtype variant of E3BD (E3BD-F166W)20, which showed that a similar pattern of F values was observed at 298K and 325K, and at high concentrations of GndCl42. As far as the sequence is concerned, E3BD consists of an additional tyrosine (Y138) in helix H1, in place of leucine, as in both hbSBD and BBL. Y138 forms hydrogen bonding interaction with D164 of helix H2. The stability of E3BD, with respect to elevation of temperatures and presence of Gnd+, can be attributed to this hydrogen bonding interaction. Arguably, this interaction keeps both the helices intact, thus allowing the protein to function at higher temperatures for B. stearothermophilus. On mutating these residues, either single mutation (Y138F and D164N) or double mutation (Y138F-D164N) and comparing free energy of interaction in both native and folding transition states, it was observed that this interaction contributed to the stability of the transition state for folding of E3BD. Control simulations of E3BD unfolding at 400 K by Pitera et al13 (similar to simulations of Ferguson et al.20 of E3BD-F166W) showed E3BD to be largely folded with extensive native-like secondary structure. Our E3BD simulations are longer than in the previous simulations, and we observed that at 400 K, though partial melting of helices is observed, we do not observe unfolding of E3BD within 100ns. Nonetheless, we do observe unfolding of E3BD during longer simulations at 330 K where unfolding begins with partial

20

ACS Paragon Plus Environment

Page 21 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

melting of the C-terminal helix H2 at around ~380ns, finally resulting in complete melting of Nterminal helix H1 at around ~450ns. The partial helical structure of H2 is however still retained. For human homolog hbSBD, unfolding is observed at elevated temperatures of 400 K within only 100ns. There is partial melting of C-terminal helix H2 observed during longer simulation runs at 330 K. Based on mapping of the melting of helices with the observed transitions in secondary structure elements at 330 K and 400 K, it is safe to suggest that for hbSBD, unfolding commences via melting of the C-terminal helix, which may lead to cooperative unfolding of the N-terminal helix H1 resulting in complete unfolding. Fluorescence equilibrium and kinetic measurements8 suggest that hbSBD folds via barrier-limited transition. The two-state folding mechanism of hbSBD was also captured by Langevin dynamics of coarse-grained Go model7. Our comparative studies suggest that BBL protein is the most sensitive of all the three proteins, with hbSBD of intermediate stability, and E3BD most stable of all three. For both BBL and hbSBD, Cterminal helix H2 is less stable than N-terminal helix H1, whereas for E3BD, though the partial melting commences with C-terminal helix H2, the complete unfolding of N-terminal helix H1 precedes that of helix H2. This was also observed by Fersht group20,41 where the pattern of Φ values determined for BBL was different from those of its homologues, E3BD and POB. While folding of E3BD and POB was observed to nucleate in helix H220, folding nucleus of BBL was less biased towards the C-terminal interaction networks and showed a greater energetic role of the N-terminal ones41. Although no Φ values based evidence is available for unfolding of the human homolog hbSBD, our results appear to suggest that for hbSBD, folding should nucleate in N-terminal helix H1, recruiting associated interaction networks later in the folding transition state, and thus folding occurs via barrier limited two-state transition. Our simulations show one of the possible folding pathways of three PSBDs, where unfolding commences via unfolding of C-terminal helix, leading to the exposure of hydrophobic core to solvent molecules, loss of native contacts and to the eventual disruption of the tertiary structure. The higher stability of N-terminal helix H1 in comparison to that of C-terminal helix H2 is compatible with the fact that it is the N-terminal helix H1 of the subunit binding domain which is reported to be interacting with either the E14 or E35 component of assembled multi-enzyme complex. BBL belongs to mesophilic organism (E. coli) with low TM value and our simulations also suggest BBL is more sensitive to thermal and chemical perturbation. While for human homolog hbSBD, we observe intermediate stability, E3BD seems to be the most stable of all three PSBDs. This is also expected since it is required by the thermophilic organism (B. stearothermophilus) for its survival. CONCLUSION

21

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 25

We have carried out simulations of three peripheral subunit binding domains: BBL, E3BD and hbSBD at different temperatures, and in water and 5 M GndCl mediums respectively. Our studies showed that thermal denaturation of these domains commences via unfolding of the C-terminal helix. While doing so, the native interactions are replaced by new non-native interactions, exposing the hydrophobic residues to the solvent, which in turn results in the unfolding of the protein. The thermal denaturation effect is more pronounced in simulations with water as solvent. High concentration of guanidinium molecules seems to produce a screening effect which keeps the helical structure intact. Here, the Gnd+ ions interact favorably with the side chains of aminoacids, notably negatively charged residues. It is observed that E3BD is the most stable of all while BBL is most sensitive to thermal and chemical denaturation. hbSBD is of intermediate stability. The thermostability of E3BD seems to be attributed to its hydrogen bonding interaction of Y138 of helix H1 with D164 of helix H2. Our results suggest that folding of PSBDs show cooperative unfolding transitions where unfolding of domain begins via melting of C-terminal helix H2, resulting in loss of interhelical interactions or native contacts, and followed by melting of N-terminal helix H1. Acknowledgements MS thanks Department of Science and Technology (DST), India for INSPIRE Award and research grant (IFA14-CH-165). Supplementary Information Supplementary document consists of figures S1-S12. Reference: (1)

(2)

(3)

(4)

22

Robien, M. A.; Clore, G. M.; Omichinski, J. G.; Perham, R. N.; Appella, E.; Sakaguchi, K.; Gronenborn, A. M. Three-Dimensional Solution Structure of the E3-Binding Domain of the Dihydrolipoamide Succinyltransferase Core from the 2-Oxoglutarate Dehydrogenase Multienzyme Complex of Escherichia Coli. Biochemistry 1992, 31, 3463–3471. Allen, M. D.; Broadhurst, R. W.; Solomon, R. G.; Perham, R. N. Interaction of the E2 and E3 Components of the Pyruvate Dehydrogenase Multienzyme Complex of Bacillus Stearothermophilus. Use of a Truncated Protein Domain in NMR Spectroscopy. FEBS J. 2005, 272, 259–268. Chang, C.-F.; Chou, H.-T.; Lin, Y.-J.; Lee, S.-J.; Chuang, J. L.; Chuang, D. T.; Huang, T. Structure of the Subunit Binding Domain and Dynamics of the Di-Domain Region from the Core of Human Branched Chain α-Ketoacid Dehydrogenase Complex. J. Biol. Chem. 2006, 281, 28345–28353. Mande, S. S.; Sarfaty, S.; Allen, M. D.; Perham, R. N.; Hol, W. G. Protein–Protein Interactions in the Pyruvate Dehydrogenase Multienzyme Complex: Dihydrolipoamide Dehydrogenase Complexed with the Binding Domain of Dihydrolipoamide Acetyltransferase. Structure 1996, 4, 277–286. ACS Paragon Plus Environment

Page 23 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15)

(16) (17) (18) (19) (20) (21) (22) (23) (24) (25)

23

The Journal of Physical Chemistry

Frank, R. A. W.; Pratap, J. V.; Pei, X. Y.; Perham, R. N.; Luisi, B. F. The Molecular Origins of Specificity in the Assembly of a Multienzyme Complex. Structure 2005, 13, 1119–1130. Spector, S.; Kuhlman, B.; Fairman, R.; Wong, E.; Boice, J. A.; Raleigh, D. P. Cooperative Folding of a Protein Mini Domain: the Peripheral Subunit-Binding Domain of the Pyruvate Dehydrogenase Multienzyme Complex J. Mol. Biol. 1998, 276, 479–489. Kouza, M.; Chang, C.-F.; Hayryan, S.; Yu, T.; Li, M. S.; Huang, T.; Hu, C.-K. Folding of the Protein Domain HbSBD. Biophys. J. 2005, 89, 3353–3361. Arbely, E.; Neuweiler, H.; Sharpe, T. D.; Johnson, C. M.; Fersht, A. R. The Human Peripheral Subunit-Binding Domain Folds Rapidly while Overcoming Repulsive Coulomb Forces. Protein Sci. 2010, 19, 1704–1713. Fan, J.; Duan, M.; Li, D.-W.; Wu, H.; Yang, H.; Han, L.; Huo, S. Observation of Two Families of Folding Pathways of BBL. Biophys. J. 2011, 100, 2457–2465. Li, P.; Oliva, F. Y.; Naganathan, A. N.; Muñoz, V. Dynamics of One-State Downhill Protein Folding. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 103–108. Lin, C.-W.; Culik, R. M.; Gai, F. Using VIPT-Jump to Distinguish Between Different Folding Mechanisms: Application to BBL and a Trpzip. J. Am. Chem. Soc. 2013, 135, 7668–7673. Yu, W.; Chung, K.; Cheon, M.; Heo, M.; Han, K.-H.; Ham, S.; Chang, I. Cooperative Folding Kinetics of BBL Protein and Peripheral Subunit-Binding Domain Homologues. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 2397–2402. Pitera, J. W.; Swope, W. C.; Abraham, F. F. Observation of Noncooperative Folding Thermodynamics in Simulations of 1BBL. Biophys. J. 2008, 94, 4837–4846. Ferguson, N.; Sharpe, T. D.; Schartau, P. J.; Sato, S.; Allen, M. D.; Johnson, C. M.; Rutherford, T. J.; Fersht, A. R. Ultra-Fast Barrier-Limited Folding in the Peripheral SubunitBinding Domain Family. J. Mol. Biol. 2005, 353, 427–446.. Liu, Y.-J.; Soumelis, V.; Watanabe, N.; Ito, T.; Wang, Y.-H.; Malefyt, R. de W.; Omori, M.; Zhou, B.; Ziegler, S. F. TSLP: An Epithelial Cell Cytokine That Regulates T Cell Differentiation by Conditioning Dendritic Cell Maturation. Annu. Rev. Immunol. 2007, 25, 193–219. Knott, M.; Chan, H. S. Criteria for Downhill Protein Folding: Calorimetry, Chevron Plot, Kinetic Relaxation, and Single-Molecule Radius of Gyration in Chain Models with Subdued Degrees of Cooperativity. Proteins 2006, 65, 373–391. Zuo, G.; Wang, J.; Wang, W. Folding with Downhill Behavior and Low Cooperativity of Proteins. Proteinss 2006, 63, 165–173. Garcia-Mira, M. M.; Sadqi, M.; Fischer, N.; Sanchez-Ruiz, J. M.; Muñoz, V. Experimental Identification of Downhill Protein Folding. Science 2002, 298, 2191–2195. Ferguson, N.; Schartau, P. J.; Sharpe, T. D.; Sato, S.; Fersht, A. R. One-State Downhill versus Conventional Protein Folding. J. Mol. Biol. 2004, 344, 295–301. Ferguson, N.; Day, R.; Johnson, C. M.; Allen, M. D.; Daggett, V.; Fersht, A. R. Simulation and Experiment at High Temperatures: Ultrafast Folding of a Thermophilic Protein by Nucleation-Condensation. J. Mol. Biol. 2005, 347, 855–870. Zhang, J.; Li, W.; Wang, J.; Qin, M.; Wang, W. All-atom Replica Exchange Molecular Simulation of Protein BBL. Proteins. 2008, 72, 1038–1047. Prakash, A.; Kumar, V.; Meena, N. K.; Lynn, A. M. Elucidation of the Structural Stability and Dynamics of Heterogeneous Intermediate Ensembles in Unfolding Pathway of the NTerminal Domain of TDP-43. RSC Adv. 2018, 8, 19835–19845. Kumar, V.; Prakash, A.; Lynn, A. M. Alterations in Local Stability and Dynamics of A4V SOD1 in the Presence of Trifluoroethanol. Biopolymers 2018, 109, e23102. Kumar, V.; Prakash, A.; Pandey, P.; Lynn, A. M.; Hassan, M. I. TFE-Induced Local Unfolding and Fibrillation of SOD1: Bridging the Experiment and Simulation Studies. Biochem. J. 2018, 475, 1701–1719. Prakash, A.; Kumar, V.; Pandey, P.; Bharti, D. R.; Vishwakarma, P.; Singh, R.; Hassan, M. I.; Lynn, A. M. Solvent Sensitivity of Protein Aggregation in Cu, Zn Superoxide Dismutase: A Molecular Dynamics Simulation Study. J. Biomol. Struct. Dyn. 2018, 36, 2605–2617. ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 25

(26) Prakash, A.; Kumar, V.; Meena, N. K.; Hassan, M. I.; Lynn, A. M. Comparative Analysis of Thermal Unfolding Simulations of RNA Recognition Motifs (RRMs) of TAR DNA-Binding Protein 43 (TDP-43). J. Biomol. Struct. Dyn. 2019, 37, 178–194. (27) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1–2, 19–25.. (28) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials for the Amber Ff99SB Protein Force Field. Proteins 2010, 78, 1950–1958. (29) Weerasinghe, S.; Smith, P. E. A Kirkwood-Buff Derived Force Field for the Simulation of Aqueous Guanidinium Chloride Solutions. J. Chem. Phys. 2004, 121, 2180–2186. (30) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33–38. (31) Kabsch, W.; Sander, C. Dictionary of Protein Secondary Structure: Pattern Recognition of Hydrogen-Bonded and Geometrical Features. Biopolymers 1983, 22 (12), 2577–2637. (32) Grant, B. J.; Rodrigues, A. P. C.; ElSawy, K. M.; McCammon, J. A.; Caves, L. S. D. Bio3d: An R Package for the Comparative Analysis of Protein Structures. Bioinformatics 2006, 22, 2695–2696. (33) Roe, D. R.; Cheatham, T. E. PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theory Comput. 2013, 9 (7), 3084–3095. (34) Harte, W. E.; Swaminathan, S.; Mansuri, M. M.; Martin, J. C.; Rosenberg, I. E.; Beveridge, D. L. Domain Communication in the Dynamical Structure of Human Immunodeficiency Virus 1 Protease. Proc. Natl. Acad. Sci. U.S.A. 1990, 87, 8864–8868. (35) Amadei, A.; Linssen, A. B.; de Groot, B. L.; van Aalten, D. M.; Berendsen, H. J. An Efficient Method for Sampling the Essential Subspace of Proteins. J. Biomol. Struct. Dyn. 1996, 13, 615–625. (36) Amadei, A.; Linssen, A. B. M.; Berendsen, H. J. C. Essential Dynamics of Proteins. Proteins 2004, 17, 412–425. (37) Girvan, M.; Newman, M. E. J. Community Structure in Social and Biological Networks. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 7821–7826. (38) Yao, X.-Q.; Malik, R. U.; Griggs, N. W.; Skjærven, L.; Traynor, J. R.; Sivaramakrishnan, S.; Grant, B. J. Dynamic Coupling and Allosteric Networks in the α Subunit of Heterotrimeric G Proteins. J. Biol. Chem. 2016, 291, 4742–4753. (39) Lim, W. K.; Rösgen, J.; Englander, S. W. Urea, but Not Guanidinium, Destabilizes Proteins by Forming Hydrogen Bonds to the Peptide Group. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 2595–2600. (40) Huerta-Viga, A.; Woutersen, S. Protein Denaturation with Guanidinium: A 2D-IR Study. J. Phys. Chem. Lett. 2013, 4, 3397–3401. (41) Neuweiler, H.; Sharpe, T. D.; Rutherford, T. J.; Johnson, C. M.; Allen, M. D.; Ferguson, N.; Fersht, A. R. The Folding Mechanism of BBL: Plasticity of Transition-State Structure Observed within an Ultrafast Folding Protein Family. J. Mol. Biol. 2009, 390, 1060–1073. (42) Ferguson, N.; Sharpe, T. D.; Johnson, C. M.; Fersht, A. R. The Transition State for Folding of a Peripheral Subunit-Binding Domain Contains Robust and Ionic-Strength Dependent Characteristics. J. Mol. Biol. 2006, 356, 1237–1247.

24

ACS Paragon Plus Environment

Page 25 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

TOC Graphic

25

ACS Paragon Plus Environment