Protein stability of p53 targets determines their temporal expression dynamics in response to p53 pulsing

The cellular response to DNA double strand breaks (DSBs) involves coordinated expression of several distinct pathways with unique temporal dynamics. The pathways enact cellular programs necessary for DNA repair and for determination of cell fate (Khanna and Jackson, 2001). As an immediate response, repair complexes are rapidly recruited to the sites of breakage (Nakamura et al., 2010; Polo and Jackson, 2011). During the repair process, cells enter a transient state of growth arrest. Depending on the outcome of the repair, cells will either reenter the cell cycle, permanently arrest via senescence, or undergo cell death (Noda et al., 2012).

The tumor suppressor p53 is a transcription factor that acts as a critical regulator of the many processes involved in the DSB response. In response to DNA damage, p53 is rapidly stabilized via phosphorylation of specific residues and alters expression of downstream target genes involved in DNA repair, cell cycle arrest, and apoptosis (Aylon and Oren, 2007; Paek et al., 2016). Depending on the extent and duration of the stress, p53 determines whether cells undergo transient cell cycle arrest, senescence, or apoptosis. Despite extensive study, how p53 determines cell fate and temporally coordinates the proper expression of its targets remains poorly understood.

In recent years, single-cell analyses of p53 expression dynamics have demonstrated that the temporal dynamics of p53 accumulation play a role in regulating the proper response to DNA damage. In response to DSBs induced through gamma irradiation or the radiomimetic drug neocarzinostatin (NCS), p53 undergoes undamped oscillations of expression with a relatively fixed pulse amplitude, duration, and frequency (Lahav et al., 2004; Geva-Zatorsky et al., 2006). p53 dynamics have been shown to be crucial for p53 function, in both transcriptional regulation and cell fate determination. Oscillatory p53 expression in response to DSBs has been shown to diversify the expression of downstream targets into a spectrum of mRNA expression patterns, which can be categorized by two extremes: (1) “pulsing” genes that have oscillatory mRNA expression, and (2) “rising” genes that have mRNA levels that steadily rise in accumulation over time (Porter et al., 2016). As the binding of p53 to target promoters is largely uniform across different targets (Hafner et al., 2017), the dynamics of target gene expression are determined by the stability of the target mRNA (Porter et al., 2016; Hafner et al., 2017). Changing p53 dynamics from repeated pulses to sustained expression through cotreatment with the MDM2 antagonist Nutlin-3 resulted in alteration of several target gene expression dynamics and a shift of cell fate from transient cell cycle arrest to senescence (Purvis et al., 2012), demonstrating a significant role for p53 dynamics in the regulation of cell fate.

While p53 oscillations generate variety in target gene mRNA expression dynamics, how the oscillatory dynamics impact target protein expression, and by extension cell stress responses, remains poorly understood. mRNA and protein correlations are relatively weak in most biological systems (de Sousa Abreu et al., 2009; Vogel and Marcotte, 2012), suggesting that the dynamic mRNA expression of p53 targets may not sufficiently explain the p53-mediated DNA damage response. Gene ontology has not been found to be a predictor of target mRNA expression dynamics, with distinct genes associated with cell cycle arrest and apoptosis having oscillatory or rising expression dynamics (Porter et al., 2016). Computational models have attempted to explain how pulsatile p53 expression can regulate cell fate decisions (Zhang et al., 2007, 2009; Sun et al., 2009). These models primarily rely on inclusion of p53 pulse-counting mechanisms to measure the duration of the DSB damage and drive cell death following prolonged damage. Such pulse counters have been proposed to involve post-translational modification of p53 (Zhang et al., 2007, 2009) or expression of specific p53 targets (Sun et al., 2009). Understanding the dynamics of downstream p53 target protein expression would likely improve the insights that can be obtained from such models and provide potential clarification to the pulse-counting mechanisms operating in the DSB response.

We examined p53 target protein expression dynamics to identify general principles governing the temporal regulation of the cellular DSB response. We found that while many genes with oscillatory expression dynamics yielded pulsing protein expression in response to DSBs, the relationship is not universal—some target genes with pulsing mRNA expression dynamics had nonpulsing protein expression dynamics. Using an ordinary differential equation (ODE) model, we identified three general classes of dynamics predicted by the mRNA and protein stabilities of p53 targets: (1) oscillating, (2) pulse counting, and (3) rising. We found that targets with oscillatory mRNA expression yielded oscillatory protein expression for proteins with short half-lives and rising protein expression for proteins with long half-lives. Temporal differences between the dynamic expression classes suggest a potential partitioning of p53 function into early and late responders based on the predicted accumulation rates of the target proteins. We further demonstrate that the dynamics of downstream p53 targets play a critical role in the DSB response, as minor shifts in the stability of MDM2 significantly reshape the p53 response and alterations to the stability of p53 up-regulated modulator of apoptosis (PUMA) significantly alter the fraction of cells that enter cell death pathways. These findings provide a potential mechanism for the temporal separation of p53-mediated cell fate decisions.

From our previous findings, we determined that p53’s oscillatory dynamics generate variety in p53 target gene mRNA expression dynamics, ranging from strong oscillations to continuously rising expression (Porter et al., 2016). As proteins are major functional effectors of the p53-mediated DNA damage response, we sought to determine how the mRNA dynamics impact target protein expression dynamics. In response to p53 oscillations induced by the radiomimetic drug NCS, we measured by Western blot analysis the protein expression dynamics for the p53 targets MDM2 and BBC3/PUMA, which we previously identified as having strongly oscillating mRNA expression dynamics (Porter et al., 2016). Despite exhibiting similar mRNA dynamics, these proteins showed qualitatively distinct protein expression dynamics (Fig. 1 and Fig. S1). At the population level, MDM2 expression showed damped pulses, correlating with the dynamics of both p53 protein and MDM2 mRNA. In contrast, PUMA protein levels increased in expression over time and did not show the damped pulses apparent in PUMA mRNA expression. Interestingly, the dynamics of PUMA protein were similar to the protein expression dynamics of the p53 target PIG3, previously identified as a target with “rising” mRNA expression dynamics (Porter et al., 2016). These results suggest that an additional filter operates downstream from p53 between target mRNA and protein expression to either maintain oscillatory dynamics or convert them to rising dynamics.

