Structural modulation of p53TAD1-TAZ2 complex upon mutations and post-translational modification

Abstract The tumour suppressing p53 is a target for genetic alterations in human cancer. Native p53, found in latent state in cells, gets activated following various intracellular or extracellular responses. It plays imperative role in cell-cycle control, via growth-arrest, DNA repair and apoptosis, mainly regulated by post-translational modifications (PTM). However, the influence of PTMs on the activity of p53 is still under extensive experimental and computational study. There are numerous PTM sites in p53, which are reported to regulate its binding affinities with other proteins. Of the many, Thr18 at transactivational domain (TAD) of p53 is reported to amplify p53 activity upon phosphorylation. To understand the molecular basis of p53 recognition by its binding partner upon mutations and PTMs, we have exploited all atom molecular dynamic (MD) simulation of p53TAD1 bound to TAZ2 domain of p300. The MD simulation inferred that phosphorylated and mutated Thr18, as a phospho-mimic, bound with TAZ2, redistributed the charge environment of the interface, thereby modulating the stronger interactions with TAZ2 to enhance the binding efficiency. The electrostatic interactions due to different charge environment together with H-bonding and hydrophobic interaction dictate diverse binding approach between the two. The results of this computational study further explain the importance of the Thr18 as a PTM site in atomistic detail, hence shedding further light to the understanding of how PTMs are imperative for p53 activity to protect the cellular world. Communicated by Ramaswamy H. Sarma

The functioning of p53 is vastly regulated by post-translational modifications (PTM), e.g., phosphorylation, acetylation, ubiquitination etc (Appella & Anderson, 2001;Ithuralde & Turjanski, 2016;Uversky et al., 2008). Although the PTM sites are mainly present on the sequence-specific DNA-Binding domain (core domain), a few are also found in N-and C-terminals of p53 (Bode & Dong, 2004). PTM has substantial importance in protein activity. Phosphorylation on serine and threonine has actively contributed to the regulation of cellular processes including signal transduction, apoptosis, cell growth etc. Also, phosphorylation at N-terminal domain of p53 has impact on the activity of its other domains. In normal cells the level of p53 is kept low by continuous degradation mediated by MDM2 binding (Ozaki & Nakagawara, 2011). Cellular stress, e.g. DNA damage, induces phosphorylation at multiple sites on the entire p53, which weakens MDM2 binding and hence increases p53 level in the cell. It has been shown earlier that the phosphorylation of p53 further plays a rather decisive role in regulating its interactions with other binding partners such as p300 (Feng et al., 2009;Ithuralde & Turjanski, 2016), which contributes to the stabilization and transactivation function of p53 following DNA damage. Mutations at specific PTM sites in p53 have been reported to modulate its secondary structure and alter binding affinity to its binding partners (Nikolova et al., 2000).
In this study we focused on the trans-activation domain 1 (TAD1, residue range: 1-39) of p53 (henceforth will be referred as p53TAD1). This domain is of immense biological importance because of its interaction with transcriptional coactivators and co-repressors. Among the main PTM sites at p53TAD1, Thr18 is the one, which may control the binding of other regulatory proteins. Importantly phosphorylated Thr18 weakens MDM2 binding to activate p53 function in cell (Ithuralde & Turjanski, 2016;Sakaguchi et al., 2000). Earlier simulation and experimental studies revealed that p Thr18 decreases the MDM2 binding affinity by strong anionic charge repulsion at the vicinity of the binding site, while keeping the helix mostly undisturbed (Brown et al., 2008;Lee et al., 2007). Feng et. al. (Feng et al., 2009) has reported that phosphorylation of p53TAD1 enhances its binding with TAZ2, the C-terminal domain of p300, by eleven times than the un-phosphorylated one. Furthermore, binding of p53TAD1 due to phosphorylation is more enhanced in p Thr18 than in p Ser15. To shed more light on investigating the imperative role of p Thr18 on TAZ2 binding, we have used molecular dynamics simulation to find out how p53TAD1 has adjusted itself to enhance the binding after phosphorylation to continue its role towards transactivation. Phosphorylation induces extra negative charge. The charged environment at the binding site upon introducing phosphate group may play pivotal role while binding. Recent MD simulation has reported the lists of newly formed interactions to stabilize the peptide (Ithuralde & Turjanski, 2016). To learn the effect of phosphorylated Thr18 to TAZ2 binding, we studied p53TAD1-TAZ2 complex with both non-phosphorylated and phosphorylated Thr18 using molecular dynamics simulation. For further understanding the role of charge environment at the binding site other than phosphorylation we mutated Thr18 by phospho-mimic Glu. Furthermore, Thr18 was mutated by non-bulky, chemically inert Ala to study the unique role of Thr or phosphorylated Thr towards the activity of p53 upon binding.

