Phosphorylation versus O-GlcNAcylation: Computational Insights into

Subscriber access provided by READING UNIV

Article

Phosphorylation versus O-GlcNAcylation: Computational Insights into the Differential Influences of the Two Competitive Post-Translational Modifications Lata Rani, and Sairam Swaroop Mallajosyula J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b08790 • Publication Date (Web): 27 Oct 2017 Downloaded from http://pubs.acs.org on October 31, 2017

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry B 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 53

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

Phosphorylation versus O-GlcNAcylation: Computational Insights into the Differential Influences of the Two Competitive PostTranslational Modifications Lata Rani and Sairam S. Mallajosyula* Department of Chemistry, Indian Institute of Technology Gandhinagar, Palaj, Gandhinagar, Gujarat, India - 382355 E-mail: [email protected] Phone: +91 - 79 – 2395 2459 Abstract Phosphorylation and O-GlcNAcylation are rapidly cycling intracellular protein post-translational modifications (PTMs) that can compete for the same serine (S) and threonine (T) sites. Limited crystal structure information is available on the direct influence of these PTMs on the underlying protein structure, especially for O-GlcNAcylation. NMR and CD studies show that these competitive-PTMs can have same or differential influence on the overall secondary structure. In Tau derived peptide fragments it was found that phosphorylation stabilized PPII conformations while O-GlcNAcylation destabilized the same. In the absence of substantial structural information we have performed a systematic computational study utilizing PDB analysis, QM calculations and MD simulations to identify key structural trends upon PTM. Our analysis of the limited PDB dataset revealed conformational shifts from PPII to α-helical geometry upon serine phosphorylation and in the opposite direction, from α-helical to PPII geometry upon threonine phosphorylation. Gas phase QM calculations covering the complete Ramachandran φ/ψ space 1 ACS Paragon Plus Environment


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

using model native, phosphorylated and O-GlcNAcylated dipeptide systems revealed preferences towards α-helical conformations. However the major structural transitions were observed in the MD simulations upon the inclusion of solvation. The model dipeptide simulations revealed a preference for PPII and α-helical conformations for phosphorylated serine and threonine, while O-GlcNAcylated dipeptides exhibited a complete shift towards extended conformations, β-sheet and PPII, disfavoring the α-helical conformation. For the Baldwin α-helix simulations it was found that both phosphorylation and O-GlcNAcylation destabilized the helix, however the destabilization was governed by H-bonding and electrostatic interactions in the former while the later was controlled by hydrophobic collapse and steric interactions. The presence of lysine in close proximity of phosphate leads to potentially stable salt-bridge interactions, which can influence the structure based on the relative placement of the lysine with respect to the PTM site. Similar strong lysine-phosphate contacts were observed in the model Tau peptides, which steers the conformations towards PPII geometries, highlighting the direct influence of the PTM on function.

2 ACS Paragon Plus Environment

Page 2 of 53

Page 3 of 53

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


Introduction Phosphorylation and O-GlcNAcylation are well-established intracellular post-translational modifications (PTMs) of Serine (Ser) and Threonine (Thr) that modify the protein at the structural and functional level. 1-2 The cycling of the phosphate group in humans is controlled by a large collection of protein kinases and phosphatases that phosphorylate and dephosphorylate the protein at specific Ser, Thr and Tyrosine (Tyr) residues. 3 In contrast, both the addition and the removal of β-O-GlcNAc to the side chains of Ser/Thr is controlled by a single β-O-GlcNAc transferase (OGT, three isoforms in humans) and a single O-GlcNAcase (OGA). Despite the stark difference in the number of enzymes involved in the cycling of the two PTMs, similar times scales have been observed for the two post-translational events.1 Since both the PTMs target the same intracellular site (Ser/Thr) significant crosstalk has been observed between the two PTMs.

1

Seminal work by Hart and co-workers have revealed that O-GlcNAcylation can

both elevate and lower phosphate stoichiometry indicating a complex interplay between the two pathways.

1, 4-5

The importance of understanding the interplay between the two PTMs becomes

crucial, as the crosstalk between the PTMs has been implicated in signaling, transcription and diseases.6-10

Even though extensive research is being carried out to understand the biochemical influence of these PTMs, the structural basis of understanding the influence of phosphorylation versus OGlcNAcylation has not been addressed in detail. It is reported that close to 30% of proteins are phosphorylated,11 while close to 4000 proteins have been identified as being O-GlcNAcylated with the numbers increasing as more sophisticated profiling and characterization tools being developed to study O-GlcNAcylation.12 However crystallographic information is available only for a very small percentage of these structures with only 285 non-redundant atomic resolution crystal structures being reported for phosphorylated proteins in the Protein Data Bank (PDB).13 This number drops drastically for the O-GlcNAcylated proteins with only 15 non-redundant structures being reported in the PDB. In fact, there are no direct X-ray crystallographic studies on the site-specific competitive influence of both the PTMs. NMR spectroscopy is the major tool used to study the structural influence of these PTMs. 14-17



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

Between the two PTMs the influence of phosphorylation on various protein secondary structure conformations, like model α-helical peptides, PPII conformations and β-hairpin structures have been studied extensively using NMR spectroscopy.

18-21

Systematic studies on position specific

phosphorylation and their influence on model α-helices have been conducted by the groups of Doig18 and Zondlo22. In their studies it was found that N-terminal phosphorylation generally stabilized the α-helix while phosphorylation in the interior of the helix generally destabilized the helix. However Doig found that phosphoserine in the interior of a helix could stabilize the helix if it was located in a i, i+4 position with respect to lysine.23 Phosphorylation and its influence on PPII and coil conformations has been specially relevant due to its implications on microtubule binding protein Tau. Tau is an intrinsically disordered neuronal protein that is directly involved in pathogenesis of Alzheimer and other neurodegenerative diseases.24 It has been observed that Tau is hyperphosphorylated when it aggregates as fibrils that are observed in Alzheimer’s disease. The longest isoform of Tau (441 residues) contains 85 potential phosphorylation sites. Of these phosphorylation atleast three sites S214, T231 and S262 have been directly implicated Alzheimer’s.