Based on the discrepancy between mRNA and protein expression dynamics, we sought to identify potential key parameters that dictate the filtering of p53 dynamical information at the level of target protein expression. To gain insight, we constructed an ODE model to capture the general features of the p53 regulatory system and generate predictions for which kinetic parameters affect target protein expression dynamics to the greatest extent (Fig. 2 A and Materials and methods). Our model of a generic p53 target contained four kinetic parameters governing (1) transcription rate, (2) translation rate, (3) mRNA stability, and (4) protein stability. We modeled p53 oscillations with a periodic function of period 5.5 h based on previously measured p53 dynamics (Lahav et al., 2004; Geva-Zatorsky et al., 2006). For simplicity, we assumed that the target gene is induced solely by p53. While more complicated models can better recapitulate dynamics for the induction of specific targets, our goal was to use this relatively simple model to capture general features of target mRNA and protein expression from an oscillating transcription factor that may be broadly applicable to a range of p53 targets. Changes to either transcription or translation rate had little effect on normalized dynamics and only altered absolute levels of molecules (Fig. S2); therefore, the production rates of both mRNA and protein were kept constant for all models to focus on factors that generate qualitative changes to protein expression dynamics. Thus, we focused strictly on mRNA and protein degradation rates as key parameters.

We ran simulations to predict the mRNA and protein expression dynamics for a generic p53 target with a range of possible mRNA and protein half-lives (Fig. 2 B). We found that mRNA decay rates affected mRNA expression patterns, as has been previously observed experimentally (Porter et al., 2016; Hafner et al., 2017). Specifically, target genes with fast decay rates had pulsatile expression dynamics; target genes with slow decay rates had dynamics in which mRNA levels rose continuously. We also found that genes with a short mRNA half-life were predicted to have a greater variety in qualitatively distinct protein expression dynamics. For targets with quickly decaying mRNAs, if the protein decay rate is also fast relative to the p53 pulse frequency, protein expression would have oscillatory dynamics. If instead the protein decay rate is slow, protein dynamics would increase monotonically during oscillatory p53 expression. In contrast, we found that target genes with slow mRNA decay rates relative to the p53 pulse frequency, and thus continuously rising mRNA expression dynamics, cannot undergo oscillatory expression at the protein level. Thus, at the level of mRNA and protein expression, the decay rate of the molecular species acts as a filter for oscillatory dynamics. The model also predicted that the filtering is irreversible: Once oscillatory dynamics are filtered, they cannot be regained, as no changes in protein stability could yield a pulsing protein from rising mRNA expression dynamics. Interestingly, the model also suggests significant differences in the temporal induction of p53 targets based on protein stability. For a given mRNA stability, shifting the protein stability results in increased time to half of the induction maximum (Fig. 2 C). As a result, genes that have identical mRNA dynamics may have substantially different protein dynamics (Fig. 2 B). In these conditions, many combinations of mRNA and protein stability also yield identical dynamics (Fig. 2 B) offering multiple mechanisms of regulation to generate similar protein expression dynamics.

Our model predicted many potential dynamics for p53 targets based on their mRNA and protein stabilities. To examine the predicted dynamics for known p53 targets, we combined our previously determined mRNA half-lives (Porter et al., 2016) with a literature-based identification of protein half-lives (Mauxion et al., 2008; Yen et al., 2008; Cambridge et al., 2011; Schwanhäusser et al., 2011; Boisvert et al., 2012; Craxton et al., 2012). Preference for half-life measurements was given to high-throughput studies that measured half-lives of multiple proteins under the same biological conditions. Additional preference was given to measurements obtained using endogenous protein levels, as opposed to fluorescently labeled proteins that may not necessarily have the same degradation kinetics as endogenous proteins. Based on these parameters, we composed a dataset of 36 well-characterized p53 targets. These targets included 17 targets previously shown to have rising mRNA expression dynamics, 18 targets with pulsing mRNA expression dynamics, and 1 unclassified target that had more complex mRNA expression dynamics (Porter et al., 2016).

Categorizing p53 targets by the previously determined mRNA decay rates and the protein decay rates resulted in four separate classifications based on the relationship between the target mRNA decay rate, the target protein decay rate, and the p53 pulse frequency (Fig. 3). We found that p53 targets fell into each of four quadrants, though quadrants I and IV (corresponding to mRNA and protein decay rates that were coherently fast or slow) contained more targets than quadrants II and III (corresponding to slow mRNA decay rates and fast protein decay rates, or vice versa). For this complete panel of genes, we observed a modest correlation between mRNA and protein stability (Pearson’s linear correlation coefficient = 0.559) as determined from our previous measurements of mRNA decay rates and protein decay rates from the literature (Mauxion et al., 2008; Cambridge et al., 2011; Schwanhäusser et al., 2011; Boisvert et al., 2012; Craxton et al., 2012; Porter et al., 2016). The correlation of this gene set was greater than previously determined based on correlations across the entire transcriptome and proteome (Schwanhäusser et al., 2011), for which little to no correlation between mRNA and protein stability was observed. Based on our model, each of these quadrants was expected to exhibit specific protein dynamics. Quadrant I, which contained rapidly degrading mRNA and proteins, was predicted to yield pulsatile protein expression. This category contained several well-characterized p53 targets, including MDM2 (supporting our results shown in Fig. 1), GDF15, p21, and WIP1. Both quadrant II and quadrant III were predicted to yield pulse counter dynamics, in which target protein expression rises with each p53 pulse. Interestingly, the pro-apoptotic p53 target PUMA was in quadrant II, consistent with our measurements of PUMA protein dynamics (Fig. 1). The magnitude and the number of pulses for proteins in quadrants II and III were predicted to be a function of the difference between the protein decay rate and the p53 pulse frequency. As the protein decay rate (or mRNA stability) increases, the relative change in expression with each pulse decreases, but the number of pulses that can be integrated before saturation occurs increases (Fig. S3). Quadrant IV contained p53 targets with stable mRNA and stable protein expression, and these targets were predicted to have smoothly rising or stable expression over time.

