DNA-binding protects p53 from interactions with cofactors involved in transcription-independent functions.

Binding-induced conformational changes of a protein at regions distant from the binding site may play crucial roles in protein function and regulation. The p53 tumour suppressor is an example of such an allosterically regulated protein. Little is known, however, about how DNA binding can affect distal sites for transcription factors. Furthermore, the molecular details of how a local perturbation is transmitted through a protein structure are generally elusive and occur on timescales hard to explore by simulations. Thus, we employed state-of-the-art enhanced sampling atomistic simulations to unveil DNA-induced effects on p53 structure and dynamics that modulate the recruitment of cofactors and the impact of phosphorylation at Ser215. We show that DNA interaction promotes a conformational change in a region 3 nm away from the DNA binding site. Specifically, binding to DNA increases the population of an occluded minor state at this distal site by more than 4-fold, whereas phosphorylation traps the protein in its major state. In the minor conformation, the interface of p53 that binds biological partners related to p53 transcription-independent functions is not accessible. Significantly, our study reveals a mechanism of DNA-mediated protection of p53 from interactions with partners involved in the p53 transcription-independent signalling. This also suggests that conformational dynamics is tightly related to p53 signalling.

INTRODUCTION

Allosteric mechanisms involve the transmission of local structural perturbations, such as those incurred by interactions with ligands, to distal protein regions (1–5). It has been proposed that allostery is an intrinsic feature of proteins (6,7) and is best understood in terms of a shift in the distribution of pre-existing states that are present in both the free and bound protein (8). The effects induced at distant sites can involve major conformational changes (9) or more subtle localized changes in either dynamics or conformation (10). The transmission of structural effects across long distances relies on the dynamic coupling between different residues (3,11–14). In this view, the network of residue contacts and the intrinsic dynamics of the protein are crucial components of long-range communication and protein regulation.

Given the key role that both conformational changes and long-range effects play in the regulation of cellular events (15–17), understanding these events in atomistic detail is an important goal (18,19). Atomic-level molecular dynamics (MD) simulations have previously been employed to elucidate the molecular mechanisms that underlie structural communication over long distances (20). MD can provide information on protein dynamics spanning time-scales from femto- to milli-seconds (21) and long-range communication paths can be inferred from MD ensembles (20,22). Enhanced sampling techniques, such as parallel-tempering metadynamics (PT-metaD) (23), can be used to determine the free energy landscape associated with conformational transitions (24–28). We have recently shown that PT-metaD, when combined with the most accurate protein force fields makes it possible to map the free energy associated with conformational changes with accuracy and precision (28).

Hub proteins, such as the p53 tumour suppressor, regulate fundamental cellular processes through protein–protein interactions and allosteric mechanisms. p53 interacts with multiple partners and integrates diverse signalling pathways to initiate distinct cellular responses (29–31). The ability to exploit the same region for interactions with different biological partners is crucial in hub proteins and is often associated with allosteric modulation (32–34). Although many experimental and computational studies have addressed important aspects of p53-mediated interactions (29,30,35,36), a detailed molecular understanding of the mechanisms underlying long-range structural communication is far from complete. Only a limited number of structures of the complexes between p53 and its biological partners have been determined, limiting the knowledge on the effects of protein–protein interactions on p53 structure and function. Furthermore, since allosteric signalling often involves shifts in the populations of dynamically interconverting conformations that are present in both the unbound and bound states, it is challenging to understand regulation based purely on static protein structures.

The DNA-binding domain (DBD) of p53 is especially important for both cytosolic and nuclear functions of the protein (29,37,38). Even though the p53 DBD does not appear to undergo a large conformational change upon binding with cofactors, a slow conformational exchange process in the proximity of the disordered N-terminal region has been identified by NMR experiments on the free protein (39). Previous studies focused mostly on local changes that occur at the DNA-binding interface and L1 loop (40–44), and much less is known about conformational changes at more distal sites (45).

To gain a deeper insight into conformational rearrangements that underlie p53 regulation and function, we have used all-atom explicit-solvent MD simulations, including PT-metaD enhanced sampling methods (23), to study p53 DBD (Figure 1A) in its DNA-bound (p53DBD-DNA) and unbound (p53DBD) forms. We identified a coupling between conformational changes at the interface for DNA-binding and changes in a loop (S6-S7 loop, residues 207–213, Figure 1B), which is located more than 3 nm away from the DNA-binding and whose conformational ensemble is modulated by DNA. In particular, we found that conformational changes in the L1 loop at the DNA-binding interface are coupled to changes in S6-S7 loop, which in turn is in proximity to the N-terminal disordered tail, for which a slow exchange was observed by NMR (39). We also observed that Aurora kinases-mediated phosphorylation at Ser215 (46–48) modulates the conformational properties of the S6-S7 loop, providing an additional layer of p53 regulation. The coupling between the L1 and S6-S7 loops provides a mechanism for DNA binding to exert long-range effects, thus linking our findings to p53 function. The mechanisms here illustrated highlight one more time the important role of loops in regulation of protein function (5). In particular, by comparing simulations in the presence and absence of DNA, we observed a population shift upon DNA binding toward p53DBD conformations that disfavour the interaction with partners involved in p53 transcription-independent apoptotic functions, such as Ku70.

Figure 1.