25-27

Structural shifts upon phosphorylation have been measured by the downfield

and up-field shifts of amide protons sometimes even distant to the phosphorylation site. In 2014 Zondlo and co-workers reported two detailed NMR studies comparing the influence of phosphorylation and O-GlcNAcylation on model α-helices and short peptide stretches selected from the Tau protein.22, 28 It was found that phosphorylation and O-GlcNAcylation have similar structural effects in α-helices while they have opposing effect on the PPII geometry observed in the Tau fragments. The positioning of the PTMs in the helix had a direct influence on the structural stability and observed effects. It was found that both the PTMs stabilized the helical conformation when placed at the N-termini of the helix while they completely destabilized the helical conformation when placed in the interior of the helix. In the Tau peptide fragments it was found that phosphorylation favored the PPII conformations while O-GlcNAcylation destabilized the same. This observation was reported to be significant, as it had been reported that increasing O-GlcNAc stabilizes tau against aggregation as opposed to hyperphosphorylation that enhances its aggregation.29 Computational studies can provide valuable role into elucidating the mechanisms of PTM at the atomistic level that can aid in bridging the structure-function gap. In fact phosphorylation is one 4 ACS Paragon Plus Environment

Page 4 of 53

Page 5 of 53

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


of the most well studied PTM computationally. Computational studies, especially molecular dynamics (MD) simulations have been used extensively to study phosphorylation and its influence (conformational changes) on the protein structures. These studies have been reviewed earlier by Johnson et. al and Narayanan et. al.13, 30 The influence of O-GlcNAcylation on protein structure has received lesser attention owing to the non-availability of X-ray crystal structures. Fernandez-Tajeda and co-workers used time-averaged restraints MD simulations (Tar-MD) simulations, with dihedral and NOE constraints derived from NMR measurements, to gain insight into the geometrical features of O-GlcNAcylation.14 The use of constraints was circumvented by Mallajosyula and co-workers by employing enhanced sampling protocols to study the influence of O-Glycosylation on model peptides.31 Zerze et. al studied the influence of glycosylation on a group of intrinsically disordered proteins by using temperature replica exchange protocols. 32 However there are no computational studies in the literature that have compared the competitive structural influence of phosphorylation and O-GlcNAcylation. We have directed this study to obtain a better picture of the conformational preferences upon phosphorylation and OGlcNAcylation. The aim of the study is to elucidate the key interactions that are involved in the PTMs, which can further be used to describe the observed functions. The study is divided into three sections. In the first section we perform a detailed structural and conformational analysis of phosphorylated serine and threonine residues reported PDB to reveal structural trends observed in crystal structures. In the second section we perform extensive QM calculations on model phosphorylated and O-GlcNAcylated peptides covering the complete Ramachandran φ/ψ space, to gain insight on the secondary structure preference upon PTM. In the third section we perform extensive MD simulations on post-translationally modified dipeptides, model α-helices and Tau fragments. The Baldwin α-helical peptides and proline rich region of the Tau peptide were based on the NMR and CD studies performed by Zondlo and co-workers, which served as primary target data for comparison and validation of our MD simulations. 22, 28



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

Experimental Section PDB Data Crystal structure data was obtained from the Protein Data Bank database (PDB, www.rcsb.org) to generate two data sets.33 The first dataset corresponded to protein structures containing phosphorylated Serine and Threonine while the second one was corresponding to OGlcNAcylated Serine and Threonine. The first dataset was created by searching for the residue identifiers SEP and TPO in the PDB database, which corresponds to phosphorylated Serine and Threonine respectively. Hereafter phosphorylated serine and threonine are referred to as SEP and TPO in the manuscript. The second dataset was created by searching for the residue identifier NAG corresponding to β-N-acetyl D-glucosamine in PDB. For the identified PDB files the link information in the PDB header was used to identify structures wherein β-N-acetyl Dglucosamine was covalently attached to Serine and Threonine. The analysis revealed 799 SEP sites in 686 structures, 605 TPO sites in 586 structures, 10 O-GlcNAc (S) sites and 5 O-GlcNAc (T) sites in 32 structures. Protein structure databases are prone to redundancy with multiple entries corresponding to a family of structurally similar proteins having a sequence similarity of > 90%. The inclusion of such structurally similar proteins in the analyses introduces undesirable biases particularly when trends have to be extracted from the data. To remove such a structural bias a three-step redundancy removal protocol was applied. In the first step identical sequences were identified using the standalone BLAST tool from NCBI by aligning the FASTA sequence of all proteins present in the dataset. A sequence similarity threshold of 90% was used to cluster the protein structures into one batch. The structure solved at the highest resolution among the batch was then used for further analysis. The resulting structures were then analyzed for repeating subunits. Only one subunit was selected from the proteins with multiple subunits to avoid multiple representations of the PTM sites from the same structure. For NMR resolved structures only the representative structure based on RMSD clustering was included in the dataset. In addition to the removal of redundancy the dataset was also analyzed for missing residues immediately before or after the PTM site. Any PTM site with missing adjacent residues (coordinates) was removed 6 ACS Paragon Plus Environment

Page 6 of 53

Page 7 of 53

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


from the dataset since these coordinates cannot be used to evaluate protein φ (Ci-1-Ni-Cαi-Ci) and ψ (Ni-Cαi-Ci-Ni+1) dihedrals, where i is the identity of the PTM site. The removal of redundancy resulted in 179 sites for SEP and 101 sites for TPO from 179 and 106 PDB entries respectively. For β-O-GlcNAc (S/T) the same protocol resulted in 10 sites for β-O-GlcNAc (S) and 5 sites for β-O-GlcNAc (T) from 10 and 5 PDB structures respectively. The results of the analysis are summarized in Table 1. Table 1. Summary of the PDB analysis. PDB Database

After applying redundancy filter

Sites

Structure

Sites

Structure

SEP

799

686

179

179

TPO

605

586

101

106

O-GlcNAc (S)

10

32

10

10

O-GlcNAc (T)

5

32

5

5

SERa