Previous studies have found that mRNA and protein half-lives define distinct protein classes based on ontology and protein function (Belle et al., 2006; Schwanhäusser et al., 2011; Sandoval et al., 2013). To determine whether mRNA and protein stabilities of p53 targets correlate with specific cellular functions, we performed a gene ontology enrichment analysis using PantherDB (Thomas et al., 2003) to identify enriched biological processes within each quadrant (Tables S1, S2, and S3). As each quadrant contained a limited number of genes and all are known to be regulated by p53, we expected that each quadrant would be enriched for genes involved in the DNA damage response and cell cycle checkpoints. As expected, we found many such pathways enriched in multiple quadrants. However, we found that different quadrants showed different enrichment for some biological processes, with quadrant I showing the highest fold enrichment for genes involved in cell cycle arrest and quadrant II showing enrichment for genes involved in apoptotic cell death. These results suggest that different protein dynamics may segregate proteins into distinct functional classes.

Based on the dynamics predicted from the ODE model and comparisons of known mRNA and protein stability, we sought to validate the predictions experimentally for our larger set of p53 targets. We quantified protein expression dynamics in lysates of cell populations by Western blot analysis at discrete time points following NCS-induced DSB formation in MCF7 cells (Fig. 4 A). Although population-level analysis can mask the true dynamics of protein expression in individual cells, it enabled us to examine a greater number of p53 targets and to monitor endogenous protein dynamics. Based on previous studies (Lahav et al., 2004; Geva-Zatorsky et al., 2006), due to desynchronization among individual cells in the population, we expected p53 targets to exhibit damped pulses in population-level assays if they track p53 oscillations. As an additional result of desynchronization, distinction between the pulse counting and rising dynamics predicted in Fig. 3 is not possible, and we make only the distinction between pulsing or non-pulsing protein expression.

We focused on protein expression of a panel of ten p53 target genes (Fig. 4 B). The panel consisted of the following subsets of genes: four genes (MDM2, GDF15, CDKN1A/p21, and PPM1D/WIP1) that have pulsing mRNA expression (Porter et al., 2016) and were predicted to also have pulsing protein expression (Fig. 3); three genes (EGFR, PIG3, and BAX) that have rising mRNA expression (Porter et al., 2016) and were predicted to also have rising protein expression (Fig. 3); and three genes (AMPKB1, TIGAR, and BBC3/PUMA) that have pulsing mRNA expression (Porter et al., 2016) but were predicted to switch to rising protein expression (Fig. 3). For the panel of p53 targets, we quantified by Western blot analysis high temporal resolution protein expression dynamics in cells treated with NCS. NCS treatment induced p53 oscillations that appeared damped at the population level, as expected (Fig. 4, B and C). Both MDM2 and GDF15 showed strong damped oscillations tracking p53 dynamics (Fig. 4, B and C), also as predicted (Fig. 3). In comparison, p21 and WIP1 showed weaker oscillatory dynamics, consistent with the fact that the decay rates of these proteins are closer to the p53 pulse frequency based on values in the literature. AMPKB1, TIGAR, and PUMA all exhibited nonpulsatile, rising protein expression despite being translated from mRNAs with pulsing expression dynamics (Fig. 3). PIG3 and BAX, which have rising mRNA expression dynamics (Porter et al., 2016), yielded stable or slowly rising protein levels (Fig. 4, B and C). EGFR showed an initial increase in expression that was not consistent with the predicted dynamics; however, the temporal induction of EGFR preceded p53 induction (Fig. 4 B comparing the first three hours), suggesting that the increase in EGFR in response to DSBs is at least in part independent of p53-mediated transcription.

Recent studies have highlighted the diversity of p53 dynamics across different cell lines (Chen et al., 2013; Stewart-Ornstein and Lahav, 2017). To examine how variations in p53 dynamics influence target protein dynamics, we performed Western blots in three additional cell lines: RPE-hTERT, U2-OS, and A549 cells (Fig. S4). With pulsatile p53 expression, both MDM2 and p21 tended to show pulses; however, as p53 shifted toward monotonic expression, targets tended to exhibit rising expression, consistent with our model.

To confirm that the proteins with pulsatile dynamics were coordinated with p53 expression, we performed a Fourier transform to estimate the dominant frequency of expression for each protein (Fig. 4 D). We measured a dominant frequency of 0.17 h⁻¹ for p53, consistent with the previously measured period of 5.5 h (Lahav et al., 2004; Geva-Zatorsky et al., 2006). GDF15 and MDM2 both exhibited a dominant frequency identical to the p53 frequency (Fig. 4 D). WIP1 and p21, which exhibit weak pulsatile dynamics, also showed frequencies at or near the p53 frequency (Fig. 4 D). As predicted, none of the proteins with rising expression dynamics showed a dominant frequency in the Fourier transform (Fig. 4 D).