Materials and methods
All-atom, explicit solvent molecular dynamics simulation has been used to study the importance of phosphorylated Thr18 of p53TAD1 towards TAZ2 binding. The NMR model of the complex (Feng et al., 2009) has been downloaded from Protein Data Bank (PDB code: 2K8F) (www.rcsb.org). Among the 10 NMR models available for the particular complex, the first model has been selected. Thr18 of p53TAD1 was phosphorylated and mutated using CHARMM GUI server (Jo et al., 2008(Jo et al., , 2014. The server has used CHARMM36m force field (Huang et al., 2017) and parameter to generate whole systems along with phosphorylated threonine. In phosphorylated form threonine exists as dianionic phosphothreonine ( Figure S1). Four constraints were prepared as follows ( Table 1).
The energy-minimized structure (ABNR algorithm) of individual p53TAD1-TAZ2 complex was solvated explicitly with TIP3P water model in a cubic box of dimension 102 Â 102 Â 102 Å 3 with periodic boundary conditions, and ionic strength of 0.05 M with KCl. The ionic concentration was kept same as found in the ITC experiments reported by Feng et. al. (Feng et al., 2009) at which the tighter binding was reported. The non-bonded long-range electrostatic interactions were calculated using Particle Mesh Ewald (PME) (Darden et al., 1993;Essmann et al., 1995). Cutoff distance was set at 12 Å. SHAKE has been used to fix the hydrogen bond distances during the simulation. Integration time step was set at 2 fs and coordinates were saved after 2 ps interval. WT, p Thr18, Thr18Ala and Thr18Glu were equilibrated for 100 ps followed by production run of 200 ns each, at 310 K. Simulation was performed using NAMD version 6.7.0 (Phillips et al., 2005) with CHARMM36m force field (C36mFF) (Huang et al., 2017) which is proven to be better at handling the simulation of disordered proteins. C36mFF is an improved version of CHARMM36 force field (Huang & MacKerell, 2013). An additional independent trajectory of each of the four systems was triggered with the equilibrated coordinates and simulated for 200 ns using different initial velocities while keeping other simulation parameters identical. Thus each system has two independent trajectories (Traj1 and Traj2) of 200 ns each, totaling 400 ns. Most of the analyses were done with CHARMM version 40 (Brooks et al., 1983) using the last 100 ns data for both the trajectories for all systems. Principal Component Analysis (PCA) was performed using the last 100 ns data of Traj1 of every set. PCA calculations and figures were done using PYTRAJ and PYPLOT of MATPLOTLIB (https://matplotlib.org) (Hunter, 2007;Roe & Cheatham, 2013). Visual Molecular Dynamics (VMD) (Humphrey et al., 1996) was used for viewing and generating the structures in figures. Residue numbers of TAZ2 were kept the same as used in PDB.

Results and discussions
The structure of TAZ2 has remained unaltered throughout the course of simulation. Thus we focused on the structural modulation of the p53TAD1.

Convergence of the simulated trajectories
Several convergence tests, including convergence between independent trajectories and self-convergence of individual trajectories were made to check the reliability of the simulation data. We have done two 200 ns independent simulations of each system. The first convergence test was made using representative structural properties between independent trajectories. TAZ2 bound p53TAD1 forms a short a-helix of the range of Ser15-Pro27 (Feng et al., 2009), whereas unbound p53, being intrinsically disordered, adopts highly flexible conformations with no trace of secondary structure ( Figure S2). Figure 2 compared the average helical propensities calculated from the two independent trajectories of WT, p Thr18, Thr18Ala and Thr18Glu.
The overlapping traces in Figure 2 indicate well-converged conformational ensembles calculated from last 100 ns of each of the two independent trajectories Traj1 and Traj2.
Representative structural properties over the course of simulations have been used to check self-convergence of the independent trajectories too. For all the four systems, residue helicity ( Figure S3) and C a RMSF ( Figure S4) were calculated using the last 120 ns of the simulation timescale.
Our all atom explicit solvent models and MD simulations produced well-converged simulation data to study the importance of Thr18 on TAZ2 binding. The trajectory has been converged after the initial 80 ns simulation. It has been seen from Figures S3 that the primary helix range for all the systems remained same for both the trajectories. Trajectories Table 1. Sequences of four constraints under study are listed. Phosphorylated and mutated residues are indicated as bold and underlined.

Constraints
Sequence MEEPQSDPSV EPPLSQEEFS DLWKLLPENN VLSPLPSQA showed occurrence of a secondary helix, but infrequently (data now shown), thus ignored. We have also calculated helical content ( Figure S5) and potential energy ( Figure S6) of p53TAD1 along the entire 200 ns simulation of both Traj1 and Traj2 and compared to check any significant fluctuations, which may indicate instability in the system. From the figures it is quite evident too that the systems are stable and converged. The fluctuations in the charged systems ( p Thr18 and Thr18Glu) are due to the charge interaction, which is present in almost equal intensities in both the trajectories. Overall simulation seems to converge well based on the calculated one-dimensional distribution of various average structural properties.
Total 200 ns of the simulation data, last 100 ns from each of the two independent trajectories, have been used to further calculation and analyses.

Structural diversion
Phosphorylation at Thr18 enhances binding affinity (Ithuralde & Turjanski, 2016). Mutation at this PTM site affects the binding affinity and/or helical fold as reported earlier (Ithuralde & Turjanski, 2016). The modulation in the binding affinity or the helical content of the peptide upon mutation is often attributed to the structural diversion (Ganguly & Chen, 2015). The secondary structure of phosphorylated and two mutated systems have deviated to some extent (Figure 3) from the WT with Thr18Ala showing almost the same helical propensity as WT while that of p Thr18 has shown reduced helical propensity (72%).
The helix (helical range: Gln16-Leu26) in Thr18Ala system is highly stabilized wherein it is found to get attached to the core of its binding partner TAZ2 by strong electrostatic and strong hydrophobic interaction. Here we find striking structural similarities with WT though the interactions leading to the structural changes varies in both the systems. Furthermore, the terminals of p53TAD1 were highly flexible and adopt different conformations, occasionally the terminal residues directly interacting with TAZ2 binding site. The key interactions responsible for the changes in the structure of the systems have been discussed thoroughly later.
The WT system successfully represented the experimentally reported p53TAD1 helicity (helical range: Ser15-Leu26). The residual range of the original helix for WT (Ser15-Leu26) was shifted by one residue in Thr18Ala and Thr18Glu systems while in p Thr18 system, one can find that the helical range has been further reduced (helical range: Gln16-Leu22). During the course of the simulation, helicity of p53TAD1 remained unchanged for all the three systems, WT, Thr18Ala, and Thr18Glu, but that in p Thr18 it has been reduced by 30% that can be justified by the electrostatic interactions which plays a very interesting role particularly in this system.
The probability of end-to-end distance and radius of gyration ( Figures 4A and 4C) were calculated to compare any structural deviation of the complexes.
The high fluctuations in the results of C a -RMSF ( Figure  4D) are mainly due to the terminals of p53TAD1 being highly flexible and can adopt different conformations. Thr18Ala system has the lowest C a fluctuation throughout the peptide chain. This is because of the strong hydrophobic and electrostatic interactions present in the system. This is also evident from the end-to-end distance, radius of gyration and C a -RMSF value of the system. Except the a-helix, rest of the p53TAD1 of WT, p Thr18 and Thr18Glu, were flanking, with Thr18Glu system has the highest fluctuation in the C-terminal region, hence have higher radius of gyration and endto-end distance. Thr18Glu system also has the highest C a fluctuation along the helical region among the other systems ( Figure 4D).
To dig deeper we have also calculated inter-peptide contacts (original and newly made). Contact maps were calculated over the course of the simulation. Last 100 ns of each of the two independent trajectories, totaling 200 ns simulation data, were used to calculate average contacts. If the heavy atoms of any two residues are within 4.2 Å distance (Ganguly & Chen, 2009;Michino & Brooks, 2009), then both residues are found in the contact map in Figure 5. Here, residues 1 to 90 are for TAZ2 while p53TAD1 are numbered from residue 91 to residue 129.
The marked regions show the difference in contacts for the p53TAD1-TAZ2 contacts in the systems. Number of contacts has changed for all systems but Thr18Glu ( Figure 5D) has fewer interactions compared to the other systems.
The fraction of native contacts (Qintra_p53) as a function of time has been calculated for all sets ( Figure S7). The snapshots of the complexes were also shown to describe the structural changes taking place within the simulation time along with to compare upon mutations. Here the native contact is defined as the intra-chain Ca-Ca contact between the region ranging from residues13 to 28 of p53 (little higher than the helical region of the WT as this region has the maximum conformational changes) and the p53 chain with cutoff distance 4.2 Å. The intra-chain contacts showed how the extended helical region of the p53 interacts with the entire p53 chain in the complex. Figure S7 shows the conformational and secondary structural changes of the four sets of the p53 complex. The diverse interactions with TAZ2 were also observed as found in Figure 5.

Hydrogen bonds
We set the cutoff distance between the donor/acceptor atoms (O-HN) at 2.4 Å to define hydrogen bond. The residues of p53TAD1 mainly formed H-Bond with other residues of the same peptide chain. The overall number of hydrogen bonds has remarkably increased in the Thr18Glu system. In case of p53TAD1-MDM2 complex H-bond formed between Asp21 and Thr18 of p53TAD1 stabilized the helix. This Hbond was broken upon phosphorylation of Thr18 thus destabilizing the complex (Sakaguchi et al., 2000). In case of p53TAD1-TAZ2 complex, one can find that Asp21 forms Hbond with residue 18 of p53TAD1 in all the systems. Additionally there is H-bond between Asp21 and Glu17 of p53TAD1 in WT, p Thr18 and Thr18Ala systems. Also for Thr18Glu system we found several H-bonds between the Nand C-terminal residues of p53TAD1 (Table S1), hence vouching for the folded structure of p53TAD1 in the Thr18Glu system ( Figure S8(c)). All the H-bonds are listed in Table S1.

Key interactions
Interaction energy of each residue of p53TAD1 ( Figure 6) reveals that p Thr18 has the highest interaction energy ($547.66 kcal/mol) not only than the other residues of the peptide but also than the other three systems.
The enhanced electrostatic interaction between p Thr18 and binding site residues of TAZ2 is due to the redistribution of the charge environment introduced by the double negative phosphate group (S1 Text). p Thr18 residue (Figure 6) showed the highest interaction energy followed by Thr18Glu. Not only p Thr18, but also the other oppositely charged residues at the binding sites and/or the vicinity of the binding site contributed to the enhanced binding of TAZ2. Arg9, Arg15 and Arg41 residues of TAZ2 (residue numbering are as in the PDB) played key role in electrostatic interaction with negatively charged residues of p53TAD1. Such interactions lead to the conformational changes of p53TAD1 in both mutated and phosphorylated systems ( Figure S8). The major electrostatic and hydrophobic interactions are depicted pictorially in Figure 7. p Thr18 shows less hydrophobic interaction than the other three systems, which is again owing to its changed orientation due to the charged interaction in the system. The detailed descriptions are mentioned in S1 Text.   The electrostatic and hydrophobic interactions are summarized below in Tables 2 and 3.

Principal component analysis
Principal component analysis (PCA) takes the trajectory of the atoms or the ensemble of snapshots of the conformations and deduces linearly independent modes of concerted motions of the atoms, arranged in decreasing order of importance to dominate the overall dynamics of a system. It reduces the huge number of degrees of freedom into the much lower number, so that most of the motions or dynamics can be well represented (David & Jacobs, 2014;Jolliffe & Cadima, 2016). The principal components are the eigenvectors with eigenvalues, calculated from the covariance matrix obtain from the MD simulation. It identifies the most prevailing motion and also the distribution of the motions accessed by a given system (David & Jacobs, 2014;Jolliffe & Cadima, 2016).
PCA of all four systems distinctly computed the nature of the dynamics of p53TAD1 while binding to TAZ2. Interestingly the conformational dynamics changes widely upon mutation. Figure 8 shows the projection of motions of the four systems along the first two eigen vectors, PC1 and PC2 where eigenvalues of the eigenvectors are plotted along the z-axis. The number of frames is color coded where frames are the conformations (snapshots) saved after every 2ps (1000 steps) during the MD Simulation. The projection of PC1 over PC2 has generated scattered plots (Figure 8), which reveal significantly diverse motions of WT, phosphorylated and mutated p53. Figures 8A-D showed particular patterns of motions. There is a similarity in the plots between WT and Thr18Ala Figure 7. Specific interactions between p53TAD1 and TAZ2 at the binding sites. Electrostatic interactions between intra and inter-chain protein segments of the four systems are shown as (A) WT, (B) p Thr18, (C) Thr18Ala, (D) Thr18Glu. Each of p53TAD1 and TAZ2 segments are shown in red and blue in color respectively. Acidic and basic residues of both p53TAD1 and TAZ2 are colored in green and purple respectively.  Thr18  Thr18Ala  Thr18Glu  p53TAD1  TAZ2  TAZ2  TAZ2  TAZ2 Glu2 systems and between p Thr18 and Thr18Glu systems. The motions in both WT and Thr18Ala are very much scattered along both the eigenvectors while for the other two systems, the motion is somewhat restricted. p Thr18 shows restriction in motion for both PC1 and PC2 while for Thr18Ala, it is restricted along PC1 while it shows some motions in PC2. These variances in motion restrictions are the result of the charge interaction as described previously. Strong charge interaction restricts the peptide motion that is reflected in the Figure 8B and little less in Figure 8D.

Conclusion
Well-converged molecular dynamics simulations of phosphorylated and mutated p53TAD1-TAZ2 complex have been used to find the role of an important PTM site Thr18. It has been revealed that phosphorylation and mutations at Thr18 have actually changed the binding interactions of the p53TAD1 by redistributing its charge environment at the binding site and its periphery. Upon introduction of phosphate group the same charge repulsion between the p Thr18 and Glu17 of p53TAD1 forced to destroy the triangular arrangement of the charge residues Glu17 and p Thr18 of p53TAD1 and Arg9 of TAZ2. The redistribution of the negative charges and rearrangements of the residues in the proximity of the binding site enhance the interactions with TAZ2 leaving the helix mostly undisturbed. This is also evident from the highest interaction energy of p Thr18 ( Figure 6). The newly modified charge environment at the interface facilitates the TAZ2 recognition through both electrostatic and hydrophobic interactions. Such role of Thr18 of p53TAD1 has also been reported in previous study (Brown et al., 2008;Lee et al., 2007) with MDM2, another binding partner of p53, in which the charge repulsion made the binding weaker.  Thr18  Thr18Ala  Thr18Glu   p53TAD1  TAZ2  p53TAD1  TAZ2  p53TAD1  TAZ2  p53TAD1  TAZ2   Phe19  Leu22  Leu25 Interacts with residues inside the helical pockets a2 (Gln37-Gly48) and a3 (Pro58-Cys74).

Disclosure statement
No potential conflict of interest was reported by the authors.
Funding DG acknowledges Department of Biotechnology, Govt. of India for granting Ramalingaswami Fellowship (BT/RLF/Re-entry/52/2012). AG is partially supported from the research grant of the RLS fellowship.