515

1035

92

92

THRa

303

966

34

34

a

Serine and Threonine sites corresponding to PTM SEP and TPO sites in datasets 1 and 2.

PDB data without PTM In order to analyze the influence of the absence of the PTM at a known PTM site a dataset was created which contained all the structures present in the PDB database sans the structures containing residue identifiers SEP and TPO, ie phosphorylated serine and threonine. This dataset was then curated to identify entries which shared > 90% sequence similarity with sequences present in the SEP and TPO datasets (datasets 1 and 2 in this study). The similarity analysis thus identifies serine and threonine sites that were not post-translationally modified. The resulting dataset contained 515 SER and 303 THR sites in 1035 and 966 PDB entries. The earlier mentioned redundancy filtering protocol was applied to arrive at a dataset independent of redundant structures, which contained 92 SER and 34 THR sites in 92 and 34 PDB entries respectively (Table 1). 7 ACS Paragon Plus Environment


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 53

To perform conformational analysis (φ/ψ analysis) on the curated dataset it is important to make sure that the residue numbers of SER and THR present in dataset align with SEP and TPO sites present in the PTM dataset. While FASTA files contain the complete amino acid sequence, PDB files do not necessarily contain the coordinate information for all amino acid residues. PDB files are numbered by omitting non-crystallized residues and thus can have discontinuous residue numbering. In order to match the exact residue number from phosphorylated and nonphosphorylated datasets, PDB files were renumbered to avoid the chances of extracting incorrect residues. Sequence alignment of phosphorylated SER/THR dataset with non-phosphorylated data gives the residue number of SER/THR aligned with the corresponding SEP/TPO site. The coordinates of these residues were then extracted for further analysis. Similar analysis was not performed for β-O-GlcNAc (S/T) owing to the statistically lower number of PDB structures containing β-O-GlcNAc (S/T). In addition to the curated dataset the PDB database sans the structures containing residue identifiers SEP and TPO was also used to prepare a database of Ser and Thr sites that were never involved in PTM. All the entries which shared > 90% sequence similarity with sequences present in the SEP and TPO datasets (datasets 1 and 2 in this study) were removed from the dataset thus creating a database of serine and threonine sites that were never involved in PTM. The same redundancy filter was applied to the database. The serine and threonine sites thus identified are referred to as Ser’ and Thr’ in the remainder of the study. DFT study Ab initio quantum mechanical (QM) calculations for model dipeptides (Figure 1) were carried out using the Gaussian 03 suite of programs.

34

The minimal dipeptide design accounts for φ (C-

N-Cα-C) and ψ (N-Cα-C-N) dihedrals by adding acetyl (-COCH3) and methylamine (-NHCH3) groups at amino and carboxyl terminal respectively. Such model compound geometries have been used in earlier studies to explore the Ramachandran ɸ/ψ space.

35

Initial geometries for the

QM ɸ/ψ scans were generated using the topology information present in the CHARMM force fields and the CHARMM simulation program.36 Dipeptide conformations were generated at 15° intervals for both ɸ and ψ dihedrals in the ɸ/ψ space amounting to 576 QM calculations per model dipeptide. In accordance with PDB survey data and NMR studies the N-Cα-Cβ-OG (χ) dihedral was fixed at 60° for the O-GlcNAc (S/T) model dipeptides, while all for SEP and TPO 8 ACS Paragon Plus Environment

Page 9 of 53

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


dipeptides the N-Cα-Cβ-OG (χ) and Cα-Cβ-O-P dihedrals were fixed at -60° and 110° respectively.14, 22, 28 Additionally in the SEP and TPO QM scans the N-H bond length was fixed to a value of 1.023 Å, obtained from averaging the N-H bond lengths from QM optimized geometries of the model compounds, to eliminate the possible shift of the acetyl or amine hydrogen’s to the phosphate oxygen’s. The QM calculation were performed at the B3LYP level employing the 6-31+G(d) functional. 37-38

Figure 1. Model Serine/Threonine and modified Serine/Threonine dipeptides used for the QM calculations Molecular Dynamics Simulations Model Dipeptides: Molecular dynamics (MD) simulation for all the dipeptides (Figure 1), were performed using the computational packages CHARMM

36

and GROMACS

. The CHARMM

39

protein, carbohydrate force fields and the TIP3P water model were used for the simulations.

40-42

The initial conformations were generated using the topology information present in the force 9 ACS Paragon Plus Environment


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 53

field. For all the systems three simulations were initiated with the N-Cα-Cβ-OG (χ) dihedral initially fixed at 60°, -60° and 180°. These dihedrals were not restrained during the simulations allowing complete conformational flexibility. A previously reported protocol was used to setup the MD simulations. 31 “ Briefly the initial geometries were subjected to a 500 step steepest descent (SD) minimization followed by 10000 steps of adopted basis Newton-Raphson (ABNR) minimization to a force gradient tolerance of 10-6 kcal/mol/Å.43 The minimized dipeptides were then immersed in a preequilibrated cubic TIP3P water box that extends 14 Å from the dipeptide molecule. Water molecules that overlapped with non-hydrogen atoms of dipeptide in a range of 2.8 Å distance cutoff were deleted. The CRYSTAL module in CHARMM was used to set up the boundary conditions for all subsequent minimizations and MD simulations. A short minimization cycle of 50 SD steps, followed by 50 ABNR steps with a harmonic restraint of 0.1 kcal/mol/Å on the non-hydrogen atoms of dipeptides was used to remove bad contacts. This was followed by a 100 ps simulation in NVT ensemble with the same harmonic restraints to equilibrate the solvent molecules around the dipeptides. A 200 ps NPT simulation at 1 atm and 298 K followed the NVT simulation, wherein all the previous restraints were removed. In the NPT simulation the center of mass of the dipeptide was restrained near the origin by using the MMFP module in CHARMM using a harmonic restraint of 1.0 kcal/mol/Å.

44

The electrostatic interactions were

treated via the particle mesh Ewald method with a real-space cutoff of 12 Å, a kappa value of 0.34 Å-1, and a sixth-order spline.

45