We next sought to verify our prediction that the observed protein expression dynamics correlated with protein decay rates, with oscillatory expression dynamics requiring rapid protein decay. Specifically, for the panel of p53 targets, we expected that the protein decay rates for the pulsing targets should have decay rates that were faster than the p53 pulse frequency of ∼0.17 h⁻¹. To validate this prediction, we directly measured target protein decay rates under conditions of DNA damage for the panel of well-characterized p53 targets. We induced DNA damage with NCS. Then, after 3 h, we treated cells with cycloheximide to block new protein synthesis. We quantified protein levels by Western blot analysis (Fig. 5 A) and determined protein decay rates (Fig. 5 B) and by extension the corresponding protein half-lives for each p53 target (Table S4). Due to rapid degradation, MDM2, WIP1, and GDF15 required higher temporal sampling than the other proteins to determine their decay rates more accurately (Fig. 5 B). Of note, the p53 target TIGAR was relatively stable during the time-scale of the first several p53 pulses (Fig. 5, B and C), and therefore we considered its degradation rate as negligible within the initial response to DSBs. As predicted, the four proteins that showed pulsatile expression dynamics had decay rates greater than the p53 pulse frequency, and the six proteins that showed rising protein expression dynamics had decay rates less than the p53 pulse frequency (Fig. 5 C).

Taken together, our computational model and experimental data suggest that protein stability acts as a downstream filter of p53 dynamics, enabling the conversion of a target gene with pulsing mRNA expression dynamics to encode protein expression dynamics that are either pulsing or rising, depending on the relative decay rate of the protein. This additional filtering enables further tuning of p53 target gene expression dynamics beyond that for mRNA expression dynamics (Porter et al., 2016).

We next sought to confirm that p53 target protein stability impacts protein expression dynamics and, by extension, target protein function during the DNA damage response. We sought to alter target stability in two ways: (1) increasing the stability of a target with pulsatile expression dynamics, and (2) decreasing the stability of a target with rising expression dynamics. We first focused on the canonical p53 target MDM2, a well-characterized target with pulsatile dynamics. Previous studies identified the S395 residue of MDM2 as a key site controlling MDM2 stability, with phosphorylation of this site by the ATM serine/threonine kinase (ATM) leading to elevated levels of MDM2 degradation during the DSB response (Fig. 6 A; Maya et al., 2001; Stommel and Wahl, 2004; Batchelor et al., 2011). We performed site-directed mutagenesis to convert S395 of MDM2 to an alanine, which resulted in an increase of MDM2 stability during the DSB response (Fig. 6 B).

In the response to DSBs, MDM2 functions to shape p53 dynamics. Since MDM2 is a target of p53, we hypothesized that the initial p53 response to DSBs would be unaltered, but the later response following induction of MDM2 would depend on MDM2 stability and dynamics. To quantify the effect on p53 dynamics, we transfected plasmids expressing either WT MDM2 or MDM2 S395A driven from the MDM2 p53 response element into MCF7 cells expressing a fluorescent p53-Venus fusion protein. We challenged cells with DSBs and tracked p53-Venus expression by live-cell microscopy for 24 h (Fig. 6 C and Fig. S5 A). We measured two characteristics of p53 dynamics: pulse frequency and pulse amplitude. We observed a significant decrease in the number of p53 pulses within 24 h upon stabilization of MDM2, indicating a decrease in pulse frequency (Fig. 6 D). Most cells expressing WT MDM2 exhibited four to five pulses of p53 expression over 24 h; in comparison, the majority of cells expressing MDM2 S395A exhibited three or fewer pulses (Fig. 6 E).

To test for effects on p53 pulse amplitude, we quantified p53-Venus expression in single cells at the peaks of the first and second pulses (Fig. S5 B). We found that cells expressing WT MDM2 showed a similar range of p53 expression between the first and second pulses, though median expression was significantly decreased at the second pulse. Expression of MDM2 S395A, however, was associated with a substantially reduced range of p53 expression at the second pulse, and an even larger discrepancy in median expression between the first and second pulses, indicating a dampening of p53 induction over time. MDM2 S395A–expressing cells also exhibited a substantial decrease in the number of cells exhibiting consecutive first and second pulses (Fig. S5 C), although we focused our analysis of p53 amplitude on cells that exhibited both pulses. Taken together, our data show that stabilizing MDM2 has the functional impact of attenuating the p53 pulsatile response to DSBs.

As a second approach, we sought to determine the functional impact of converting the protein expression dynamics of a target from rising to pulsatile. Previous studies have suggested that alterations in p53 dynamics play a crucial role in the induction of apoptosis (Chen et al., 2013; Paek et al., 2016). We focused on the canonical p53 target PUMA, a well-characterized protein that functions to induce apoptosis in response to cellular stress. PUMA has pulsatile mRNA expression but rising protein expression dynamics (Fig. 1), suggesting that destabilizing PUMA protein should enable its protein expression to shift from rising to pulsatile. To destabilize PUMA, we added a C-terminal proline, glutamic acid, serine, and threonine (PEST) sequence from ornithine decarboxylase known to increase degradation of fusion proteins (Murakami et al., 1992). In comparison to an analogous WT PUMA construct, we predicted that PUMA-PEST dynamics would shift from a rising dynamic toward a pulsing dynamic and potentially impact rates of apoptosis (Fig. 7 A).