(A) Structural features of the p53 DBD. The 3D structure of p53 DBD in complex with DNA (PDB entry 1TSR) is shown as a cartoon of different shades of colours from the N- (yellow) to the C-terminus (green). (B) Substates of the S6-S7 loop in the 2D subspace described by the first and second PCs. The different substates are indicated by capital letters (A-D). The average 3D structure identified for each substate is represented as a cartoon and the conformations of the S6-S7 loop are highlighted in red, orange, yellow, green and cyan for A, B, C and D, respectively.

MATERIALS AND METHODS

Preparation of molecular dynamics (MD) simulations

We used the X-ray structure of the p53 DBD (PDB entry 1TSR chain B, residues 95–289, in complex with DNA (PDB entry 1TSR chains E and F) as the initial structure for classical MD p53DBD-DNA simulations (49). We employed the same structure (after removing the DNA) as starting structure for classical MD p53DBD simulations. We also performed additional simulations (p53DBD-DNA(91-289) and p53DBD(91-289)) using both MD and PT-metaD starting from the X-ray structure of the p53 DBD (50) with four additional residues at the N-terminus (PDB entry 2XWR chain A, residues 91–289). The atomic coordinates of the DNA were modelled from 1TSR upon structural alignment of the two p53 DBDs. Further details are provided in Supplementary Text S1. MD simulations were performed with Gromacs 4.5 (51) coupled to the PLUMED 1.3 plugin for PT-metaD (52). The systems were solvated in a dodecahedral box of Tip3p water molecules (53) at 150 mM NaCl, applying periodic boundary conditions with a minimum distance of 1.3 nm between protein and box edges. We used the CHARMM22/CMAP protein force field (54) and CHARMM27 parameters for DNA (55). In both systems the interactions between zinc (Zn2+) and its coordinating residues were treated using a non-bonded model with van der Waals and electrostatic terms available in the CHARMM22/CMAP force field. To further confirm our results, we also performed PT-metaD simulations using Zn2+ parameters developed by Merz et al. (56) combined with the CHARMM22/CMAP force field for p53DBD-DNA(91–289), as well as simulations with the protein CHARMM22* force field (57) for p53DBD-DNA(91–289) and p53DBD(91–289).

Classical MD simulations

We collected a total of 25 independent classical MD simulations (replicates) of 50–100 ns each (Supplementary Figure S1) at 300 K and 1 atm using the isothermal-isobaric ensemble (NPT) and an external Berendsen bath with thermal and pressure coupling of 0.1 and 1 ps, respectively. The list of non-bonded interactions was updated every 10 steps and conformations stored every 4 ps. Other details on the preparatory steps and the parameters used for MD and PT-metaD are provided in Supplementary Text S1.

Analyses of MD trajectories with TimeScapes

We used the TimeScapes package version 1.3 (58) to identify conformational changes associated with DNA interactions with the p53 DBD during the MD simulations. The TimeScapes algorithm identifies conformational changes through analysis of the formation and breaking of atomic contacts in the system under investigation. In particular, we employed two different metrics available in Timescapes, i.e. (i) a coarse-grained model to calculate contacts between residue side chains and to measure system activity that is defined as the total number of changes in contact per unit time and (ii) a finite difference approximation of the first derivative of the total Generalized Masked Delaunay (GMD) tetrahedralization-based activity to define local minima (basins) and local maxima (transitions) (58). We used the same residue subset used for the principal component analysis (PCA) (see Supplementary Text S1), considering all protein atoms. The half width median filter was set to 2 ns. To separate contacts/non-contacts between the side chains of each pair of residues we set a buffer region with edge order of 2 (contacts) and 3 (buffer) for the GMD graphs, as suggested in the original publication (58) to avoid re-crossing events in the MD simulations.

Enhanced sampling MD simulations

We used PT-metaD (23) in the well-tempered ensemble (59) to enhance sampling. In these simulations, sampling of the free energy surface is enhanced by adding a history-dependent potential to a set of collective variables (CVs). We also further enhanced the sampling by allowing exchanges between different temperatures through a replica exchange approach. We employed eight replicas (at 296K, 298K, 300K, 308K, 320K, 332K, 345K, 358K) where the width of the energy distribution (of all but the ‘neutral’ 298K replica) was increased as previously described (24). All replicas are further subject to an additional biasing force through metaD simulations in which a Gaussian of width 0.1 nm in all the four dimensions is deposited in the collective variable space every 4 ps (initial height of 0.12 kcal/mol and bias factor of 6). The set of CVs are Cα–Cα distances that are described in Supplementary Text S1 and Supplementary Figure S2. PT-metaD simulations were carried out for 0.3 μs per replica (aggregate time 2.4 μs). Convergence of the simulations is discussed in Supplementary Figure S2B.

Paths of communication

The shortest paths of communication between Lys120 and residues of S6-S7 loop were calculated from the unbiased MD simulations using the Protein Structure Network/Linear Mutual Information (PSN/LMI) method (60). Briefly, the Cα LMI matrix was calculated by averaging over non-overlapping windows of 1 ns, and using a cut-off of 0.4 to reduce noise. The Icrit value employed in our simulations was 7, calculated as previously described (61). The PSN was calculated for each structure and only edges present in at least half of the simulation frames were considered. The PSN, LMI and PSN-LMI calculations were performed using WORDOM (62). We mapped the paths on the 3D structures by the xPyder (63) plugin for PyMOL.

Protein Interactions by Structural Matching (PRISM)

PRISM is a motif-based software to predict protein-protein interactions (64). It includes a rigid-body structural comparison of target proteins to known templates of protein–protein interfaces and a flexible refinement using a docking energy function. We employed PRISM (64,65) to screen 54 different potential interaction partners for the p53 DBD. We used 40 conformations of p53 DBD from our simulations, including both the DNA-bound and -unbound states, along with the average structures from each PCA substate and they were subjected to energy minimization with the steepest descent algorithm. The potential interaction partners selected in the first round step were then re-submitted to PRISM2.0 and evaluated on the basis of their docking energy scores considering an energy value cutoff of −2.39 kcal/mol for favourable interactions, as previously suggested by the developers of the method (64). To benchmark our PRISM analyses, we performed calculations using four crystallographic structures of p53 DBD obtained in presence or absence of DNA and known p53 partners for which experimental data of the complexes with p53 DBD are already available (see Supplementary Text S1 for details on the selected structures). We also calculated the Interface Similarity score (IS-score), which is often used to assess models of protein-protein complexes (66,67). IS-score was calculated for those models for which a reference experimental structure of the complex was available. This score evaluates geometric and side-chain contact similarity of the predicted interface with respect to a reference complex interface and it varies in a range from 0 to 1 for divergent and perfect models, respectively.

RESULTS AND DISCUSSION

Prediction of long-range coupling between the DNA-binding L1 loop and the S6-S7 loop

To identify the regions of p53DBD that are influenced by the interaction with DNA, we first performed unbiased MD simulations and a joint PCA of p53DBD and p53DBD-DNA, inspired by a previously published work (43). Comparing the subspace described by the two first principal components (Supplementary Figure S3), we observe (as expected) noticeable differences in the dynamics of the L1 loop that is in direct contact with DNA (Supplementary Figure S3). In the simulations, in the absence of DNA, we find that the L1 loop also samples conformations known as ‘recessed states’ (Supplementary Figure S4, Figure 2A and B). These states are involved in the early events of DNA recognition (43) as experimental studies have shown that the L1 loop populates two different conformations, one extended and elongated toward the DNA and one ‘recessed’ (Supplementary Figure S4) (42).

Figure 2.

(A) Average structure of each conformational substate of the principal component analysis (PCA) subspace. The average conformations from A-D are shown as red, orange, yellow and green cartoons, respectively. The DNA-interacting residues (Lys120, Arg248, Arg273, Arg280) are shown as sticks. We here show the coupling between the different structural states of S6-S7 loop and conformations of L1 loop ranging from states not competent for DNA interaction as A (red) and D (green) and DNA-bound-like states as B (orange) and C (yellow). The DNA is taken from the initial structure (PDB entry 1TSR, chain B) and shown as a reference in grey. (B) Conformations of loop L1 in extended (pink) or recessed (blue) state in the X-ray structure of p53 tetramer (PDB entry 3Q06 monomer B and C). The L1 loop and the S6-S7 loop of a conformation from a generalized masked delaunay (GMD) minimum are shown as red cartoons for reference. (C–L) Activity level during p53DBD-DNA (C–H) and p53DBD (I–L) simulations. We used an approach that employs a contact metric based on higher-order statistics (GMD) to describe conformational changes in the p53 DBD, measured as activity (see Materials and Methods). The GMD activity plot is reported for each simulation. We identified basins of minimal activity (GMD minima) of the trajectories, and the corresponding structures are reported as cartoons in Supplementary Figure S6. We showed the conformations of the S6-S7 and L1 loops in each GMD minimum compared to the states observed by (A) PCA and the recessed and extended states observed in X-ray structure of (B) p53 tetramer (PDB entry 3Q06), respectively. Structural states of the S6-S7 loop in GMD minima are represented by spheres coloured according to the structural similarity with average conformations of substates identified by PCA (A) (A red, B orange, C yellow, D green) measured by Cα RMSD of the S6-S7 loop. Structural states of the L1 loop in GMD minima are represented by rectangular boxes coloured according to their structural similarity with extended (pink) or recessed (blue) state in the X-ray structure of p53 tetramer (B) (PDB entry 3Q06, monomers B and C, respectively) measured by Cα RMSD of L1 loop.

More intriguingly, we also observed differences in the dynamics of the S6-S7 loop when comparing simulations of p53DBD and p53DBD-DNA (Supplementary Figure S3A). This loop connects the S6 and S7 β-strands (Supplementary Figure S3A) and is located 3 nm away from the DNA-binding site. We identified four distinct conformational substates for this loop in the PCA subspace, labelled A-D (Supplementary Figure S3B, Figure 2A). The starting structure for the simulations has the S6-S7 loop in an A-like conformation and L1 in an extended conformation. Simulations of the DNA-bound protein suggest that interaction with DNA promotes B- and C-like structures (occluded states for the S6-S7 loop), in which the S6-S7 loop gets closer to the flanking S4 β-strand (Supplementary Figure S4B). Similar observations on structural heterogeneity of the S6-S7 loop were also pointed out in a recent simulation study of a p53 mutant variant (45).

To further examine the potential coupling between conformational changes at the DNA-binding interface and those in the S6-S7 loop, we selected a set of key residues for DNA interaction (Lys120 in L1, Arg248 in L3, Arg273 in loop S10-H2 and Arg280 in H2 α-helix) and analysed their conformations in the different PCA regions (Supplementary Figure S5). Lys120, Cys277 and Arg280 are residues known to be important for the interaction with the DNA major groove, whereas Arg248, Ala276, Arg273 and Ser241 mostly mediate the interaction with the minor groove. Among those, we focused on Lys120, Arg248, Arg273 and Arg280 since their substitutions are associated with severe abrogation of p53 DNA binding (68,69). Moreover, Arg248 and Arg273 are also altered in many human tumours (68,69).

Although our definition of the four substates A–D was solely based on the conformation of the S6-S7 loop, we observed that the different conformations of the S6-S7 loop also correlate with distinct conformational preferences in L1, which ranges from substates not optimal for DNA interaction (in A and D conformations) and DNA-bound-like substates (in C conformations) (Figure 2A and B). In p53DBD B-like states, we observed both Lys120 conformations. In particular, in D and A, the side chain of Lys120 is displaced outward by 1 nm compared with the states sampled in p53DBD-DNA simulations (Figure 2A, Supplementary Figure S5). These structures also resemble experimentally observed ‘recessed states’ (42) (Figure 2B, Supplementary Figure S4B). Recently, Verma et al. also showed by MD simulations of the p53 DBD that the L1 loop populates recessed and extended conformations in the absence of the DNA and that these states are likely to be associated with motions in the S6-S7 loop (43). In our simulations, Lys120 is characterized by conformations optimal for DNA interaction in structures from state C (Figure 2A, Supplementary Figure S5). In contrast, arginine residues at 248, 273 and 280 positions are less perturbed by DNA interaction and decoupled from the S6-S7 conformations (Figure 2A, Supplementary Figure S5), in agreement with previous simulations (40).

For an additional characterization of the structural transitions occurring in p53 simulations, especially considering the limitations associated to PCA analyses of conventional MD simulations (5,58,70–72), we also employed a method based on the so-called ‘higher-order statistics’ using the TimeScapes software (58). This analysis is based on monitoring the time-dependent formation and breaking of side-chain contacts in MD trajectories to pinpoint conformational rearrangements. We performed the analysis on each of the p53DBD and p53DBD-DNA simulations, measuring the GMD-based ‘system activity’ considered as the total number of contact changes over the unit of time. From this analysis, we identified basins of minimal activity (GMD minima) of the trajectory, corresponding to stable states of the protein. We then analysed in details the structures of each minima identified by GMD analyses with attention to the conformations of the S6-S7 and L1 loops (Figure 2C–L and Supplementary Figure S6). We calculated the Cα RMSD of the S6-S7 loop between structures of each GMD minimum and the average conformations of the different substates identified by the PCA analysis (A, B, C, D). Moreover, we calculated the Cα RMSD of the L1 loop among the structures in each GMD minimum and the monomers B and C of the X-ray structure (42) of the p53 tetramer (PDB ID 3Q06) that we can use as references of both the extended and the recessed states of the L1 loop. We then classified each GMD minimum according to its RMSD to the aforementioned L1 loop conformations.

The TimeScapes analysis confirmed that the S6-S7 loop and the L1 loop have different conformational propensities in p53DBD and simulations (Figure 2C–H and I–L for p53DBD-DNA and for p53DBD simulations, respectively). Indeed, the GMD graph of p53DBD simulations showed that the events of high activity are associated with rearrangements in the L1 loop, promoting the loss of extended states in favour of recessed-like states (Figure 2I–L), along with changes in the L3, L2 and S7-S8 loops. Furthermore, in GMD minima of p53DBD simulations, the loop S6-S7 did not undergo large rearrangements and did not significantly populate C-like conformations. On the contrary, high-activity events were not generally associated to changes in the L1 loop in the p53DBD-DNA simulations (Figure 2C–H). They are, however, often associated with conformational changes in the loop S6-S7 (Figure 2C–H), and C-like DNA-bound conformations could be observed. Moreover, each GMD minimum characterized by the S6-S7 loop in a C-like DNA-bound conformations featured the L1 loop always in extended conformations (Figure 2E and H). On the contrary, A, B and D-like conformations of the S6-S7 loop were associated both with extended and recessed conformations of the L1 loop (Figure 2I–L). These results further suggest the long-range coupling between conformations of the L1 DNA-binding loop and the S6-S7 distal loop that was observed by PCA analyses (Figure 2A, Supplementary Figure S3B).

Both X-ray crystal structures (50) and NMR relaxation dispersion experiments (39) of a longer p53DBD(91–289) variant are suggestive of a role of the N-terminal tail in modulating the conformational ensemble of the p53 DBD (Figure 3A). Thus, we performed an additional set of p53DBD(91–289) classical MD simulations to examine whether the inclusion of the N-terminal residues influences the dynamics of the S6-S7 loop. In agreement with previous experimental data (39,50,73), our simulations show that the N-terminal residues can interact with the p53 DBD through a network of hydrogen bonds, including residues in the proximity of the S6-S7 loop (such as Thr211 and Arg213) (Figure 3A). However, we do not expect classical unbiased MD simulations lasting a few hundreds of ns to sample the longer timescale dynamics revealed by the NMR experiments (39), as we also previously observed for another protein (28). Indeed, our classical MD simulations of p53DBD(91–289) do not show substantial differences between the DNA-bound and -unbound forms of the S6-S7 loop (Supplementary Figure S7). This result suggests that the N-terminal residues modulate the dynamics of the S6-S7 loop since in its absence MD simulations of the p53DBD were characterized by a larger amplitude motions of the loop on the ns-timescale. The introduction of the N-terminal residues in the study of the p53 DBD dynamical properties is thus crucially important for modelling and simulation purposes and it is here explored for the first time.

Figure 3.

(A) Sub-networks of hydrogen bonds in MD p53 DBD(91–289) simulations. We here show the cluster of hydrogen bonds that involve Thr211 and Arg213 in the S6-S7 loop, Ser94 and Ser90 of the N-terminal tail, along with Thr170 and Glu171 in the L2 loop. The Cα atoms of the residues involved in hydrogen bonds between the S6-S7 loop, the N-terminal tail and other regions are shown as spheres. The Cα atoms of the interacting residues are connected by sticks, whose thickness is proportional to the persistence of the interaction. The analysis was carried out by PyInteraph (96). (B and C) Paths of long-range communication from the L1 to S6-S7 loop. The communication paths are shown as orange (occurrence probability > 0.25) or yellow (occurrence probability < 0.25) sticks. The terminal nodes of the paths are highlighted as blue and red spheres centred on the Cα atoms for the S6-S7 and L1 loops, respectively. The paths calculated from both (B) DNA-unbound simulations and (C) DNA-bound simulations are shown. S6-S7 loop is highlighted in blue. The paths are listed in Supplementary Table S1.

Since p53 exists in solution as a tetramer (39,74) it is important to verify that the conformations of the S6-S7 loop observed in our study are still relevant in the oligomeric assembly. Thus, we compared the structures sampled by our simulations with the known experimental quaternary structures of p53 (Supplementary Figure S8). The different S6-S7 conformations nicely fit in the tetrameric assembly without structural clashes. The loop residues, apart from Phe212, are generally solvent exposed and available for interaction with biological partners.

Lastly, we employed methods inspired by graph theory to analyse our MD simulations and to provide residue-level insight into the long-range coupling between the L1 DNA-binding loop and the S6-S7 loop. Even in the unbound protein, we found evidence of structural communication between the L1 and the S6-S7 loops in both A-like and C-like structures (Figure 3B and C, Supplementary Table S1). The interaction with DNA strengthened the coupling between the two regions and also promoted new communication paths that were not observed in the free protein. In more detail, the communication between the two loops propagated from the DNA-binding site across the external face of the core β-sheet and through the L3 loop and the N-terminal tail in the p53DBD-DNA ensemble.

DNA-binding alters the free energy landscape of p53 DBD and promotes a population shift in the S6-S7 loop

The analyses collected above with the p53DBD(91–289) construct and the available experimental NMR studies (39) suggests that the conformational changes in the regions in proximity of the N-terminal tail and the S6-S7 loop are occurring on long timescales that would be difficult to sample accurately by standard MD simulations.

The data collected above allowed uniquely to identify the S6-S7 loop as potentially influenced long-range by the DNA binding but they cannot be used for an accurate estimate of the populations of the different p53 sub-states in presence or absence of DNA. Thus, to overcome the intrinsic limitations in the conformational sampling of classical MD and to obtain quantitative information on the population of the A-like/D-like states versus C-like/B-like states and the effect of DNA binding on the free energy landscape, we turned to enhanced-sampling techniques. In particular, we employed PT-metaD with a set of four CVs that we chose based on the distances between key residues of the S6-S7 loop (Asp208 and Phe212) and its surroundings (Glu258, Glu221, Arg156 and Arg158). Figure 4 reports the 2D free energy surface (FES) for one pair of collective variables as an example. The full set of 2D FES plots are reported in Supplementary Figure S2A.

Figure 4.

Interaction with DNA alters the free energy landscape of p53 DBD. The 2D free energy surface (FES) profiles are shown for two out of the four collective variables employed in PT-metaD simulations. p53DBD(91–289) and p53DBD–DNA(91–289) are shown on the left and right panels, respectively. B and C-like occluded states (the minima corresponding to lower distances in the plots) are only a minor population in the free state, whereas their population increases by more than four fold upon DNA binding. The conformation of the loop and the disordered N-terminal tail in each basin of the FES are shown in red, orange, yellow and green, for A, B, C and D states, respectively. We show only the structure of the p53 DBD for sake of clarity. The average structure for each p53 state from the unbiased MD simulations is shown as a reference in white cartoons. The other 2D FES profiles are reported in Supplementary Figure S2B.

The PT-metaD simulations clearly show that the 2D FES of p53 DBD is affected by DNA-interaction (Figure 4, Supplementary Figure S2). In particular, we observed a shift in the populations of the S6-S7 loop toward C-like states with a concomitant decrease in the populations of A-like and D-like states from more than 95% (p53DBD(91–289)) to less than 50–60% (p53DBD–DNA(91–289)). When the domain is free in solution, A and D are the dominant conformations, whereas B states are mostly absent and C states constitute a minor population of the ensemble characterized by a higher free energy with respect to the major state (ΔG ≈ 2 kcal/mol, Supplementary Figure S2B). DNA binding decreases the free-energy barriers between the different substates of the S6-S7 loop and at the same time shifts the populations so that the B/C occluded states are roughly equally populated compared to the A/D states (Figure 4, Supplementary Figure S2). The results are also robust with respect to different force-fields, as illustrated by the p53DBD–DNA(91–289) and p53DBD(91–289) metaD simulations carried out with the CHARMM22* force field (Figure 5A) and also p53DBD–DNA(91–289) with recently published protein–Zn2+ parameters (56) that we here combined with the CHARMM22/CMAP force field (Supplementary Figure S2B). There are some differences in the description of A-like and D-like structures changing the force field, such as for example with CHARMM22* they are not separated by high energy barriers (Figure 5A) resulting in populating one single minimum including either A-like and D-like structures. Another difference is that with the new Zn2+ parameters (56) we did not observe D states and A states are split in two minima, i.e. A’ and A’′ (Supplementary Figure S2B). Nevertheless, the differences induced upon DNA binding are still observable with an increase of the population of the occluded states of the S6-S7 loop.

Figure 5.

Minor and major states are consistently described by different MD protein force fields and phosphorylation of Ser215 traps the p53 DBD in the major conformation of S6-S7 loop. The 2D FES profiles are shown for two out of the four collective variables employed in PT-metaD simulations. p53DBD(91–289) and p53DBD–DNA(91–289) with CHARMM22* force field are shown on the upper left and right panels, respectively. p53DBD(91–289) with phosphorylation at Ser215 is shown in the bottom panel. A different protein force field (CHARMM22*) provides results consistent with a conformational shift upon DNA binding in loop S6-S7 towards more occluded states of the loop. Phosphorylation at Ser215 traps the protein in its major state A, also observed in the X-ray available structures of the protein.

Interestingly, a well-known phosphorylation site at Ser215, which has been suggested to inhibit transcription and destabilize the protein (46–48,75), is capable of trapping p53 in its A/D-like states in the absence of DNA. This decreases the propensity for the minor population of the S6-S7 loop when the protein in in the DNA-unbound state (Figure 5B), suggesting an extra level of regulation at the conformational level.

In addition, we calculated Cα and Cβ NMR chemical shifts and compared them to previously published values (38). In particular, we predicted the average chemical shifts from an ensemble of conformations that included all the structures in the minima of the 2D FES reported in Figure 4 for p53DBD–DNA(91–289) using PPM (76). We then also predicted chemical shifts using uniquely the experimental X-ray structures that we used as starting structures for our simulations (Figure 6). Notably, we found a substantially improved agreement with the experimental data when the metaD-derived ensemble is used (which includes either the major and minor states of the S6-S7 loop) with a χ2red of ∼2.5 compared to χ2red of ∼4.5–8 when the two experimental structures were used. A better agreement with the experimental data was also locally observed for the S6-S7 loop, where for most of the residues Cα (Figure 6A) and Cβ (Figure 6B) chemical shifts featured lower χ2 values. The comparison of our simulated ensemble with the NMR data support the view that the conformation of the S6-S7 loop is heterogeneous and that it can exist in multiple substates in solution.

Figure 6.

Comparison of calculated and experimental Cα and Cβ chemical shifts. We predicted Cα and Cβ chemical shifts using PPM (76) and compared them to published NMR chemical shifts of p53DBD (38) calculating a χ2 for each atom, accounting also for the errors associated to the prediction of that specific class of atoms as reported in (76). In the figure, the results for the residues in the S6-S7 loop and the residue Leu206 in position -1 with respect to the loop are reported using as reference the metaD ensemble and the two different experimental structures used as starting structure of p53 DBD simulations. We can observe a drop in χ2 for most of the residues when the metaD ensemble is used.

The S6-S7 loop is close to the disordered N-terminal region that is known to undergo a slow conformational exchange in the free protein (39), whereas it is not clear whether DNA binding modulates this dynamics. The exchange between these two experimentally observed states has previously been related to the disruption of the interactions between Trp91 and Arg174 (39). Indeed, the authors showed that the p53 DBD is characterized by a major state (<90% of the ensemble), in which the disordered tail tightly interacts with the DBD, and a minor population where these interactions are lost. Our results suggest that changes in the hydrogen bonds mediated by the residues in S6-S7 (Figure 3A) could contribute to the process. According to our simulations, the loss of all these interactions not only causes the opening of the p53-hinge region located N-terminal to the DNA-binding domain (Figure 4), but also promotes the conformational changes in the S6-S7 loop. Our finding that the DNA alters the free energy landscape of p53 DBD in the S6-S7 loop suggests that the slow dynamics observed by NMR (39) might be regulated via DNA binding.

The S6-S7 loop is a recruitment site for biological partners

To correlate the conformational changes observed in the S6-S7 loop with functional properties, we wondered whether the different states of the S6-S7 loop could act as hotspots for protein–protein interactions. We employed the PRISM method to model protein–protein interactions (64). We used a subset of 40 conformations from our MD simulations of p53 DBD (95–289), both in the A-like (DNA-bound) and C-like (DNA-unbound) forms, along with the average structures from each PCA substate (Figure 1B) and available experimental structures of p53 DBD (PDB entries: 1TSR chain B (49), 2XWR chain A (50), 2GEQ chain B (77), 1YCS chain A (78)). We selected 54 different proteins with a known 3D structure that play a role in p53-mediated pathways (79) or have been demonstrated to physically interact with the full-length p53 (80) or have been deposited in the Molecular Interaction database (81) (Supplementary Table S2). We used these proteins as interaction partners to be screened for their capability of binding to p53 by the PRISM method. We evaluated the PRISM results on the basis of the docking energy scores for the predicted complexes, using a cutoff value of energy score of −2.39 kcal/mol, as previously suggested (64).

At first, we benchmarked PRISM by evaluating its capability to predict known protein–protein interfaces in the complexes of p53DBD and its partners for which experimental structures are already available (38,78,82–84). These complexes included the apoptosis-stimulating p53 protein 2 (53BP2) (PDB entry 1YCS), RelA-Associated inhibitor (iASPP) (PDB entry 2VGE), Oncoprotein SV40 Large T-antigen (PDB entry 2H1L), 53BP1 (PDB entry 1GZH) and BCL-xL (PDB entry 2MEJ). We used the IS-score (see Materials and Methods) to compare the predicted interface with the experimental reference (66). An IS-score of one would be associated with a model that perfectly matches the experimentally known interface. IS-scores in our comparisons were generally higher than 0.6 with the exception of the complex with SV40 Large T-antigen, which was poorly predicted and with an IS-score lower than 0.2. In these complexes, the docking scores were in the range of −13 to −5 kcal/mol, confirming the validity of the cutoff proposed by the PRISM developers (64). Nevertheless, we should notice that for those proteins for which a template for the interface is available in the PRISM database (such as 53BP2) the IS-scores are higher than 0.8, whereas other known complexes with p53 but currently absent in the PRISM database generally provide lower IS-scores.

Once validated, we continued to apply PRISM to other proteins that have been reported to interact with p53. Our PRISM modelling of complexes with p53DBD or p53DBD–DNA conformations suggests several interaction partners for p53 DBD (53BP2, iASPP, Chk1, Casp3, Cdk7, Crk, Cdk2, Plk1, E2f1, Ku70, Ark1, Supplementary Table S2, Figure S9), some of which have also been previously identified (36).

In addition to binding partners that we predicted would bind to both classes of conformations, we also found proteins that appeared to prefer to interact specifically with either A-like or C-like conformations. Some of their interaction interfaces on p53 are predicted in the region of the DBD that comprises the S6-S7, S9-S10 loops and the β-sheet composed of the S4, S9 and S10 β-strands (Supplementary Table S2, Figure S9, Figure 7). This observation suggests that different conformations of the S6-S7 loop affect the spectrum of cofactors for p53 DBD and its conformation thus plays a role in modulating the recruitment of binding partners. In agreement with this hypothesis, changes in the S6-S7 loop from A-like to C-like conformations alter locally the solvent accessible surface, along with the electrostatic surface potential (Supplementary Figure S10). The C-like states, which are the ones more populated in p53DBD–DNA, feature less accessible residues in the S6-S7 loop (such as Asp208 and Thr211) and in the S4 β-strand (such as Arg158) compared to the A-like structures (Supplementary Figure S10).

Figure 7.

Complex predicted by PRISM for (A and B) Ku70 and (E and F) Ark-1 p53 DBD. The complex predicted for interaction of Ku70 (PDB entry 1JEY, (97)) and Ark-1 (PDB entry 1MUO (98)) and p53 DBD is depicted as an example of the predicted interfaces reported in Supplementary Table S2. Moreover the crystallographic structure of the complex between (C and D) Nb139 and p53 DBD is shown for comparison. Structures are shown as a cartoon. We calculated the residues of the two proteins in each complex that have at least one atom within 5 Å of distance from the atoms of the binding partner and we highlighted their Cα atoms with coloured spheres. p53 DBD is highlighted in blue, Ku70 in yellow, Ark-1 in cyan and Nb139 in green. The other interactors identified by PRISM and modulated by different S6-S7 conformations are described in details in the Supplementary Table S2.

Among p53 partners that are likely to interact exclusively with A/D-like states of p53 DBD, we identified protein kinases, some of which are involved in p53 regulation, especially Cdk7 and the Aurora-A kinase (Ark1, Figure 7E and F). The results obtained for other kinases are described in detail in Supplementary Table S2. Interestingly, we showed by metaD simulations that phosphorylation occurring at Ser215, which is one of the target site for Ark1 (46,47,85), traps p53 DBD in its major A-like state (Figure 5B). DNA interaction, on the contrary, disfavour the sampling of p53 A/D-like open states (Figure 4) and, according to our PRISM analyses, also protects the DBD by the interaction with Ark1, a kinase that can inhibit p53 transcriptional activity and promote p53 unfolding (46,47).

Interestingly, PRISM also predicts Ku70 as a cofactor that interacts with p53DBD structures only when the S6-S7 loop is in the A-like conformations with energy docking scores lower than −2.39 kcal/mol (Supplementary Table S2 and Figure 7A and B). Ku70 is involved in control of apoptosis and in transcription-independent functions of p53, mainly through formation of an inhibitory complex with Bax, which is a pro-apoptotic factor (86). The Ku70-Bax complex is cytoprotective and impairs initiation of apoptosis (86). In addition, it has been shown that acetylated p53 in complex with Ku70 can promote the dissociation of Bax from Ku70 (87,88), contributing to transcription-independent induction of apoptosis (89). The interaction between Ku70 and p53 is, however, p53-acetylation independent (87). In support of our prediction, mass spectrometry and pull down assays have pointed out that p53 interacts with Ku70 in vitro (87). In our predicted complex, the p53 DBD interacts with the N-terminal domain of Ku70 while experimental data show that Bax interacts with the C-terminal region of Ku70 (88), making thus possible the simultaneous binding of the two proteins with Ku70. Our predictions suggest that the S6-S7 loop and surrounding residues recruit the N-terminal domain of Ku70 and this is possible only when the S6-S7 loop is in A/D-like states, which are disfavoured by DNA-binding according to our simulations. We also showed by simulations that the interaction with DNA shifts the population of the S6-S7 loop towards occluded conformations that protect from the interaction with Ku70, thus down-regulating this transcription-independent apoptotic mechanism.

Recently a novel nanobody, the monoclonal single-chain antibody Nb139, was developed to specifically interact with the wild-type p53 DBD, without altering its structure and function (90). Although Nb139 does not affect p53 DNA binding activity, it alters p53 transcriptional activities and perturbs the transactivation of p53 target genes. The crystallographic structure of the complex of p53DBD and Nb139 (Figure 7C and D) shows that the interaction interface of Nb139 on p53DBD includes many van der Waals interactions and hydrogen bonds in the β-sheet region of p53, which is flanked by both the S6-S7 loop and the N-terminal disordered tail (90). The occluded conformations of the S6-S7 loop promoted by DNA interaction in our simulations reduce the accessibility of a region that overlaps not only with Ku70 or Ark1 binding interfaces but also where Nb139 is known to bind (90).

In particular, using a 5 Å distance cutoff in the protein–protein interfaces of the aforementioned p53-mediated complexes (p53-Ku70 and p53-Nb139), we found p53 residues such as Ser96, Arg156, Glu204, Tyr205, Asp207, Asp208, Arg209, Phe212, Thr256, Asp259, Ser260, Ser261 for complexes with ku70 and Tyr103, Pro152, Met160, Leu206, Asp208, Arg209, Arg158, Asn210, Thr211, Ile254, Glu258, Asp259, Ser260, Ser261, Asn263, Arg267 for complexes with Nb139.

Altogether our simulations and PRISM analyses suggest a novel interaction interface for cofactor recruitment on p53 DBD, which could be tested experimentally by mutagenesis. For example, candidates for experimental mutagenesis are either residues of the S6-S7 loop (such as Asp207, Asp208, Arg209 and Thr211), in the S4 β-strand (such as Arg158 and Arg156) and in the S9 β-strand, such as (Thr256 and Glu258). All these residues feature changes in their accessibility upon conformational changes of the S6-S7 loop and are also involved in the interfaces predicted for Ku70, Nb139 and Ark-1.

In summary, our results overall suggest that the S6-S7 loop of p53 DBD is modulated by the interaction with DNA and it is a regulatory site for p53 transcription-independent functions, possibly involved in recruitment of binding partners and further modulated by post-transcriptional modifications.

CONCLUSIONS

The available X-ray structures of the p53 DBD show very similar conformations for both DNA-bound and -unbound states (91). Nevertheless, high flexibility of p53 DBD has been suggested as crucial for its biological functions (29). p53 structural plasticity can also be exploited to design inhibitory small molecules (92). Conformational changes may occur not only at the interface for DNA binding (40,41,43) but also at distal sites, such as the interface between the p53 DBD core and the disordered N-terminal region (39,74,93) or the C-terminal domain (44). Until now, many MD studies of p53 have focused on conformational changes in the immediate proximity of the DNA binding site (40–43), instead of changes that might potentially occur at distal sites. p53 is central in cellular signalling and interacts with many different partners to activate and regulate diverse activities and biological pathways (29,30). A hub protein such as p53 can either exploit different binding interfaces or can often use the same binding region to interact with multiple partners (32–34).

Here, we propose that the S6-S7 loop in the p53 DBD populates different states in solution. In our simulations of p53DBD(91–289), we were only able to sample these at equilibrium using enhanced sampling techniques, suggesting that they occur on longer timescales than accessible by our unbiased simulations. Exposed or occluded states of the loop S6-S7 were also observed in recent classical MD simulations of a p53 mutant variant with a different force field than the one here employed (45). The enhanced sampling techniques that we used makes it possible, in principle, to describe conformational changes on the micro-millisecond time scale (5,25,28), such as the ones experimentally observed for the regions surrounding the N-terminal disordered tail of p53 DBD (39), where also the S6-S7 loop is located. We found that DNA binding modulates the populations of different conformations of the S6-S7 loop. The structural communication between the two distal sites occurs through the L1 DNA-binding loop and the N-terminal disordered region. We also showed that post-translational modifications, such as the phosphorylation of Ser215, could modulate the conformational properties of the S6-S7 loop and disfavour the occluded states for transcription-dependent activity, suggesting a complex network of regulation of p53 function at the S6-S7 loop. We note that mutations of residues in the S6-S7 loop are also cancer-related (94,95). Those mutations are known to affect p53 transcriptional regulation (94,95), suggesting that the DNA-regulated conformational changes that we observed in the S6-S7 loop have an even broader functional relevance which might be compromised by cancer-related mutations.

We found that DNA binding promotes states in which residues of the S6-S7 loop lose most of the interactions with the disordered N-terminal region and get involved in interactions with residues in the S4 and S9 β-strands, as well as in the L2 loop. In those conformations, certain residues of the S6-S7 loop and S4 β-strand are shielded from the solvent, altering the accessibility of an interface region that affect the recruitment of biological cofactors, such as Ku70, Ark1 and Nb139. We predict that they interact more strongly when p53 DBD is not bound to DNA. We thus speculate that the population shift of p53DBD towards conformations less favourable for Ku70 binding could affect the signalling mediated by p53. Ku70 interaction with p53 is indeed necessary to release and activate Bax, playing a role in the multiple steps involved in initiation of the apoptotic pathways (89). The important role of the p53 DBD in apoptotic regulatory functions is also emphasized by the fact that BCL-xL, in the cytoplasm, interacts on the same surface of the DBD devoted to DNA-binding (38).

In conclusion, our study points out a new potential regulatory mechanism of p53 based on distal effects induced upon DNA binding. DNA binding promotes an interface conformation that is not compatible with binding to certain interaction partners, thus protecting p53 from interaction that mediate transcription-independent functions, such as transcription-independent apoptosis. As long as p53 is DNA-bound, the DBD is protected from interactions that mediate transcription-independent pathways and is able to carry out its transcription-dependent functions. If experimentally tested and validated, our findings will open up new avenues for mechanistic exploration towards a better understanding of the complex regulatory network mediated by p53.