Nonbond interaction lists were updated heuristically out to

16 Å with a force switch smoothing function from 10 to 12 Å used for the Lennard-Jones interactions.

46

The dimensions of the systems are reported in Table S1 of the supporting

information (SI) file. ” The equilibrated structures were used to initiate 100 ns production simulations performed with GROMACS version 5.1.2 simulation package.

39

For the GROMACS simulations protein and

non-protein atoms were coupled to their own temperature baths, which were kept constant at 303.15 K using Nose-Hoover algorithm. the Parrinello-Rahman barostat. algorithm.

49

48

47

Pressure was maintained isotropically at 1 atm using

All hydrogen bond lengths were constrained using the LINCS

A 12 Å cutoff was used for Lennard-Jones interactions. Electrostatic interactions

were calculated using the Particle-Mesh Ewald (PME) method with 12 Å cutoff. 10 ACS Paragon Plus Environment

Page 11 of 53

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


Baldwin helix and Tau peptide simulations: To examine the effects of PTMs on the structure of compact α-helical geometries and extended PPII structures MD simulations were performed for a selection of peptides studied by Zondlo. et. al.

22

In

the

experimental

study

the

14

residue

long

Baldwin

peptide

(Ac-

AKAAAAKAAAAKAA-NH2) was used as a prototype for the α-helical structure. PTM were introduced into the structure at the N-terminus end (position 1), interior (position 5) and near the C-terminus (position 10) end to study their influence on the structure. In the experimental study structural changes were characterized using CD and NMR spectroscopy. Significant structural changes were observed for PTM introduced into the interior of the helix at the 5th position. Thus the modification at the 5th position was targeted for the computational investigation to gain insight into the atomistic level changes upon post-translational modification (Table 2). The Baldwin peptide sequence and the Ser/Thr mutated sequences (Ser/Thr at position 5) were submitted to the I-TASSER server to generate the starting conformations.

50

The subsequent

2-

mutations at the PTM site (-PO3 (SEP/TPO) and O-GlcNAc (S/T)) were introduced using the topology information present in the CHARMM topology file onto the unmodified Baldwin helix structures (Ser and Thr) generated from I-TASSER. For all the modified peptides, Ser, Thr and their phosphorylated and glycosylated counterparts, the N-Cα-Cβ-OG (χ) dihedral was initially constrained to 60° in accordance with PDB survey and NMR studies of model dipeptides. 14, 22, 28 Table 2. Model amino acid sequences used to study the influence of PTM on α-helical and PPII geometries. Sequence

Position

Modification

AKAAAAKAAAAKAA

---

Unmodified

AKAAXAKAAAAKAA

5

Ser, Thr, -PO32- (SEP/TPO), -O-GlcNAc (S/T)

KTPP

2

-PO32- (SEP/TPO), -O-GlcNAc (S/T)

KSPP

2

-PO32- (SEP/TPO), -O-GlcNAc (S/T)

Similarly to examine the influence of the PTM on extended PPII geometries, model peptides from the proline rich domain of tau protein (microtubule binding protein), KTPP and KSPP units



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

present in tau174-183 (K174TPPAPKTPP183) were studied.

28

Page 12 of 53

The proline rich domain of tau links

the N-terminal sequence and the tubulin-binding domain. It has been reported that hyperphosphorylation in the proline rich regions leads to tau aggregation. .24 The experimental study by Zondlo et. al. targeted four units in the proline rich domain, namely tau174-183, tau196-209, tau211219

and tau229-238 for their experimental study but we have restricted the study only to KTPP and

KSPP units from the tau174-183 due to the availability of NMR 3J-coupling data for these sequences. Three Karplus equations were used to estimate the φ (C-N-Cα-C), χ1 (N-Cα-Cβ-OG) and χ2 (Cα-Cβ-OG-P) dihedrals from the J-coupling data, 51-53 3

JαN = 6.51 cos2 (φ-60) - 1.76 cos (φ-60) +1.6

(1)

3

(2)

3

(3)

JHαHβ = 4.37 – 1.86 cos (χ1-120) + 3.81 cos2 (χ1-120) – 0.37 sin (χ1-120) JPH = 15.3 cos2 (χ2-120) – 6.1 cos (χ2-120) + 1.6

For generating the initial conformations the dihedral values back calculated from the J-coupling data available in the experimental literature were used to constrain the respective dihedrals. The residue specific dihedral values used to constrain φ are provided in Table S2 of the SI file. χ1 (NCα-Cβ-OG) and χ2 (Cα-Cβ-OG-P) were constrained near -60° and 120o in accordance with the Jcoupling data available for the KTPP and KSPP peptides. The remainder of the topology information present in the force field was used to generate the staring conformations. The previously described simulation protocol to study the model dipeptide was used for minimization and equilibration. Similar to the model dipeptide systems the equilibrated structures were used to initiate 100 ns production simulations performed with GROMACS version 5.1.2 using the same protocols as described earlier. The dimensions of the systems are reported in Table S3 of the SI file. Result and Discussion PDB analysis φ/ψ Analysis: To determine the influence of the PTM on the local conformational preference we analyzed the φ/ψ distributions from the curated PDB datasets. For each PTM site i, φ (Ci-1-NiCαi-Ci) and ψ (Ni-Cαi-Ci-Ni+1) dihedrals were evaluated, where i corresponds to the modified 12 ACS Paragon Plus Environment

Page 13 of 53

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


(SEP/TPO) or unmodified (Ser/Thr) amino acid based of the dataset (Figure 1). The obtained φ/ψ distributions were classified into the well-defined Ramachandran regions, αR-helix (φ = 160o to -20o, ψ = -120o to 50o), β-sheet (φ = -180o to -90o, ψ = 110o to 180o), PPII (Polyproline II helix) (φ = -90o to -20o, ψ = 110o to 180o) and αL helix (φ = 20o to 160o, ψ = -50o to 120o).

54

The

populations corresponding to these regions are tabulated in Table 3. Table 3. φ/ψ distribution statistics from PDB data. αR-helix β-sheet PPII αL-helix Outliers Total SEP