To assess how shifts in PUMA dynamics affect cell death, we used the RPE-hTERT immortalized retinal pigmented epithelial cell line, as these cells retain functional apoptotic pathways, unlike caspase-3–deficient MCF7 cells. These cells have also been shown to exhibit pulsatile p53 dynamics (Loewer et al., 2013; Chen et al., 2016; Fig. S4), which we confirmed using single-cell imaging of RPE-hTERT cells expressing a p53-Venus fusion protein (Fig. 7 B). To confirm that WT PUMA and PUMA-PEST had different protein expression dynamics, we performed Western blot analysis of RPE-hTERT cells transfected with each expression construct in response to DSBs. As predicted, WT PUMA showed rising protein expression dynamics, while PUMA-PEST showed more pulsatile expression dynamics (Fig. 7 C).

To determine whether PUMA expression dynamics affected cell fate in response to DSBs, we measured the percentage of cells that underwent cell death in response to NCS. To avoid confounding effects from endogenous PUMA expression, we knocked down endogenous PUMA using siRNA targeting its 3′ UTR. We cotransfected cells with an empty control, PUMA, or PUMA-PEST expression plasmid in which the transgene lacks the endogenous 3′-UTR, rendering transgene transcripts immune to the knockdown. We found that NCS failed to induce cell death in the control RPE-hTERT cells, consistent with previous studies showing that pulsatile p53 expression tends to induce cell cycle arrest rather than apoptosis in response to DSBs (Purvis et al., 2012; Chen et al., 2013). Transfection with the WT PUMA plasmid resulted in a significant increase in cell death overall, likely due to increased basal expression of PUMA; however, induction of DSBs did not significantly increase cell death in these cells. Transfection of cells with the unstable PUMA-PEST did not significantly alter apoptosis under steady-state (Fig. 7 D), likely due to low expression levels based on the unstable nature of the protein. Interestingly, induction of the DSB stress response resulted in a significant increase in the percentage of apoptotic cells (Fig. 7 D), suggesting that the dynamics of PUMA play a significant role in determination of whether cells undergo cell cycle arrest or activate cell death pathways.

With greater use of single cell approaches, a growing number of cellular signaling molecules have been observed to undergo oscillatory or repeated pulsatile dynamics (Levine et al., 2013). Using the p53 system’s response to DNA DSBs as a model, we have identified a general principle that may be broadly applicable to other oscillatory systems. We have previously shown that p53 pulses generate diversity in the dynamics of target mRNA expression (Porter et al., 2016), and we now have shown that target expression dynamics are further tuned to generate distinct dynamics at the level of target protein expression. In examining 10 p53 targets with a range of mRNA expression dynamics (ranging from rising to strongly pulsing expression), we determined that mRNA expression dynamics are not universally predictive of protein expression dynamics. Instead, the combination of mRNA stability and protein stability for each target defines its own specific temporal pattern of expression for mRNA and protein levels. This mechanism potentially enables distinct p53 targets to cross functional thresholds in expression at different times even when they are responding to the same environmental cue and activated by the same regulator. These findings are consistent with previous studies that have demonstrated mRNA and protein stability are critical factors regulating the dynamics of downstream targets in other signaling systems (Learn et al., 2000; Belle et al., 2006; Hao and Baltimore, 2009; Elkon et al., 2010).

Using an ODE modeling approach, we predicted that protein expression driven by a pulsatile regulator exhibits a broad spectrum of qualitatively distinct dynamic expression patterns. For simplicity, we characterized this spectrum into three major modes: (1) pulsing, (2) pulse counting, and (3) stable/rising expression. These dynamics have significant implications for regulation of p53 dynamics and p53’s function as a tumor suppressor, especially for those targets with rapidly decaying mRNA and therefore pulsatile mRNA expression dynamics that can potentially give rise to a range of protein expression dynamics. In this study, we focused on the canonical p53 targets MDM2 and PUMA, which are known to have important functions for the DNA damage response. We demonstrated that modest shifts in the stability of MDM2 had functional consequences in the shaping of p53 dynamics by decreasing the pulse frequency and attenuating pulse amplitude, which are likely to impact several other downstream effects. More directly, we found that destabilizing PUMA and switching its protein expression dynamics increased the responsiveness of cells to DSBs in their activation of cell death. We expect that mutations that alter the stability of other pulsing p53 targets may have a profound impact on target gene function. Characterized mutations in the p53 pulsing target MLH1 (Porter et al., 2016), a gene whose product plays a critical role in DNA repair processes and which is associated with colorectal cancer, have been shown to alter the stability of the protein product (Perera and Bapat, 2008). Based on our findings, we expect that these mutations would alter the protein dynamics of MLH1 and are likely to disrupt normal MLH1 function and DNA repair processes in cancer cells.

Previous work examining global mRNA and protein stabilities across the entire transcriptome and proteome in a human cell line suggested that broad functional classes of proteins may be segregated based on mRNA and protein decay rates (Schwanhäusser et al., 2011). From our analysis focused solely on p53 targets, gene ontology analysis suggested that targets clustered by mRNA and protein decay rates correlate with discrete functions in cell fate decisions, including cell cycle arrest and apoptosis. This correlation suggests that decay rates may be a potential mechanism for temporally segregating distinct pathways downstream of p53. For example, p53 targets with high mRNA and protein decay rates, and thus with pulsatile expression dynamics, were enriched for cell cycle arrest–related genes, including p21, BTG2, and GADD45A. The predicted dynamics from our study suggest that these proteins rapidly accumulate, thereby promoting immediate cell cycle arrest in response to DNA DSBs. However, due to the rapid decay rates of these proteins, upon repair of DNA damage and cessation of p53 pulsing, the levels of the cell cycle arrest mediators rapidly decrease, allowing cells to reenter the cell cycle. Interestingly, work examining the role of p53 dynamics in regulating cell fate found that sustaining p53 activity resulted in permanent cell cycle arrest, as opposed to apoptosis (Purvis et al., 2012). Our work suggests that high sustained levels of cell cycle arrest regulators, in addition to early activation of senescent-associated genes (Purvis et al., 2012), may be a contributing factor in the permanent cell cycle arrest.