62

31

31

2

52

179

SER

32

11

38

0

11

92

SER’a

551

154

203

30

40

978

TPO

18

16

50

2

18

106

THR

17

8

4

2

3

34

91

59

4

38

415

THR’a 223 a

SER’ and THR’ corresponds to all the serine and threonine sites present in the PDB database.

The analysis reveals significant differences between the conformation preferences of phosphorylated serine and threonine (Figure 2a and 2b). SEP is found to prefer the compact αhelical conformation (35%) while TPO prefers the extended PPII conformation (47%). The propensity of SEP to adopt the extended β-sheet and PPII conformations turns out to be 17% each, while the propensity of TPO to adopt the compact α-helical conformation and the extended β-sheet conformation turn out to be 17% and 15% respectively. In contrast analysis of the unmodified sites reveals that serine prefers both the compact α-helical conformation (35%) and extended PPII conformations (41%), while threonine predominantly prefers the compact αhelical conformation (50%) to the extended β-sheet conformation (24%). The results and the classification must be taken into account with the knowledge that there is significant reduction in the size of the datasets while comparing the modified and unmodified sites. The reduction in the dataset is about 50% on going from SEP to serine, while the reduction is close to 68% on going from TPO to threonine. To quantify the conformation preference of serine and threonine we also analyzed all the serine and threonine residues present in the PDB database. This resulted in 978 13 ACS Paragon Plus Environment


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 53

serine and 415 threonine unique sites. Our analysis reveals that both serine (56%) and threonine (53%) prefer the compact α-helical conformation. This observation is in line with earlier observations of intermediate α-helix propensities for Ser and Thr. 55 PDB analysis reveals a clear structural transition upon phosphorylation with the effects being significant for threonine when compared to serine. A clear structural transition from compact αhelical structures for unmodified threonine residues (50%) to extended PPII geometry for modified TPO sites (47%) are seen (Figure 2c). To gain insight into the intra and inter residue interactions responsible for the structural transitions we performed close contact analysis for all the structures present in the PDB datasets. Close Contact Analysis: The introduction of the charged phosphate group (PO32-) at serine/threonine residues upon PTM introduces potential sites for nonbonding interactions. The oxygen atoms of the phosphate group can be involved in various interactions like; H-bonding and salt bridges.

13

To gain insight into the non-bonded interactions a distance cutoff of 3.0 Å

(between heteroatoms) from the oxygen’s of the phosphates was used to identify the H-bonding partners while a distance cutoff of 4.0 Å was used to identify salt bridge interactions. Three nonbonded interactions were identified from the analysis, (i) intra-residue H-bonding between the phosphate oxygen’s and the amide hydrogen of SEP/TPO, (ii) salt bridge interactions with proximal Arg/Lys, and (iii) intra-residue H-bonding interaction between the phosphate oxygen’s and the alcoholic side chain of Thr two amino acids (i-2) from the PTM site. The interactions were further classified based on the structural class of the underlying peptide, αR-helix, β-sheet or PPII. The results of the close contact analysis are presented in Table 4. Table 4. Close contact analysis statistics for SEP and TPO. Structures are classified based on Intra-residue H-bonding, Salt-bridges and Thr (i-2) H-bond. SEP

α-helix

β-sheet

PPII

Total Structures

62

31

31

Intra-residue Hbonding

17

1

4

Salt-bridge

12

1

4


Page 15 of 53

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


Thr (i-2) H-bond

3

1

1

TPO

α-helix

β-sheet

PPII

Total Structures

18

16

50

Intra-residue Hbonding

9

4

34

Salt-bridge

5

4

24

Thr (i-2) H-bond

2

1

7



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 53

Figure 2: Ramachandran φ/ψ analysis for (a) SEP, Ser and Ser’ and (b) TPO, Thr and Thr’ from the redundancy removed datasets. (c) Relative population distributions in the α-helix, β-sheet and PPII-helical regions for SEP/Ser/Ser’ and TPO/Thr/Thr’. Ser’ and Thr’ corresponds to all the


Page 17 of 53

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


serine and threonine sites present in the PDB database. The φ/ψ plots were created using MolProbity. 56 Intra-residue hydrogen bonding: The di-anionic phosphate group in SEP/TPO is involved in intra-residue H-bonding with the amide proton of the underlying serine/threonine. Interestingly we find that the probability of H-bonding is more in the case of TPO (43%) when compared to SEP (12%). An analysis of the underlying structural classes (α-helix, β-sheet and PPII helix) reveals that TPO in the PPII form favors the formation of intra-residue H-bond, with the H-bond being present in 68% of the PPII structures. In contrast the H-bond is not very prevalent in the SEP structures. The H-bond locks the conformation of the phosphate with respect to the protein backbone (Figure 3a). Cho and co-workers reported differential intra-residue phosphate-amide H-bonding preference between SEP and TPO using a combination of NMR and vibrational circular dichroism (VCD) studies. 57-58 Using phosphorylated model di-peptides constructs they found that the phosphorylated constructs favored the PPII and β-strand conformations. These observations are reflected in the PDB statistics with the extended conformation favoring the formation of the intra-residue H-bond. Cho and co-workers also found that the phosphate group could form a H-bond with the ith and the i+1th amide hydrogen for phosphorylated threonine while only with the ith amide hydrogen for phosphorylated serine. However, PDB analysis reveals a strong preference for the intra-residue H-bond with the ith amide hydrogen for both phosphorylated serine and threonine. Arg/Lys salt bridge interaction: In addition to the intra-residue H-bond a salt bridge was observed between Arg/Lys side chains and the di-anionic phosphate side chain of SEP/TPO. In fact stabilization by phosphorylation, utilizing salt bridge interactions has been used earlier in the design of stable four-helix bundles, leucine zippers and improving α-helix stability.

20, 23, 59

Interestingly a tailored phosphorserine-lysine salt bridge could stabilize the α-helical geometry by about -2.0 kcal/mol.

23

The PDB analysis revealed that salt-bridges are principally observed in

the PPII form of TPO (48%) and to a lesser extent in the α-helical form of SEP (19%). A closer inspection of the salt bridge interaction revealed no sequence specificity but rather a regiospecificity in the interaction. All the TPO-Arg/Lys salt bridge interactions were observed in the activation loop of kinases. In Figure 3b we present the overlay of all the kinase activation loops exhibiting the Arg/Lys-Phosphate interaction. This interaction has been recognized to be crucial 17 ACS Paragon Plus Environment


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 53

towards the kinase activity with mutation of the Thr197 in the activation loop to Ala resulting in a significant loss of function.30,

60

In fact the conformational changes brought about by

phosphorylation in the activation loop of CDK2 (cyclin-dependent kinase protein) were found to be crucial towards the enzymatic activity.61 An analysis of the PDB dataset containing TPO and SEP residues reveals that 48% of the structures belonged to the protein kinase family of different classes such as MST3, IRE1, Aurora, PRK, MAP, AMPK, PK, CDK, Y81, PIK1 and PKN2 (Figure 3b). Two modes of phosphate-Arg/Lys interactions were primarily observed in the crystal structures, the first one being a bi-dentate interaction between Arg and the phosphate with an additional supporting mono-dentate interaction between the Lys and the phosphate (Figure 3c). In the second case no assisted Lys interaction was observed in the crystal structures (Figure 3d). i-2 Thr interaction: Interestingly the PDB analysis also revealed a previously unrecognized and unreported sequence dependent (i, i-2) H-bonding interaction between the alcoholic side group of Thr and the phosphate side chain of TPO (Figure 3e). This interaction was structurally observed in only 14% of the TPO PPII structures. However it is important to note that out of 52567 experimentally verified phosphorylated threonine sites (dbPTM statistics,

62

described

later in the text) structural information is available for only 605 phosphorylated threonine sites (PDB statistics), which amounts to only 1% of the experimentally verified sites. Since the interaction was sequence dependent (Ti-2XTip), where Tip corresponds to the phosphorylated threonine site and X is any amino acid, we attempted to explore the presence of TXTP and TXSP (TP and SP corresponding to phosphorylated threonine and serine sites) motifs within the dbPTM database which is a repository of post-translational modification sites. dbPTM contains 62346 and 52567 experimentally verified phosphorylated serine and phosphorylated threonine sites respectively. The TXTP motif was found in 3548 PTM entries while the TXSP motif was identified in 3652 entries. This corresponds to an occurrence rate of 8% and 6% respectively, which is significant corresponding to the number of crystallized structures, which is only 1%. Analysis of the intervening residue X revealed a preference for serine, with TSTP and TSSP being observed in 401 and 475 entries out of 3548 and 3652 entries respectively.


Page 19 of 53

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


Figure 3: (a) Representative intra-residue hydrogen bonding between TPO phosphate and the amide backbone from PDB id: 1APM,63 (b) Left: Overlay of the kinase activation loops exhibiting Arg/Lys-Phosphate interactions, Right: Classification of the protein structures containing SEP and TPO residues. Representative phosphate-Arg/Lys binding modes in kinases (c) Bi-dentate interaction between Arg-phosphate supported by mono-dentate interaction between Lys-phosphate (PDB id: 1APM),63 (d) Bi-dentate interaction between Arg-phosphate without any additional Lys interaction (PDB id: 4W8D).

64

(e) Representative inter-molecular

hydrogen bond between TPO phosphate (i) and threonine (i-2) side chain hydroxyl (PDB id: 1APM).63 DFT Studies The analysis of the PDB datasets revealed clear structural trends for the phosphorylated systems, with TPO preferring extended conformations as against compact geometries. While for SEP an opposite trend was observed with the compact α-helical geometry being preferred over the extended conformations. An analysis of all the serine and threonine sites in the PDB revealed a



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 53

sight preference for the compact α-helical structures for both serine and threonine. We note that the PDB datasets only account for a small percentage (1%) of the experimentally verified PTM sites. Additionally the PDB database does not contain statistically significant O-GlcNAcylated structures to comment on the structural trends upon glycosylation. Thus to gain insight into the structural consequences upon phosphorylation and glycosylation on an equal footing we performed extensive QM calculations at the model dipeptide level (Figure 1). QM scans using the model dipeptides, namely serine and threonine along with their phosphorylated and glycosylated counterparts were performed to scan the entire Ramachandran φ/ψ space. The 2D energy profiles obtained from the scans are presented in Figure 4. In total 3456 (576 x 6) conformations were used to generate the QM energy maps. The statistics of the global minimum conformation and the relative energy difference between the various major conformations (αhelix, β-sheet, PPII) and the global minimum are presented in Table 5. A preview of the QM scans reveals that for all the systems other than O-GlcNAc (S) there are two well-defined energy minimums. In Table 5 we also present the φ/ψ values and the relative energies corresponding to the 2nd energy minimum and the unrestrained minimum obtained when the φ/ψ dihedral constraints were removed. Our results are in agreement with previous gas phase QM studies for the Ser and Thr dipeptides that find the global energy minimum to reside in the γ-turn region, an intermediate space between the β-sheet and α-helix regions.

35, 43

Additionally previous QM studies have shown that

the global minimum shifts to Ramachandran allowed regions in the presence of solvation effects.35,

65

We have explored solvent effects and their influence on the φ/ψ distribution using

MD simulations, which are discussed later in the manuscript. The discussion on the QM energy scans allows us to compare the conformational energetics of all the dipeptides on the same footing. Table 5. φ/ψ values corresponding to the global minimum, 2nd minimum and the unrestrained minimum from the QM energy maps. Relative energy difference between the various major conformations (α-helix, β-sheet, PPII) and the global minimum. Energies are in kcal/mol. Dipeptide

φ/ψ

φ/ψ

(Global

(2nd

Relative free energy for canonical conformations with respect to the

φ/ψ


Page 21 of 53

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


minimum) minimum)

(unrestrained) global minimum α-helix

β-sheet

PPII-helix

(-600/450)

(-1350/1350)

(-750/1500)

-1650/1800

-900/750 (0.3)

-1650/1800

12.8