Previous computational models examining the role of p53 pulses in cell fate determination have invoked a p53 pulse-counting mechanism to distinguish mechanisms for activating cell cycle arrest versus apoptosis (Zhang et al., 2007, 2009; Sun et al., 2009). Our model predicts the existence of several pulse-counting targets of p53. Interestingly, the group of pulse-counting targets we identified were enriched for regulators of apoptosis, including the well-characterized pro-apoptotic target PUMA. Computational models examining p53-mediated cell fate decisions have previously proposed that PUMA may function as a pulse counter (Sun et al., 2009), and our results support this hypothesis. Other studies have speculated that p53 undergoes pulse counting through accumulation of specific post-translational modifications, such as phosphorylation of serine 46 (Zhang et al., 2007, 2009). Our results suggest that temporal ordering of downstream pathway activation can arise as a property of the relationship between p53 pulse frequency, target mRNA decay rates, and target protein decay rates. While it is likely that post-translational modification of p53 can fine-tune specific gene regulatory activities, the requirement of distinct temporal steps in the activation of downstream events can occur solely by p53 pulsing.

Interestingly, we found that shifting the dynamics of PUMA from rising toward pulsing resulted in increased apoptosis in response to NCS, despite reducing the stability of a pro-apoptotic factor. Studies of mitochondrial priming have indicated that it is the relative levels of apoptotic activators and inhibitors that may be important for inducing cell death, rather than absolute thresholds of any given factor (Certo et al., 2006), and such priming has been shown to be important for p53-dependent apoptosis (Liu et al., 2013). Other studies have highlighted that bimodal switching of p53 dynamics from pulsing to monotonic rising expression can induce apoptosis, and often this depends on the rate of p53 accumulation (Chen et al., 2013; Paek et al., 2016). These changes are likely associated with more rapid accumulation of downstream p53 targets. Based on our results, less stable proteins show more rapid fluctuations in expression level (pulsing dynamics as opposed to rising), potentially promoting apoptosis by altering the balance between pro- and anti-apoptotic factors and allowing pro-apoptotic factors to accumulate faster than apoptotic inhibitors at certain times following DNA damage, in agreement with previous findings. However, further studies will be necessary to fully elucidate the mechanisms involved.

We identified a potential group of p53 targets predicted to either slowly rise or remain largely stable in expression in response to DNA damage. This group contains known metabolic targets of p53 such as AMPKB1 and TIGAR (Bensaad et al., 2006; Feng et al., 2007), for which we experimentally validated the predicted protein expression dynamics. We observed little to no change in expression of AMPKB1 and TIGAR even 24 h after DSB induction. As it would likely be beneficial for cells to avoid significant changes in metabolism in response to transient stimuli, these results make intuitive sense.

Our results highlight a general principle for how biological oscillators can either filter or retain oscillatory dynamics in downstream pathway activation, potentially generating temporal ordering of major stress responses. These findings provide insight into the mechanisms behind not only the p53-mediated response to DNA damage but also a growing number of cell-signaling regulators found to undergo oscillatory or repeated pulsatile dynamics (Tay et al., 2010; Albeck et al., 2013). Given the relative ease with which mRNA and protein decay rates can be altered through relatively few mutations, our results are likely to provide novel insight into understanding the deregulation of pathways downstream of p53 in cancer contexts. Additionally, the ability to specifically perturb mRNA and protein stability to establish a distinct temporal ordering of events may prove useful in the construction of complex outputs from synthetic circuits driven from a biological oscillator.

MCF7 breast carcinoma cells were maintained at 37°C and 5% CO₂ in a base medium of RPMI containing 10% FBS, 100 U/ml penicillin, 100 µg/ml streptomycin, and 250 ng/ml amphotericin B (Corning). MCF7 cells expressing fluorescently tagged p53 (p53-Venus; Batchelor et al., 2008) were grown in base medium supplemented with a G418 at a maintenance concentration of 400 µg/ml. RPE-hTERT cells were maintained in DMEM/F12 medium containing 10% FBS, 100 U/ml penicillin, 100 µg/ml streptomycin, and 250 ng/ml amphotericin B (Corning). RPE p53-Venus cells were generated by transfection of RPE-hTERT parental cells with the pDESTN-MTp-p53-Venus plasmid used to establish the MCF7 p53-Venus cell line (Batchelor et al., 2008). Cells were transfected using Fugene6 (Promega) at a 3:1 ratio following the manufacturer’s protocol, and stable cells were selected with 400 µg/ml G418. A clonal cell line was isolated by growth from a single-cell limiting dilution, and p53-Venus expression was confirmed by fluorescence microscopy using the microscope system described below. RPE-hTERT cells expressing p53-Venus were further supplemented with 400 µg/ml G418. U2-OS cells were maintained with McCoy’s 5a medium supplemented with 10% FBS, 100 U/ml penicillin, 100 µg/ml streptomycin, and 250 ng/ml amphotericin B. A549 cells were grown in F12K medium containing 10% FBS, 100 U/ml penicillin, 100 µg/ml streptomycin, and 250 ng/ml amphotericin B.