3.3

5.6

Threonine -1500/1800

-900/600

-1550/1780

13.8

3.5

5.7

Serine

(1.9)

SEP

-750/-300

-750/750 (1.5)

-770/-270

0.6

13.8

7.1

TPO

-750/-300

-750/-600 (1.4)

-1150/-450

0.2

13.1

7.0

OGlcNAc (S)

-900/600

--

-820/640

10.1

5.8

6.0

OGlcNAc (T)

-1050/00

-1200/1650 (3.0)

-1120/110

6.9

4.7

4.7

Serine and Threonine: In Figure 4a and 4d we present the QM energy maps corresponding to the serine and threonine dipeptides. An analysis of the global energy minimum regions reveals that both serine and threonine prefer the extended β-sheet conformation with the φ/ψ values corresponding to the energy minimum being -1650/1800 and -1500/-1800 for serine and threonine respectively. While this corresponds to a canonical Ramachandran region the second minimum lies at -900/750 and -900/600 which corresponds to a 2.27 ribbon (γ-turn, φ/ψ = -780/590) conformation, a non-canonical region in the Ramachandran space.

66

The energy surface is found

to be relatively flat for serine with the relative energy difference between the two minimums being 0.3 kcal/mol, while for threonine this energy difference is 1.9 kcal/mol. These differences between serine and threonine correlate with the steric influence of the γ-methyl group in threonine, which restricts its conformational flexibility.



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 53

SEP and TPO: In Figure 4b and 4e we present the QM energy maps corresponding to SEP and TPO dipeptides. In contrast to serine and threonine dipeptides SEP and TPO dipeptides prefer the canonical α-helical region with an energy minimum located at -750/-300. This observation is partially in line with the PDB survey data wherein SEP was found to prefer the compact αhelical (35%) geometry to the extended β-sheet (17%) geometry. Similar to the serine and threonine dipeptides the second energy minimum for SEP (-750/750) and TPO (-750/600) lies in the γ-turn region with the relative energies being 1.5 kcal/mol and 1.4 kcal/mol respectively. O-GlcNAc (S) and O-GlcNAc (T): In Figure 4c and 4f we present the QM energy maps corresponding to O-GlcNAc (S) and O-GlcNAc (T) dipeptides. As mentioned earlier for OGlcNAc (S) we observe only one energy minimum that is centered at -900/600 corresponding to the γ-turn region of the Ramchandran surface. The energy surface is relatively flat towards the extended conformation regions of the φ/ψ space with the canonical β-sheet and the PPII helix geometries being 5.8 and 6.0 kcal/mol higher in energy compared to the global minimum (Table 5). In contrast the global minimum is shifted towards compact conformations for O-GlcNAc (T). The minimum lies at -1050/00, which in principle does not correspond to canonical regions on the φ/ψ space. The second minimum at -1200/1650 lies closer to the extended β-sheet regions.


Page 23 of 53

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


Figure 4: 2D-dihedral potential energy scans about the φ/ψ dihedrals for the model dipeptides. Top panel: Ser, SEP and O-GlcNAc dipeptides; Bottom panel: Thr, TPO and O-GlcNAc (T) dipeptides. Energies are in kcal/mol, with contours every 1.75 kcal/mol. Only energies below 14 kcal/mol have been plotted for the sake of clarity. Intra-residue H-bonding: As observed in the PDB analysis, the closely spaced H-bonding donors and acceptors in the studied dipeptides can be involved in intra-residue H-bonding. To explore the possible H-bonding interactions and in turn their influence on the overall energy profile we analyze the distance distribution between the H-bonding partners. The most prevalent H-bond interaction among all the dipeptides was found between the i-1 carbonyl oxygen and the i+1 amide hydrogen (OYi-1…HNTi+1) of the dipeptide. This H-bond is known to stabilize the C7 conformation of the dipeptides and has been observed in earlier gas phase theoretical studies using model dipeptides.35 Analysis of the distance distribution corresponding to this H-bond reveals that the H-bond is relatively stronger for the phosphorylated dipeptides with a larger spread (1.50 Å – 2.50 Å) in the φ/ψ surface when compared to other dipeptides (Figure S1 of the SI file). However on comparison with the overall energy profile we observe that the global minimum and the low energy surfaces do not correlate with the C-Oi-1…H-Ni+1 H-bonding profile indicating that this might not be the significant H-bonding interaction. For Ser and Thr dipeptides the side chain alcoholic hydroxyl is involved in an intra-residue Hbond with the ith carbonyl of the amino acid (HG1…O) (Figure S2(a) and Figure S2(b) of the SI file). The H-bonding interaction is observed in regions that correspond to the low energy conformation for Ser and Thr in the γ-turn regions of the φ/ψ space indicating that this is a stabilizing interaction. However we note that owing to the geometrical constraints of the side chain the H-bond is relatively weak with the HG1...O angle being in the range 1300-1400 for both Ser and Thr dipeptides. Interestingly such conformers have been detected in the conformational landscape of serine in the gas phase by Fourier transform microwave spectroscopy. 67 For SEP and TPO there is a very strong H-bond between the phosphate oxygen’s and the ith amide hydrogen (OP…HN) (Figure S2(c) and Figure S2(d) of the SI file). The occurrence of the H-bond also corresponds to the observed minimum in the energy surface. However we note that in our DFT scans we constrained the Cα-Cβ-O-P dihedral to 1100 in accordance with NMR data,



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

57-58

Page 24 of 53

which also geometrically constrains the phosphate group relative to the amide hydrogen that

also reflects in the really large spread of the H-bond. We comment more on the influence of flexibility of the phosphate group and the corresponding H-bond in detail from the unrestrained MD simulations of the dipeptides. For O-GlcNAc (S) and O-GlcNAc (T) the presence of the acetlyamide (-NHCOCH3) group on the carbohydrate introduces additional H-bonding sites. H-bonding between the amide hydrogen of the acetlyamide group and the backbone carbonyl oxygen has been observed earlier in targeted MD studies, QM calculations and NMR studies.

14

We observe two different patterns to