Plasmids expressing WT or S395A MDM2 were generated as previously described (Batchelor et al., 2011). Generation of the mCherry nuclear marker plasmid was performed by Gateway cloning (12537023; Thermo Fisher Scientific) of pDONR221 containing a nuclear localization signal, pDONRP2RP3 containing mCherry, pDONRP4P1R containing the ubiquitin C promoter, and the pDESTP destination vector. Generation of PUMA and PUMA-PEST plasmids were also performed using Gateway cloning. The coding sequence for PUMA without a stop codon was contained with a gblock (IDT) and cloned into pDONR221. This plasmid was then combined with pDONRP4P1R containing a p53 response element, and either a stop codon or the PEST sequence in pDONRP2RP3. All plasmids were confirmed by sequencing before transfections. Plasmid transfections were performed using 1 µg of MDM2 plasmid and mCherry plasmid. Identical conditions were used for the transfection of PUMA and PUMA-PEST in RPE-hTERT cells. Transfections were performed using Fugene6 (Promega) at a 3:1 ratio in six-well plates.

Cells were cultured in either 35-mm or 60-mm plates for Western blot analysis. Culture medium was removed from cells by aspiration, and cells were washed once with PBS. Cells were then scraped into 1 ml of PBS and the plate rinsed with an additional 1 ml of PBS added to the collected cells. Cells were pelleted by centrifugation for 10 min at 13,000 rpm and flash-frozen in a dry ice:EtOH bath. Cell pellets were lysed in lysis buffer (50 mM Tris, pH 7.5, 100 mM NaCl, 0.5% sodium deoxycholate, 1% Triton X-100, and 0.1% SDS) supplemented with 50 mM NaF, 1 mM PMSF, and 1:100 phosphatase inhibitor I (Sigma-Aldrich). Lysates were incubated on ice for 30 min and insoluble debris pelleted by centrifugation. Protein concentration was determined by Bradford assay. 20 µg of total lysate as run per lane and separated using 4–20% Tris-Glycine gels (Bio-Rad). Proteins were transferred to Immobilon-FL polyvinylidene fluoride membrane. To facilitate probing for multiple targets, membranes were cut into sections and individual membrane sections probed for p53 and downstream targets based on the expected molecular weight. Membrane sections were blocked for 1 h in Licor Blocking Buffer (Licor) and incubated with primary antibodies overnight. Membranes were washed three times with PBS/0.1% Tween before incubation with Licor IrDye secondary antibodies (1:5,000) for 1 h. Use of IrDye-based secondaries allows for separation of multiple targets into distinct detection channels based on host species of the antibody. All blots that were probed for multiple targets were done so sequentially to ensure adequate resolution between targets and nonspecific bands did not confound quantification. Membranes were washed three times with PBS/0.1% Tween and imaged using an Odyssey Imager. Tubulin served as the internal loading control and was probed on each membrane. Stripping of blots was avoided except for instances of proteins with similar molecular weights, as in the case of p53 and tubulin. Stripping was performed using NewBlot polyvinylidene fluoride Stripping Buffer (Licor) according to the manufacturer’s protocol. Fig. S1 contains uncropped images demonstrating blotting strategy and efficacy of stripping for p53 quantification.

For quantification, fluorescence intensity data were acquired by the Li-Cor Odyssey imaging system and then analyzed using ImageJ software (National Institutes of Health). Integrated band intensity was calculated using the ImageJ Analyze Gels tool to obtain local background-subtracted band intensity for regions of interest that are manually defined by the user on each Western blot image. Integrated band intensity was measured for each individual target protein and normalized by the intensity of the internal tubulin loading control performed on each individual membrane. Relative expression for each protein was then calculated by dividing tubulin-normalized values for each time point by the normalized band intensity for the starting point (time equal to zero) for each membrane. Final quantification of results represents the mean of multiple biological replicates along with SEM.

All primary antibodies were used at a concentration of 1 µg/ml unless otherwise indicated. The following antibodies were used in our study: mouse monoclonal anti-p53 DO-1 (1:10,000; sc-126; Santa Cruz Biotechnology), mouse monoclonal anti-MDM2 SMP14 (1:500; sc-965; Santa Cruz Biotechnology), mouse monoclonal antiGDF15 G-5 (sc-377195; Santa Cruz Biotechnology), mouse monoclonal anti-TIGAR G-2 (sc-74577; Santa Cruz Biotechnology), mouse monoclonal anti-PIG3 A-5 (sc-166664; Santa Cruz Biotechnology), mouse monoclonal anti-PUMAα/β G-3 (sc-374223; Santa Cruz Biotechnology), mouse monoclonal anti-EGFR A-10 (sc-373746; Santa Cruz Biotechnology), rabbit polyclonal anti-WIP1 H-300 (sc-20712; Santa Cruz Biotechnology), mouse monoclonal anti-BAX 6A7 (sc-23959; Santa Cruz Biotechnology), mouse monoclonal anti-AMPKB1 Z14 (sc-100357; Santa Cruz Biotechnology), and mouse monoclonal anti-tubulin β E7 (Developmental Studies Hybridoma Bank; 1:3,000). Secondary antibodies used were anti-mouse 680 RD and anti-rabbit 800 CW (Licor) and were used at a dilution of 1:5,000.

To estimate protein half-life during the DNA damage response, MCF7 cells were treated with 50 µg/ml cycloheximide and 400 ng/ml NCS. Cells were then harvested according to the Western blot protocol for each indicated time point. Protein expression was quantified by Western blot analysis and measured using ImageJ software. The expression data were normalized to expression at the zero time point before log transformation. Results were fit to a linear model with a slope representing the degradation rate α. Half-life was calculated according to the equation T_1/2 = ln(2)/α. The protein TIGAR was found to be stable over the measured time, and therefore the degradation rate was set to 0. For plotting on logarithmic scales, this value was redefined to 0.001. For p53 targets that were not experimentally measured, half-lives were identified using a literature-based approach (Mauxion et al., 2008; Yen et al., 2008; Cambridge et al., 2011; Schwanhäusser et al., 2011; Boisvert et al., 2012; Craxton et al., 2012). We preferentially selected half-lives that were measured using high-throughput systems, such as stable isotope labeling by amino acids in cell culture (SILAC), as these studies could measure multiple proteins under identical conditions. Proteins that exhibited large discrepancies (measurements suggesting both very short and very long half-lives) between multiple studies were omitted from analysis.