H-bonding in O-GlcNAc (S) and O-GlcNAc (T) from the QM calculations (Figure S2(e) and Figure S2(f) of the SI file). For O-GlcNAc (S) the HNcarb…O H-bonding (subscript carb indicates amide hydrogen from carbohydrate) profile corresponds to the observed QM energy surface while it is not the case for O-GlcNAc (T). Interestingly the free hydroxyl at the C6 carbon on O-GlcNAc (T) is involved in H-bonding with the carbonyl oxygen of the underlying peptide (HO6…O) which significantly influences the energy profile (Figure 5). The H-bonding is observed in the unusual minimum at -1050/00, which in principle does not correspond to a canonical Ramachandran region. We note that the propensity of this H-bond is significantly lost in the solvated MD simulations as the free hydroxyl is solvated. However it is interesting to note that the presence of the H-bond was observed in targeted MD simulations. 14

Figure 5: 2D distribution of the HO6…O H-bond distances as a function of φ/ψ dihedrals for OGlcNAc (T). The minimum energy geometry of O-GlcNAc (T) corresponding φ/ψ value of -


Page 25 of 53

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


1050/00. The distances are presented in Å, with contours every 0.25 Å. Only distances between the range 1.25 Å - 3.0 Å have been plotted for the sake of clarity. MD Studies To study the influence of solvation on the structural propensities of phosphorylated and glycosylated Ser and Thr we performed MD simulations for model dipeptides, Baldwin helices (α-helical structures) and segments of the Tau protein (PPII structures) as described in the experimental section. First we discuss the results from the dipeptide simulations. Dipeptide Simulations φ/ψ distributions: In Figure 6 we present the relative free energy as a function of φ/ψ for each dipeptide system calculated from the extended simulations (3 x 100 ns for each system). In Table 6 we also present the φ/ψ values corresponding to the global minimum and the second relative minimum in the free energy surface. The relative energies of the canonical conformations with respect to the global minimum along with the percentage populations in the corresponding canonical regions are also tabulated in Table 6, the data is also presented in Figure S3 of the SI file. We observe a clear shift towards canonical structures in the MD simulations when compared to the DFT studies, an effect that has been observed earlier. The shift towards canonical structures is primarily due to solvation affects that are captured by MD simulations.

35, 65

Both Ser

and Thr show a strong preference for both the α-helical and PPII geometries with the α-helical conformations being slightly preferred. The relative energy difference between the two being only 0.7 kcal/mol and 0.5 kcal/mol for Ser and Thr respectively, which is close to the Boltzmann’s thermal energy at room temperature (0.59 kcal/mol). Interestingly for the phosphorylated systems, SEP preferred the extended PPII conformation to the compact α-helical conformation, while TPO preferred the α-helical conformation to the PPII conformation. It is to be noted that this conformational preference is contrary to the conclusion arrived at from the PDB statistics, wherein SEP preferred the α-helical conformation while TPO preferred the PPII geometry. As noted earlier one of the factors that might be responsible for this observed discrepancy is insufficient structural information in the PDB database as structural information is available for only 1% of experimental verified PTM sites, which is statistically low. Additionally we note that the discrepancies could also arise from the influence of the neighboring residues in 25 ACS Paragon Plus Environment


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 26 of 53

the PDB structures and biases in classical force fields. For both O-GlcNAc (S) and O-GlcNAc (T) we see a remarkable shift in the population towards extended structures with the PPII conformation being most favorable. A significant decrease in the overall population in the αhelical region (11%, 17%) is seen for the O-glycosylated dipeptides when compared to the native dipeptides (46%, 47%) and phosphorylated dipeptides (40%, 54%) indicating a significant structural bias upon glycosylation. To gain insight into the observed structural trends we analyze the simulation trajectories for close contacts. Table 6. φ/ψ values corresponding to the global minimum and the 2nd minimum from the MD free energy maps. Relative energy difference between the various major conformations (α-helix, β-sheet, PPII) and the global minimum. Percentage populations of the major conformational regions are provided in parentheses. Energies are in kcal/mol. Dipeptide

φ/ψ

φ/ψ

(Global minimum)

(2nd minimum)

Relative free energy for canonical conformations with respect to the global minimum (% population)

α-helix

β-sheet

PPII

(-600/-450)

(-1350/1350)

(-750/1500)

Ser

-900/-500

-900/1650

0.0 (46)

2.4 (20)

0.7 (19)

Thr

-900/-500

-900/1650

0.0 (47)

2.9 (22)

0.5 (23)

SEP

-900/1650

-900/-750

1.5 (40)

3.9 (15)

0.0 (26)

TPO

-900/-750

-900/1650

0.0 (54)

2.7 (14)

1.3 (23)

OGlcNAc (S)

-450/1650

-1650/1650

2.1 (11)

2.5 (28)

0.0 (50)

OGlcNAc

-450/1650

-750/-450

2.0 (17)

2.6 (22)

0.0 (59)


Page 27 of 53

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


(T)

Figure 6: 2D relative free energy surfaces for φ/ψ dihedrals, given in degrees, calculated from MD simulations for all the dipeptide systems. Relative free energies are calculated from the natural logarithm of the relative probability and are given in kcal/mol. Energies are in kcal/mol, with contours every 0.5 kcal/mol. Only energies below 4 kcal/mol have been plotted for the sake of clarity. Intra-residue H-bonding: Similar to the PDB analysis and the QM studies we analyze the MD simulation trajectories for significant intra-residue H-bonding interactions. The most significant H-bonding interactions were found in the phosphorylated di-peptide systems. For both SEP and TPO a strong H-bond was observed between the phosphate oxygen’s and the underlying amide hydrogen’s i(HN) and i+1(HNT) (Table 7). For SEP the H-bond was found to be predominantly between PO32- and HN, while for TPO the H-bond was found to be prevalent between both PO32-HN and PO32--HNT (Figure 7a). The average H-bonding distance was found to be around 1.8 Å. The φ/ψ distribution corresponding to the H-bonded (dH-bond

Phosphorylation versus O-GlcNAcylation: Computational Insights into

Recommend Documents