To measure dynamics of p53 and its targets in response to DNA DSBs, cells were treated with 400 ng/ml of the radiomimetic NCS (Sigma-Aldrich). This drug has previously been demonstrated to rapidly induce DNA DSBs (Shiloh et al., 1983). Cells were harvested at indicated time points, and lysates were used for Western blot analysis.

The ODE model of p53 target dynamics was generated using MATLAB software (MathWorks). In this model, we defined changes in mRNA and protein levels using simple equations dependent on p53 expression:

For these systems, we defined p53 expression using a periodic function whose values fluctuated between 1 and 5 with a frequency of 5.5 h, mirroring the dynamics of p53 induction observed within single cells:

Simulations were performed using MATLAB for ranges of mRNA and protein degradation rates.

For live-cell imaging, cells were plated into six-well glass-bottom plates (Mattek). For MCF7 cells expressing a p53-Venus fluorescent fusion protein (Batchelor et al., 2008), the cells were cotransfected with plasmids expressing a nuclear-localized mCherry and either WT MDM2 or MDM2 S395A expressed from the human MDM2 promoter (Batchelor et al., 2011) and imaged 48 h after transfection. The nuclear-localized mCherry served as a marker of transfection, and for subsequent analysis, we examined strictly cells that showed mCherry expression. For RPE-hTERT cells expressing p53-Venus, cells were imaged 24 h after plating. For experiments with either cell line, 1 to 2 h before imaging, medium of all samples was replaced with RPMI medium lacking riboflavin and phenol red (Gibco) and supplemented with 2% FBS, 100 U/ml penicillin G, 100 mg/ml streptomycin sulfate, and 250 ng/ml amphotericin B (Corning). Cells were imaged with a Nikon Eclipse Ti-inverted fluorescence microscope equipped with an automated stage (Prior) and a custom chamber to maintain constant 37°C temperature, high humidity, and 5% CO₂. DSBs were induced with NCS as described to induce the p53 stress response. Multiple positions were analyzed per experiment with images acquired every 20 min for 24 h using a YFP filter set (Chroma; 488–512-nm excitation filter, 520-nm dichroic beam splitter, and 532–554-nm emission filter for both cell lines analyzed) and mCherry filter set (Chroma; 540–580-nm excitation filter, 585-nm dichroic beam splitter, and 593–668-nm emission filter for MCF7 cell line only). Images were collected using a 20× CFI Plan Apochromat Lambda (NA 0.75) objective (Nikon). For each condition, at least 50 individual cells were tracked. Images were segmented manually, and average fluorescence levels for nuclear p53 were determined in an automated manner using custom MATLAB software.

Cell death was quantified in RPE-hTERT cells using the FITC–Annexin-V and propidium iodide apoptosis kit (640914; Biolegend). Briefly, endogenous PUMA was knocked down by transfection with Accell siRNA targeting the 3′ UTR of the endogenous PUMA (A-004380-17-0005; Dharmacon). Cells were cotransfected with control, PUMA, or PUMA-PEST constructs lacking the 3′ UTR to enable protein expression. After 24 h, cells were treated with 400 ng/ml NCS. 48 h after NCS treatment, cells were collected and analyzed by flow cytometry and FlowJo software (version 10). Gates were assigned using single positive as well as negative controls. All conditions were repeated three times with biological replicates. Statistical analysis was performed using MATLAB and a one-way ANOVA for multiple comparisons.

Genes identified within each quadrant were entered into PantherDB (Thomas et al., 2003) and analyzed for enrichment in biological processes. Tests were run against the Homo sapiens reference list and statistics performed using the binomial statistical test with Bonferroni correction for multiple hypothesis testing. Data were sorted according to fold enrichment.

Fig. S1 contains uncropped Western blots related to Fig. 1 demonstrating internal loading control, stripping efficacy, and distinct migration patterns for p53 targets. Fig. S2 shows the impact of shifting mRNA and protein production rates in the ODE model. Fig. S3 models the impact of protein stability on pulse-counting dynamics predicted by the ODE model. Fig. S4 shows how variation in the p53 response across additional cell lines correlates with dynamics of downstream targets. Fig. S5 shows single cell traces for all cells analyzed in Fig. 6 and additional analysis of pulsatile dynamics. Tables S1, S2, and S3 contain significantly enriched biological pathways based on gene ontology analysis of genes in Fig. 3. Table S4 contains calculated protein half-lives and decay rates from Fig. 5.

We thank members of the Batchelor laboratory, as well as D. Levens and the members of the Levens laboratory. We also thank J. Stommel and J. Liu for helpful discussions, and K. Wolcott and the National Cancer Institute Flow Cytometry Core Facility for help with flow cytometry.

This work was supported by the Intramural Research Program of the Center for Cancer Research, National Cancer Institute, National Institutes of Health (Intramural Research Project 1ZIABC011382).

The authors declare no competing financial interests.

Author contributions: R.L. Hanson and E. Batchelor conceived the study and designed experiments. R.L. Hanson performed experiments and analyzed data. J.R. Porter contributed to the design of the ODE model. R.L. Hanson and E. Batchelor wrote the paper.