Advertisement

Open Access

Source Data

Transparent Process

Combined experimental and computational analysis of DNA damage signaling reveals context‐dependent roles for Erk in apoptosis and G1/S arrest after genotoxic stress

Andrea R Tentner, Michael J Lee, Gerry J Ostheimer, Leona D Samson, Douglas A Lauffenburger, Michael B Yaffe

Author Affiliations

  1. Andrea R Tentner1,
  2. Michael J Lee1,
  3. Gerry J Ostheimer1,
  4. Leona D Samson1,
  5. Douglas A Lauffenburger1 and
  6. Michael B Yaffe*,1
  1. 1 Departments of Biology and Biological Engineering, David H Koch Institute for Integrative Cancer Research at MIT, Massachusetts Institute of Technology, Cambridge, MA, USA
  1. *Corresponding author. Departments of Biology and Biological Engineering, David H Koch Institute for Integrative Cancer Research at MIT, Massachusetts Institute of Technology, 77 Massachusetts Avenue, 76‐353, Cambridge, MA 02139, USA. Tel.: +1 617 452 2103; Fax: +1 617 452 4978; E-mail: myaffe{at}mit.edu
View Abstract

Abstract

Following DNA damage, cells display complex multi‐pathway signaling dynamics that connect cell‐cycle arrest and DNA repair in G1, S, or G2/M phase with phenotypic fate decisions made between survival, cell‐cycle re‐entry and proliferation, permanent cell‐cycle arrest, or cell death. How these phenotypic fate decisions are determined remains poorly understood, but must derive from integrating genotoxic stress signals together with inputs from the local microenvironment. To investigate this in a systematic manner, we undertook a quantitative time‐resolved cell signaling and phenotypic response study in U2OS cells receiving doxorubicin‐induced DNA damage in the presence or absence of TNFα co‐treatment; we measured key nodes in a broad set of DNA damage signal transduction pathways along with apoptotic death and cell‐cycle regulatory responses. Two relational modeling approaches were then used to identify network‐level relationships between signals and cell phenotypic events: a partial least squares regression approach and a complementary new technique which we term ‘time‐interval stepwise regression.’ Taken together, the results from these analysis methods revealed complex, cytokine‐modulated inter‐relationships among multiple signaling pathways following DNA damage, and identified an unexpected context‐dependent role for Erk in both G1/S arrest and apoptotic cell death following treatment with this commonly used clinical chemotherapeutic drug.

Synopsis

Data‐driven modeling was used to analyze the complex signaling dynamics that connect DNA repair with cell survival, cell‐cycle arrest, or apoptosis. This analysis revealed an unexpected role for Erk in G1/S arrest and apoptotic cell death following

Embedded Image

  • The cell‐fate choice between apoptosis and cell‐cycle arrest after doxorubicin treatment is dose dependent and modulated by TNF.

  • An early, transient G1/S arrest following doxorubicin treatment acts as a cell‐fate decision point between re‐entry into S‐phase and apoptotic cell death.

  • Two complementary relational modeling methods identify the proliferation and survival signal, ERK, as potentially causal for the G1/S arrest and death phenotypes.

  • Inhibitor experiments validate a role for ERK in maintenance of the early, transient arrest and in promotion of apoptosis from this arrested state following doxorubicin‐induced DNA damage.

Introduction

The DNA damage response (DDR) constitutes an evolutionarily conserved set of signaling pathways that are essential for maintaining genomic integrity in all eukaryotic cells (Hoeijmakers, 2001; Jackson and Bartek, 2009). Mutations and/or acquired defects that compromise the function of these DDR pathways results in enhanced mutagenesis, and are thought to underlie the development and progression of cancer. (Kastan and Bartek, 2004; Halazonetis et al, 2008). Paradoxically, the defective DDR of tumor cells facilitates their killing by genotoxins, rendering these types of DNA damaging agents useful for cancer chemotherapy (Kennedy and D'Andrea, 2006; Helleday et al, 2008; Litman et al, 2008; Powell and Kachnic, 2008).

A variety of chemotherapeutic agents in common clinical use, including camptothecin, etoposide, and doxorubicin, inhibit Topoisimerases I and II (Topo I/II), respectively, creating large numbers of DNA double‐strand breaks (DSBs). These lesions are thought to be particularly cytotoxic, and initiate a canonical DNA damage signaling pathway through activating the upstream PI3‐kinase‐like kinase (PIKK) ATM (Hoeijmakers, 2001; Jackson and Bartek, 2009). Recruitment of Mre11, Rad50, and Nbs1 to the DSB results in an ATM autoactivation loop that increases and sustains ATM activity, and together with the PIKKs ATR and DNA‐PK, drives the accumulation of the phosphorylated form of the histone variant H2AX (so‐called γH2AX) in the area surrounding the break (Furuta et al, 2003; Pommier, 2006; Harper and Elledge, 2007; Jackson and Bartek, 2009). γH2AX is required for the efficient retention of additional signaling molecules including MDC1, RNF8, RNF168, BRCA1, along with repair factors into foci near the break, and facilitates the activation of downstream signaling pathways that mediate phenotypic responses to damage—that is, cell‐cycle arrest, apoptosis, or repair and cell‐cycle re‐entry (Figure 1A). ATM effectors involved in DNA damage‐induced cell‐cycle arrest and induction of programmed cell death include the checkpoint kinase Chk2 and the multi‐functional transcription factor p53. Chk2, together with the kinases Chk1 and MK2, phosphorylates and functionally inactivates Cdc25 family members, thereby blocking activation of CyclinE/A‐Cdk2 and CyclinA/B‐Cdk1 complexes and establishing G1/S and G2/M checkpoints (Bartek and Lukas, 2003; Donzelli and Draetta, 2003; Reinhardt et al, 2007). p53 assists in maintenance of DNA damage‐mediated cell‐cycle arrest through upregulation of the cyclin‐dependent kinase inhibitor p21 and/or initiates programmed cell death through transactivation of pro‐apoptotic Bcl‐2 protein family members (Levine and Oren, 2009).

Figure 1.

A strategy for systematic analysis of signal transduction pathways involved in the DDR. (A) An expanded DNA damage signaling network that reports on internal cellular states and external microenvironment cues indicating the intersection and integration of multiple distinct signaling pathways that respond to or modulate chromatin integrity, cellular stress, survival, apoptosis, and cell‐cycle progression or arrest. Indicated nodes and responses in this expanded DNA damage‐response network were interrogated by quantitative western blotting for the total (W) or phospho‐forms (pW) of signaling molecules, ELISA‐based activity assay (E), or flow cytometry (F). (B) Experimental approach: cells were treated +/− doxorubicin (Dox; 0, 2, or 10 μM) and +/− TNFα (0 or 100 ng/ml) for 4 h. Four hours following treatment, the media was replaced with fresh media containing 1% FBS. Cell extracts for signal measurement were collected at 0, 0.25, 0.5, 1, 1.5, 2, 4, 8, 12, 16, and 24 h following treatment. Whole cell samples were collected and fixed for cellular‐response measurements at 6, 12, 24, and 48 h following treatment.

While cell fate after DNA damage is known to be regulated by these core DDR pathways, additional signaling pathways governing general stress and survival responses such as the PI3‐kinase pathway, the NF‐κB pathway, and the three major MAPK pathways are also likely to contribute (Figure 1A). Since extracellular signals are transduced through these pathways, they may serve as information processing junctions whereby signals from the microenvironment are integrated with the DDR. A common component of the tumor microenvironment is the inflammatory cytokine TNFα, which is present within many tumor types and appear to contribute to tumor progression (Ben‐Baruch, 2003). In addition, exogenous administration of TNFα has been used to induce tumor cell death, particularly in combination with genotoxic agents (Mocellin and Nitti, 2008). The relative importance of these additional signal transduction pathways, and the manner in which their signals are integrated together with those from the canonical DDR pathways to control the fate of DNA‐damaged cells, is poorly understood. Therapeutic manipulation of these pathways, however, could significantly alter and potentially improve the clinical response of tumors to DNA damaging therapies (Samuels and Ericson, 2006; Ahmed and Li, 2007; Reinhardt et al, 2007; Roberts and Der, 2007).

To pursue this integrative understanding in a systematic manner, we collected quantitative time‐resolved signaling data on the activation status of key nodes in several molecular signaling pathways that could potentially impact cell‐cycle control and induction of cell death, together with similar data on selected components of the core DDR pathway, in an osteosarcoma cell line treated with varying doses of doxorubicin in the presence or absence of the inflammatory tumor cytokine TNFα (Szlosarek et al, 2006). In parallel, we collected quantitative time‐resolved phenotypic data on the apoptotic and cell‐cycle regulatory response of the cell population following DNA damage. We used computational modeling approaches to identify critical context‐dependent molecular signals that modulate damage‐induced cell‐cycle arrest and the initiation of programmed cell death.

Results

An experimental treatment‐response protocol for quantifying phenotypic changes and information flow within an expanded DNA damage signaling network

Following DNA damage, cells display complex dynamic phenotypes that connect cell‐cycle arrest in G1, S, or G2/M, and DNA repair with decisions related to survival, cell‐cycle re‐entry, permanent cell‐cycle arrest, or cell death. How these cellular decisions are made remains poorly understood, but must involve the integration of large amounts of molecular signaling information distributed across pathways that reflect both the current internal state of the cell and the state of the surrounding microenvironment (Figure 1A). To explore this molecular signaling–phenotypic response relationship, we treated U2OS osteosarcoma cells with varying concentrations of doxorubicin, a clinically relevant anti‐cancer agent that generates DNA DSBs by inhibiting Topo II, in the presence or absence of TNFα, a cytokine commonly found in tumor microenvironments (Szlosarek et al, 2006). Six specific treatment conditions were examined encompassing exposure to no (0 μM), low (2 μM), or high (10 μM) doxorubicin, in the presence or absence of 100 ng/ml of TNFα, for 4 h followed by transfer to drug‐free medium without doxorubicin or TNFα. Cell lysates were prepared at 11 time points between 0 and 24 h after exposure to doxorubicin, and used in quantitative western blotting or ELISA assays to measure 17 molecular signals (Figure 1B; Supplementary Table S1). Cells were fixed at four time points between 0 and 48 h after damage for flow cytometry‐based assays of cell‐cycle progression and cell death. Together, these measurements provide a high‐dimensional quantitative time‐dependent data set that captures dynamic changes in signaling pathways and cellular responses.

The cell‐fate choice between apoptosis and cell‐cycle arrest after doxorubicin treatment is dose dependent and modulated by TNFα

In the absence of DNA damage, minimal levels of apoptosis were observed over the 48h course of the experiment, as evidenced by <5% of the cells staining positively for the apoptotic markers cleaved caspase‐3 and cleaved PARP (Figure 2A, upper panels, and 2B). The addition of TNFα to undamaged cells resulted in only a minimal increase in basal levels of apoptosis. Treatment with either 2 or 10 μM doxorubicin caused a marked dose‐ and time‐dependent increase in apoptotic cell death. The presence of TNFα dramatically increased the fraction of cells undergoing DNA damage‐induced cell death, and shortened the time to apoptosis, particularly within 12 h after exposure to doxorubicin. For example, 6 h after treatment, the amount of apoptotic cell death observed in the TNFα+2 μM doxorubicin‐treated samples was greater than threefold that observed in samples treated with 2 μM doxorubicin alone (Figure 2B, gray and black squares; 14 versus 4%, respectively), and instead was comparable to the amount of apoptosis observed in the doxorubicin alone‐treated samples after an additional 6 h of incubation (12 h time point, 17% death). These effects were even more pronounced following treatment with a higher concentration of doxorubicin. The amount of apoptotic cell death observed 6 h following treatment with TNFα+10 μM doxorubicin was nearly sixfold greater than that observed with 10 μM doxorubicin alone (Figure 2B, gray and black diamonds, 7 versus 40%, respectively) and was comparable to the amount of apoptosis observed in the 10‐μM doxorubicin‐treated samples an additional 18 h later (24 h time point, 36%).

Figure 2.

TNF enhances dose‐dependent cell death following doxorubicin‐induced DNA damage with minimal affect on dose‐dependent cell‐cycle arrest. (A) Representative scatter plots for cleaved caspase‐3 and cleaved PARP as measured simultaneously by flow cytometry in single cells, for all six treatment conditions, at 12 h following treatment. (B) Apoptosis in treated cell populations (percentage of the cell‐population staining positively for both apoptotic markers cleaved caspase‐3 and cleaved PARP), for all six treatment conditions investigated at 6, 12, 24, and 48 h following treatment. (C) Representative cell‐cycle profiles based on PI staining of DNA content for cells exposed to 0, 2, or 10 μM doxorubicin +/− TNFα at 6, 12, 24, and 48 h following treatment. (D) Quantification of the percentage of the treated cell populations in the G1 (a), S (b), G2/M (c), and pHH3+/M (d) phase of the cell‐cycle, and the sub‐G1 population (e) at 6, 12, 24, and 48 h following each treatment. (B, D) Mean values from n⩾4 independent experiments; error bars denote standard error of the mean. All measurements of cell‐cycle stage and apoptosis represent ‘window’ measurements that do not reflect the cumulative fraction of the initial population, but rather the fraction of the population that is still measurable at the indicated time point. For example, apoptosis measurements at later time points will not include cells that have died and disintegrated between the time of treatment and the time of measurement. Source data is available for this figure in the Supplementary Information.

Source Data for Figure 2B [msb20121-sup-0001-SourceData-S1.xls]

Source Data for Figure 2C [msb20121-sup-0002-SourceData-S2.xls]

Source Data for Figure 2D [msb20121-sup-0003-SourceData-S3.xls]

In contrast to the ability of TNFα to enhance DNA damage‐induced cell death, this cytokine had little effect on the cell‐cycle arrest observed in response to DNA damage by either low‐ or high‐dose doxorubicin (Figure 2C, compare left and right columns, and 2D). In untreated control cells, with or without TNFα, we observed an initial slight rise in the S and G2 population, followed by a slow progressive increase in the percentage of G1 cells, and a corresponding drop in S and G2, as the cells approached confluence. After DNA damage, dramatic cell‐cycle changes were observed in the surviving cell populations treated with low‐dose doxorubicin, whereas more subtle cell‐cycle effects were seen in the surviving cells after high‐dose treatment.

Both 2 and 10 μM doxorubicin caused a very slight increase in the percentage of G1 cells 6 h after treatment (Figure 2D). In cells exposed to low‐dose doxorubicin, this was followed by a significant accumulation of the S‐phase population at 24 h, along with a significant depletion of the G1 population (Figure 2Da and Db, black squares). By 48 h, both the G1 and S‐phase populations of these 2 μM doxorubicin‐treated cells were reduced from their levels at 24 h, and a significant increase in the population of G2‐arrested cells was observed (Figure 2C and Dc). The addition of TNFα reduced the S‐phase accumulation seen following 2 μM doxorubicin at 24 h, and reduced the percentage of G2‐arrested cells at 48 h (Figure 2D, gray squares). Importantly, this was accompanied by an increase in the sub‐G1 population at both of these time points. In cells treated with high‐dose doxorubicin, we also observed a steady and significant depletion in the G1 population from 6 to 48 h. However, this decline was muted compared with the decline observed after low‐dose doxorubicin treatment (Figure 2D, black diamonds). Similarly, the time‐dependent rise in the percentage of S‐phase cells at 24 h was also greatly reduced after 10 μM treatment, and no meaningful increase in the percentage of G2‐arrested cells was observed over the time course of the experiment.

Taken together, these observations reflect alternative cell‐fate choices in response to low versus high‐dose doxorubicin treatment. Less damaged cells preferentially trigger cell‐cycle arrest, while cells with more extensive DNA damage activate programmed cell death. This cell‐fate choice is further modulated away from cell‐cycle arrest and toward cell death by the presence of TNFα.

An early G1/S checkpoint constitutes a cell‐fate decision point

The pattern of cell‐cycle progression observed in response to DNA damage induced by doxorubicin, particularly under lower‐dose conditions, suggests that the damaged cells may undergo an early transient arrest at the G1/S boundary between 0–12 h, followed by a more or less synchronous cell fate transition of either release into S‐phase or initiation of apoptotic cell death by 24 h. To probe this process in more detail, asynchronous U2OS cells were exposed to 2 μM doxorubicin, and cell progression analyzed by flow cytometry at more frequent intervals between 0 and 42 h. As shown in Figure 3Aa–Ae, an accumulation of G1‐arrested cells was evident as early as 3 h after damage, followed by a drop in, and broadening of, the G1/S peak at 15 h and the appearance of a clear shoulder off of the G1/S peak by 18 h, indicating S‐phase entry. Under these DNA damage conditions, the greatest increase in cell death occurred between 6 and 12 h (Figure 2B), when the cells appeared to be arrested in G1 (Figure 3A). Direct flow cytometry analysis of cells stained for both DNA content and cleaved caspase‐3 (cC3) further confirmed that the bulk of cC3‐positive cells had 2 N DNA content, suggesting an apoptotic G1 population (Figure 3B). This finding was unexpected, since it is generally believed that DSB generation, as well as sensing and processing of the protein‐linked DSBs induced by Topo II poisons such as doxorubicin is strongest in S‐phase, when Topo II levels are high (Goswami et al, 1996) and replication forks are likely to encounter these impassable DNA lesions (Howard et al, 1994; Hong and Kreuzer, 2000; Bartek et al, 2004; Nitiss, 2009). In agreement with this, there is also evidence that doxorubicin‐induced cell death may preferentially occur during S‐phase (Grdina et al, 1980; Meyn et al, 1980)

Figure 3.

Cell death following DNA damage occurs from cells arrested in G1 or at the G1/S boundary. (A) Representative flow cytometry profiles showing DNA content for cells exposed to 0 (black) or 2 μM (gray) doxorubicin, at time points between 3 and 42 h following treatment. Treatment with low‐dose doxorubicin (gray traces) results in an accumulation in the G1 population (2 N DNA content) as compared with untreated cells (black traces) by as early as 3 h following treatment. This accumulation persists through at least 15 h. By 24 h following treatment, the accumulation in G1 is completely abrogated and an accumulation in the S‐phase population is apparent, suggesting synchronous release into S‐phase. (B) Representative scatter plots showing DNA content and cleaved caspase‐3 in cells treated with 2 μM Dox (left) or 2 μM Dox+TNFα (right) at 24 h following treatment. The majority of cells staining positively for cleaved caspase‐3 have 2 N DNA content (red box), consistent with apoptosis in G1 or at the G1/S boundary. (C) Cells were synchronized in late G1 and at the G1/S boundary via a double thymidine block and released directly into aphidicolin‐containing media to block progression through S‐phase. Cells were treated with or without 2 μM Dox 1 h later and analyzed 8 h later for apoptosis by flow cytometry (top). The percentage of cells staining positively for both apoptotic markers is indicated from two biological replicates in each of two independent experiments (bottom). (D) Cells treated with either 0 or 2 μM Dox were pulsed with BrdU for 4 h beginning 15 h following treatment. At 48 h following treatment, cells were analyzed by flow cytometry for DNA content and BrdU incorporation (top). Overlays of the cell cycle/BrdU distributions from cells treated with either 0 (black) or 2 μM (gray) Dox reveal an accumulation of BrdU‐positive cells with 4 N DNA content in populations treated with 2 μM doxorubicin, indicating that the population of cells synchronously released from a G1 or G1/S arrest into S between 12 and 24 h following treatment is the same population that accumulates in G2/M by 48 h following treatment. No BrdU‐positive cells re‐entered the G1 population by 48 h after 2 μM Dox, consistent with a persistent G2/M arrest. Left and right panels are identical, with the Dox‐treated cells shown either behind (left) or in front of (right) the untreated controls. (E) Cells treated with 0 or 2 μM Dox as in (D) were collected at 48 h following treatment and analyzed for DNA content and the mitotic marker, pHH3. Mitotic cells are boxed, indicating stable G2 arrest after doxorubicin treatment. (F) A model for the early G1/S decision point. Source data is available for this figure in the Supplementary Information.

Source Data for Figure 3A [msb20121-sup-0004-SourceData-S4.xls]

Source Data for Figure 3B [msb20121-sup-0005-SourceData-S5.xls]

Source Data for Figure 3C [msb20121-sup-0006-SourceData-S6.xls]

There is strong evidence that the primary mechanism of doxorubicin action is via Topo II inhibition (Nitiss, 2009). Levels of Topo II in tumors in vivo are correlated with the apoptotic response to doxorubicin, while downregulation of Topo IIa is associated with acquired resistance of tumors to this drug (Ogiso et al, 2000; Tanner et al, 2006). Importantly, however, doxorubicin is well known to cause additional types of DNA damage through Topo II‐independent effects, including free radical‐induced single‐strand breaks, direct Doxorubicin‐DNA adducts, and inter‐strand crosslinks (Swift et al, 2006 and references therein). To directly address whether doxorubicin‐induced DSBs were present in the G1 population, we treated cells with 2 μM doxorubicin and used flow cytometry to measure the percentage of cells staining positively for γH2AX as a function of cell‐cycle stage. As shown in Supplementary Figure S1, elevated levels of γH2AX could be detected in roughly 50% of the G1 cells, a value that while large, was somewhat lower than that observed in S and G2 cells where the response rate was ∼90%. Thus, while other types of doxorubicin‐induced DNA damage almost certainly exist in G1 cells, our data does indicate the presence of drug‐induced DSBs during this phase of the cell cycle. Interestingly, Ju et al (2006) have reported topoisomerase‐mediated generation of transient DSBs during normal gene transcription, providing a possible mechanism for enhanced DSB formation by Topo II inhibition even in the absence of DNA replication.

To further explore whether cell death could occur directly from G1 or at the G1/S boundary without a requirement for replicative DNA synthesis, cells were synchronized at the G1/S border using a double thymidine block, and released into media containing the DNA polymerase inhibitor aphidicolin to block progression into S‐phase. The cells were then treated with 2 μM doxorubicin 1 h later. Under these conditions, the extent of apoptosis was comparable to, or higher than that observed in the DNA‐damaged asynchronous population (Figure 3C), consistent with doxorubicin‐induced cell death occurring in G1 or at the G1/S border in the absence of continued DNA replication. As shown in Figure 3Ag–Aj, the subpopulation of cells in the asynchronous population treated with 2 μM doxorubicin that had not undergone apoptosis had largely progressed into mid‐S‐phase by 24 h, ultimately accumulating in G2 between 30 and 42 h after damage. To ensure that the increase in the G2 population at late times observed by flow cytometry represented the same population of cells that had transited through S‐phase, rather than the appearance of an artifactual G2 peak from selective death and drop‐out of the S‐phase cells at these late time points, we pulsed the asynchronous damaged cells with BrdU at 15 h, when a portion of this population is synchronously entering S‐phase, and followed the subsequent cell‐cycle progression of these BrdU‐labeled cells. As shown in Figure 3D, at 48 h following treatment, the 4 N DNA‐containing cells are strongly BrdU+, indicating that the cell population that entered S between 15 and 19 h following treatment transited through S and into G2 by 48 h. These cells stain negatively for phospho‐Histone H3 (pHH3; Figure 3E), indicative of G2 arrest. Furthermore, we found that these cells were unable to form proliferative colonies in culture, suggesting that the G2 arrest is persistent (data not shown, manuscript in preparation). Together, these data indicate that DNA‐damaged cells undergo a transient arrest in G1 or at the G1/S border. This arrest serves as a decision point from which cells either die or re‐enter the cell cycle by releasing synchronously into S and subsequently arresting in G2 (c.f. Figure 3F), perhaps as a consequence of defective G1 and/or S‐phase checkpoints in this cancer cell line.

Cells show complex patterns of signaling network activation in response to low and high amounts of DNA damage in the presence or absence of cytokine cues

To explore specific signaling pathways and molecules responsible for mediating the different phenotypic responses observed in response to low‐ and high‐dose doxorubicin, in the presence or absence of TNFα, we examined the time‐dependent changes in levels and phosphorylation states of a subset of molecules that reflect the status of the core DNA damage signaling network, along with signaling pathways that modulate stress, survival, apoptosis, and cell‐cycle progression. Western blotting and ELISA assays used for this analysis were validated to be both strictly quantitative and highly reproducible, with signaling measurements always performed in the linear regime of the assay (Janes et al, 2005) (Figure 4A; see Materials and methods; Supplementary Table S1). For several signaling molecules of interest, including Chk1 and Chk2, we were unable to develop appropriate assays due to the lack of phospho‐site specific antibodies that passed this critical validation step.

Figure 4.

DNA damage induces complex signal activation dynamics across many signaling pathways. (A) Representative quantitative immunoblotting assays, shown for phospho‐ERK and phospho‐JNK. (a, b, d, e) Assay validation and linearity. Cells were starved and serum stimulated to activate ERK, or UV irradiated (50 J/m2) to activate JNK. Lysates from these positive controls (pc) corresponding to 2.5–30 μg of total protein from these positive controls (pc), along with 10 μg of total protein from untreated negative control (nc) cells were analyzed by western blotting for active ERK and JNK using phospho‐specific antibodies. Quantified band intensity was plotted against protein loaded to verify signal linearity. Gray diamonds indicate the signal intensity obtained from 10 μg of negative control lysate, and indicate high signal to noise ratio. (c, f) Representative western blots for phosphorylated ERK and JNK from lysates (10 μg total protein/lane) at the indicated time points following each of the six treatment conditions. Positive and negative control (pc and nc) lysates (10 μg/lane) were included in all gels as internal controls for normalization. (B) Quantification of western blot data for all signaling measurements of activity under the six treatment conditions were investigated. Measured signals quantified include p‐ATM (S1981), p‐Nbs1 (S343), p‐H2AX (S139), p‐p53 (S15), p‐p38 (T180/Y182), p‐JNK (T183/Y185), p‐Akt (S473), p‐ERK1/2 (T202/Y204), and total p53 and quantification of ELISA data for activation status of NF‐κB. Data are the mean of two to three independent experimental replicates and normalized to the signal measured at 0 h (prior to treatment) to yield a measure of ‘Signal Intensity.’ Source data is available for this figure in the Supplementary Information.

Source Data for Figure 4B [msb20121-sup-0007-SourceData-S7.xls]

For the majority of signaling molecules that were examined, we observed large quantitative changes, and distinct dynamic activation profiles (Figure 4A and B; Supplementary Figures S2 and S3). MAPKs and other non‐canonical DNA damage responsive signals, including phosphorylated forms of p38MAPK, JNK, and ERK (denoted p‐p38, p‐JNK, and p‐ERK1/2, respectively) along with NF‐κB, were robustly activated across treatment conditions.

While some of the signaling responses we observed were consistent with the general understanding of the expanded DNA damage signaling network (Figure 1A), others were surprising and counterintuitive, particularly with respect to their dynamics. For example, phosphorylation of ATM in response to high‐dose doxorubicin‐induced DNA damage showed a brief early peak followed by a more prolonged late peak (Figure 4Ba). In response to low‐dose doxorubicin treatment, the early p‐ATM peak was slightly delayed and broader, while the late peak was reduced in amplitude. Paradoxically, however, two direct downstream substrates of ATM, Nbs1, and H2AX, showed greater phospho‐activation in the cells treated with low‐dose doxorubicin compared with high‐dose doxorubicin (Figure 4Bb and Bc), as well as altered kinetics, while a third ATM substrate, p53, showed more phosphorylation in the cells treated with high‐dose doxorubicin (Figure 4Bd). Phosphorylation/activation of Akt was increased by high‐dose treatment (Figure 4Bg), despite the greater amount of observed cell death under this treatment condition (Figure 2B), while NF‐κB activity was increased to a similar extent by both treatments, but with delayed kinetics in the high‐ versus low‐dose treatment group (Figure 4Bi).

The addition of TNFα caused significant changes in both the magnitudes and activation kinetics of many of the molecular signals that we measured. For example, in the presence of TNFα, the levels of ATM phosphorylation in both the early and late peaks were now higher in response to low‐dose doxorubicin treatment than to high‐dose treatment, while H2AX phosphorylation showed the opposite behavior (Figure 4Ba and Bc, compare red and green dash lines). Similarly, the presence of TNFα enhanced the phosphorylation of p53 on Ser‐15, a known ATM site, in response to low‐dose DNA damage but suppressed phosphorylation of this site by high‐dose damage, and shifted the peak activity to earlier times (Figure 4Bd). Phosphorylation/activation of p38MAPK over most of the time course was dramatically enhanced by the combination of TNFα and high‐dose doxorubicin, but was essentially unchanged by the presence of TNFα in response to low‐dose doxorubicin (Figure 4Be). NF‐κB activity was robustly and equally enhanced by TNFα in the low‐ and high‐dose doxorubicin‐treated cells, with an increased rate of rise and an earlier peak of NF‐κB activity in the high‐dose case (Figure 4Bi).

Partial least squares regression identifies two signaling axes that correlate with survival–death and cell‐cycle regulation

The divergent phenotypes and complex patterns of cell signaling that we observed in response to variable amounts of DNA damage in the presence or absence of a cytokine cue obscured a facile intuitive establishment of causal relationships between signals and responses by simple visual inspection of the data. Therefore, we applied two distinct computational approaches to better infer signal–response relationships toward ultimate causation hypotheses (Figure 5A). To correlate the data set of time‐dependent molecular signals across time and treatment conditions with the resulting cell‐cycle and cell‐death cellular phenotypes, we constructed a comprehensive relational model using partial least squares regression (PLSR), a dimensionality reduction and relational modeling technique that we and others have previously applied to similar biological data sets (Janes et al, 2004, 2005; Janes and Yaffe, 2006). PLSR reduces the original explanatory variable set (170 molecular signals encompassing ∼2000 time‐dependent measurements) by calculating an effective explanatory variable set composed of orthogonal latent variables, or principal components (PCs). PCs are linear combinations of the original independent variables that maximally capture the co‐variation within the signaling space while also simultaneously maximizing co‐variance between the signals and the responses. The first principal component (PC1) maximally captures the variance in the independent variable space. The second principal component (PC2) maximally captures the residual variance in the independent variable space not captured by PC1, and so on for any additional PCs, such that PCs are orthogonal and uncorrelated. These relatively few uncorrelated PCs serve as ‘super‐variables,’ defining a low‐dimensional ‘PC space’ in which independent and dependent variables can be plotted according to their calculated weights in each PC.

Figure 5.

PLSR maps signaling events onto cell death versus cell cycle progression axes. (A) Signaling and cellular‐response heatmaps for all six treatment conditions investigated. In the heatmap representation of the signaling data, the mean values for each indicated signal, measured at each of the 10 time points, are time ordered from top to bottom in each column. For each signal (t denotes total level, p denotes the phospho‐form, and a denotes DNA‐binding activity), the measured values across all six conditions were normalized to the maximal value ever observed under any condition, and color coded from low (blue) to high (red) activity. The heatmap of the cell response data shows mean values for the four indicated parameters of cell survival/death (% of cells that are cC3+/cPARP+; cC3+/cPARP−; cC3−/cPARP–; and cC3+ regardless of cPARP status), and three measurements of cell‐cycle state (G1, S, and G2/M), color coded by absolute percentage as indicated in the color bar at right. (B) Loadings of cellular responses at the indicated time points plotted in the PC space defined by the PLSR model. The cellular responses of apoptosis (cC3+/cPARP+) at 6 and 12 h following treatment, survival (cC3−/cPARP−) at 6 and 12 h following treatment and measures of the G1, S, and G2–M cell‐cycle populations at 24 h following treatment are plotted in the PC space based on the weights of each of the cellular responses in PC1 and PC2. Apoptotic responses measured at 6 and 12 h following treatment (red and maroon diamonds, respectively) are positively correlated with PC1, while survival responses measured at 6 and 12 h following treatment (green and dark green diamonds, respectively) are negatively correlated with PC1; PC1 can therefore be interpreted as a ‘Survival–Death’ axis. Measurements of S and G2–M cell‐cycle populations at 24 h following treatment (blue and dark blue circles, respectively) are positively correlated with PC2, while measurement of the G1 population at 24 h following treatment (purple circle) is negatively correlated with PC2, consistent with PC2 as a ‘Cell‐cycle regulation’ axis. (C) Loadings of signaling measurements (gray circles) for the indicated molecules are plotted in PC space, along with the cellular responses. Increasing symbol size denotes increasing time from 0.25 to 24 h following treatment. Source data is available for this figure in the Supplementary Information.

Source Data for Figure 5A [msb20121-sup-0008-SourceData-S8.xls]

We found that 93% of the variance in the cell‐death and cell‐cycle responses could be captured by three PCs, with the first two PCs capturing 87% of the variance and having clear biological correlates. A two‐component PLSR model often has the advantage of being more intuitive and interpretable for analysis. However, if retention of a third PC substantially enhances the ability of the model to predict the dependent variable (i.e., cellular responses in our case), then a three‐component model may be justified. To address the marginal predictive capability gained by retaining additional PCs, we evaluated the ability of one‐, two‐, and three‐component PLSR models to predict the change in apoptosis between 6 and 12 h, or the accumulation of cells in G1 and G2 at 24 h, following doxorubicin administration. As shown in Supplementary Figures S4A and S5, the cross‐validated cumulative predictive capability Q2cum of the models for these responses did not substantially increase when moving from two to three components.

While signaling from core DDR components is clearly important in controlling cell‐fate decisions following damage, we were interested in how these signals are integrated with non‐core DDR signals, such as the MAP kinases (p38, ERK, JNK), Akt and NF‐κB, that regularly monitor the cellular environment, and are known to influence life–death decisions following other types of stimuli. We hypothesized that both core and non‐core DDR signaling responses would be information rich and have predictive capability with regard to cell fate following DNA damage. To test this, we evaluated the ability of a two‐component apoptosis PLSR model built with either all of the measured signals, or built with a reduced set of only core DDR signals or non‐core DDR signals as independent variables, to predict cell death between 6 and 12 h following doxorubicin administration. As shown in Supplementary Figure S4B, a model built with only the core DDR signals as independent variables had a cross‐validated predictive capability that was as good or better than either a model built with all of the signals, or a model built with only non‐core DDR signals. However, the predictive capability of all three models was quite close (Q2=0.79, 0.78, and 0.70, respectively), and substantially higher than that of a model built using only the total protein levels of signaling molecules as the independent variable set (i.e., neglecting all phosphorylation/activity data; Q2=0.36). We interpret these results as evidence that while the core DDR signals are most informative for the prediction of apoptotic response, non‐core DDR signals carry a significant amount of complementary predictive capability on their own. We therefore focused further analysis on two‐component PLSR models built with all signals in the independent variable set in order to gain insight into the contributions from both core and non‐core signaling pathways.

In the two‐component model built to investigate relationships between signaling and cell‐death and cell‐cycle phenotypic responses to DNA damage, the measured apoptosis responses (cleaved caspase‐3 and cleaved PARP) were highly positively correlated with PC1 while the measured survival responses (uncleaved caspase‐3 and PARP) were highly negatively correlated with PC1 (Figure 5B). This clear separation along a single PC axis allows the assignment of PC1 as a ‘Survival–Death Axis,’ defining a direction in multivariate signaling space along which cells may move from life to death in response to DNA damage.

In contrast to measurements of the apoptotic response, we observed that the percentage of cells in both S and G2/M stages of the cell cycle at 24 h following treatment were highly positively correlated with PC2 while the percentage of cells in G1 at this time point were highly negatively correlated with PC2. The observation that the S and G2/M populations co‐correlate with each other along PC2, while being anti‐correlated with the G1 population is consistent with our observations and preliminary interpretation of an early, transient arrest in G1 or at the G1/S boundary, followed by release and transit through S by 24 h and into G2/M following doxorubicin treatment of U2OS cells. This clear separation of the measured cell‐cycle phenotypes along PC2 indicates that that the second PC is a ‘Cell‐cycle Arrest/Release Axis,’ defining the direction in signaling space along which cells are released from a DNA damage‐induced G1/S checkpoint to cell‐cycle progression.

PLS regression mapping of specific molecular signals confirms a role for p53 and identifies an unexpected role for ERK in promoting cell death and enforcing cell‐cycle arrest following DNA damage

The physical interpretation of the first two PCs as a survival–death and cell‐cycle arrest/release axis define a vector space in which plots of the signal measurements might reveal their contributions to each of these cellular decision‐making processes. p53, for example, has well‐established roles in promoting both cell‐death and cell‐cycle arrest following DNA damage in a manner that depends on the activation status of the ATM‐Chk2 pathway (Jiang et al, 2009; Levine and Oren, 2009). Plots of both total and ATM‐phosphorylated p53 in PC space revealed that the intermediate time point measurements were highly positively correlated with PC1 (Figure 5Cd and Cj), consistent with a causal role for ATM phosphorylation of p53 in cell death (Banin et al, 1998; Powers et al, 2004; Jiang et al, 2009). In contrast, early and late time point measurements of total p53 were highly correlated with PC2, with early p53 elevations correlating with the G1 arrest phenotype, and late p53 elevations consistent with accumulation of cells in G2/M (Figure 5Cj). Intriguingly, these correlations between p53 and G1 and G2/M cell‐cycle arrest along PC2 were less pronounced when ATM‐phosphorylated p53 was used as the independent variable (Figure 5Cd), consistent with recent data that induction of p53 in the absence of ATM phosphorylation is primarily involved in cell‐cycle control rather than in promoting apoptosis (Jiang et al, 2009).

In contrast to p53, plots of NF‐κB activity measurements in PC space showed that intermediate and late time point measurements were highly correlated with measured apoptosis, while NF‐κB activity at other times mapped close to the center of the plot, and were not correlated with any measured cellular response or with either PC (Figure 5Ci). Plots of p‐p38 in the PC space, on the other hand, revealed that at early and intermediate times, its activity was best correlated with the cell‐cycle arrest phenotype, while at late times, nearly all of the measurements were highly correlated with measured apoptosis along PC1 (Figure 5Ce). These findings are consistent with known roles for NF‐κB in life–death decisions, and the p38 pathway in both cell‐cycle arrest and cell death (Karin and Greten, 2005; Manke et al, 2005; Cuadrado and Nebreda, 2010).

Surprisingly, when p‐ERK signals were plotted in PC space (Figure 5Ch), we found that ERK activity at early and intermediate times showed a strong correlation with both the G1 arrest phenotype on PC2, and with apoptosis along PC1. ERK activity at intermediate and late time points was correlated along PC2 with the accumulation of cells in S‐phase by 24 h. This mapping for phospho‐ERK signals in PC space is strikingly similar to that observed for p53 signals, and was intriguingly paradoxical given the known roles of ERK in promoting cell‐cycle progression through G1/S in undamaged cells, and in promoting survival following general types of stress (Weber et al, 1997; Ley et al, 2003; Luciano et al, 2003; Chambard et al, 2007; Ewings et al, 2007; Meloche and Pouyssegur, 2007).

A complementary time‐interval stepwise regression technique implicates a specific role for ERK in the G1/S cell‐fate decision point following DNA damage

Because PLSR is constrained to incorporate all of the signaling and response data into a final model, large amounts of co‐variation between measurements, and ‘noise’ from non‐contributory signals, can obscure or over‐predict signal–response relationships (Janes et al, 2008). Given these possibilities, the paradoxical role for ERK in G1 cell‐cycle arrest and induction of programmed cell death suggested by PLSR was suspect. We therefore developed time‐interval stepwise regression (TI‐SWR), a complementary computational method that builds on the established variable selection method, stepwise regression with forward selection (Efroymson, 1960), to identify specific signaling events, or combinations of events, that were highly correlated (and likely causal) for specific phenotypic outcomes. The PLSR method starts with all time‐dependent signals as independent variables in a model that seeks to explain one or more cellular responses, and is inherently a dimensionality reduction method whose end‐state must still contain all of the signals but with appropriately disparate weightings. TI‐SWR starts with an empty model and sequentially builds a model by incorporating as independent variables only signals that are significantly correlated with a specific cellular response (see Supplementary Methods).

In the TI‐SWR protocol, the time‐dependent change in a specific cellular response (e.g., the change in percentage of cells in S‐phase between 12 and 24 h following DNA damage) serves as the dependent variable in a linear regression model. The set of all signals measured at a specific time point serve as candidate independent variables (e.g., all signals measured at 4 h following DNA damage; see Figure 6A; Materials and methods). The forward stepwise regression algorithm begins with an empty model and considers each signal in turn, such that all possible univariate regression models are sampled. The single most substantively correlated signal is then retained into a ‘first generation’ univariate model. Next, the effect of addition of any second signal to the initial univariate model is considered, such that all possible bi‐variate models based on the original univariate model are sampled. The next signal whose addition to the model most substantively improves its correlation with the cellular response is retained in the ‘next generation’ model. This is repeated until the addition of no other signal significantly improves the overall signal–response correlation. The final model thus contains a subset of the signals measured at a specific time point as independent variables. These signals are closely coupled in a correlative sense to the cellular response used as the dependent variable, and potentially causal drivers of it. This model‐building process is repeated for each time point that precedes the dependent variable cellular response. The end result is a set of models, one for each time point, which captures time‐dependent correlations between signals and the dependent variable cellular response. This can then be repeated for every cellular response measured in the data set, or for selected subsets of cellular responses.

Figure 6.

Time‐interval stepwise regression identifies a role for early and intermediate ERK signaling in G1/S arrest and G1/S induced cell death after DNA damage. (A) Illustration of the TI‐SWR method. TI‐SWR posits that a subset of molecular signaling events measured at a given time point are related in a causal fashion to the change in a cellular phenotypic response measured over a later time interval. In this example, the change in the S‐phase population between 12 and 24 h is the dependent variable, and the independent variable set includes all signaling measurements taken at the 4‐h time point. (B) Sparkline representation, grouped by signal, of correlations identified using TI‐SWR between time‐dependent molecular events and cell‐death responses (a–c) or cell‐cycle responses (d–f). Bar direction reflects the sign of the correlation coefficient, and bar height indicates the strength of the correlation (calculated as [(R2 × 100) × (5(regression degree−1)]), where R denotes the initial regression coefficient (primary correlate) or the improvement in regression coefficient (secondary correlates) from the prior model. Red bars indicate primary correlates (regression degree=1), blue bars indicate secondary correlates (regression degree 2 or 3). Signaling measurements at 4 h implicated in S‐phase progression between 12 and 24 h are shown in (e) beneath the asterisk. The final regression model is shown in (C).

Figure 6B shows an application of TI‐SWR to model cell survival and death between 6 and 12 h, and cell‐cycle progression between 12 and 24 h, after doxorubicin treatment. This analysis revealed a strong correlation between ERK activity at 4 and 8 h following treatment and G1 arrest at 12 h followed by release by 24 h (Figure 6Bd) as revealed by models using ΔG112−24 h as the dependent variable. A similarly strong correlation was observed between ERK activation at 4 or 8 h and S‐phase entry by 24 h in model using ΔS12−24h as the dependent variable (Figure 6Be). These findings are consistent with a relationship between ERK activity and a G1/S arrest at 12 h followed by synchronous release and S‐phase entry by 24 h following treatment. In addition, there was also a positive correlation between synchronous S‐phase entry by 24 h and ERK activity at 2 h after doxorubicin treatment in a model where the primary and secondary correlating variables were γH2AX and total H2AX levels, respectively (Figure 6Be). Together, these data indicated a conceptual interpretation whereby ERK activity beginning at 2 h after DNA damage played a causal role in enforcing an early G1 checkpoint that terminates by 24 h. In addition to the cell‐cycle phenotype, there was a strong positive correlation in asynchronous cells between cell death measured between 6 and 12 h following treatment (as quantified by ΔCC3+ cells6−12 h), and both total ERK1/2 at early times and ERK1/2 activity at 12h following doxorubicin treatment (Figure 6Bb and Bc), suggesting that ERK at these times after DNA damage may invoke a pro‐death response.

To directly test these cell cycle and apoptosis predictions, we inhibited ERK activation using the MEK inhibitor PD98059, at either the same time as doxorubicin treatment, or 2 h later. As shown in Figure 7A and B, inhibition of ERK activity at the same time as doxorubicin treatment, or 2 h following doxorubicin treatment, resulted in a premature release of these cells from G1 into S‐phase, indicating that ERK activity plays a critical role in maintenance of a DNA damage‐induced arrest in G1 or at the G1/S boundary. Figure 7B quantifies this early synchronous release into S‐phase when ERK activity is inhibited, revealing that inhibition of ERK activity does not change the size of the population released into S‐phase but results in its premature release from a damage‐induced G1/S arrest ∼12 h earlier than when ERK is active. Furthermore, inhibition of ERK activity resulted in an ∼40% decrease in doxorubicin‐induced apoptotic cell death (Figure 7C).

Figure 7.

Early and intermediate ERK signaling is involved in maintenance of a cell‐cycle arrest in G1 or G1/S and in mediating cell death after DNA damage. (A) Time‐dependent cell‐cycle profiles for cells exposed to 2 μM Dox in the absence (gray) or presence (black) of PD98059 co‐treatment. Comparison of the gray trace to the black trace reveals that addition of PD98059 at 0 h following treatment with doxorubicin results in early release of a population from G1 or G1/S arrest into S‐phase and later into G2/M, revealing a role for ERK activity in maintenance of the arrest in G1 or G1/S. (B) Percentage of S‐phase cells in untreated or 2 μM Dox‐treated asynchronous cultures in the presence or absence of PD98059 added immediately or 2 h following administration of doxorubicin or vehicle control. Data are representative of two to four independent experimental replicates. (C) Percentage of apoptotic cells 12 h following treatment with or without 2 μM Dox in the presence or absence of co‐treatment with PD98059. (D) Cells were synchronized in late G1 via a double thymidine block and released directly into aphidicolin to prevent progression through S‐phase, followed by treatment with vehicle or 2 μM Dox in the presence or absence of PD98059 1 h later. Apoptotic cell death was assayed 8 h following treatment by flow cytometry. (C, D) Mean values from n⩾4 independent replicates; error bars denote standard error of the mean and * indicates P<0.05 (Student's t‐test, two‐tailed). (E) A model for ERK1/2 in cell‐fate decisions following DNA damage. ERK1/2 sustains an early, transient arrest in G1/S following DNA damage. This arrest may serves as a decision point from which cells either undergo apoptosis, or are released as a synchronous population into S by 24 h following damage. G1/S‐directed cell death is enhanced by ERK1/2 activity, either directly (dashed arrow) or indirectly by prolonging a cell‐cycle arrest during which other death‐promoting factors are active. Source data is available for this figure in the Supplementary Information.

Source Data for Figure 7A [msb20121-sup-0009-SourceData-S9.xls]

Source Data for Figure 7B [msb20121-sup-0010-SourceData-S10.xls]

Source Data for Figure 7C [msb20121-sup-0011-SourceData-S11.xls]

Source Data for Figure 7D [msb20121-sup-0012-SourceData-S12.xls]

We had previously observed that a substantial percentage of cell death following low doxorubicin treatment could occur directly from the G1 or G1/S‐arrested state (Figure 3B). Furthermore, inhibition of the major DNA replication polymerases in these doxorubicin‐treated cells resulted in similar or greater levels of cell death as that seen in doxorubicin‐treated asynchronous cells capable of progressing into S‐phase (Figure 3C). We therefore investigated whether ERK1/2 activity was specifically involved in this process. As shown in Figure 7D, addition of PD98059 to doxorubicin‐treated G1 cells in the presence of aphidocolin (to block DNA synthesis) resulted in a significant reduction in the percentage of cells staining positively for the apoptotic markers cC3 and cPARP. These data suggest that the DNA damage‐induced apoptotic cell death that occurs out of a damage‐induced G1 or G1/S‐arrest is ERK1/2 dependent. Together, these data establish a causal role for the correlations revealed from both relational modeling approaches, and indicate specific roles for ERK in promotion of a G1 or G1/S cell‐cycle arrest and apoptotic cell death from this arrested state following DNA damage.

Discussion

In this study, we generated a quantitative signal–response data set for human U2OS osteosarcoma cells in response to doxorubicin‐induced DNA damage, and used relational modeling techniques to examine the flow of information between signaling pathways and cellular phenotypic responses. Our data set revealed a series of complex and non‐obvious dynamic cellular phenotypic and molecular signaling events whose magnitudes, kinetics, and inter‐relationships were highly dependent on both the extent of DNA damage and the surrounding cytokine microenvironment.

We observed that TNFα synergized with doxorubicin to increase cell death after DNA damage, tipping the balance between cell‐cycle arrest and apoptosis. This effect was most pronounced at lower levels of DNA damage, and has potential implications for understanding both the response of individual tumor cells to intrinsic DNA damage, and the heterogeneous response of whole tumors to exogenous DNA damaging therapies, depending on local cytokine profiles within and immediately adjacent to the tumor. There have been several reports that exogenously administered TNFα can synergize with genotoxic chemotherapy to improve tumor regression in both human patients and animal models of cancer, but the effect has been largely attributed to the effects of TNFα on tumor vasculature and on the innate and adaptive immune system (Mocellin and Nitti, 2008). Our data, however, indicate that TNFα can dramatically enhance chemotherapy‐induced cell death through modulation of signaling pathways within the tumor cells themselves. The detailed molecular mechanism underlying this enhancement of doxorubicin‐induced cell death by TNFα remains unclear, but since TNFα was found to elicit an increase in early JNK and p53 signaling and a loss of late AKT signaling, interplay between these pathways may be implicated in the observed enhancement of cell death. Consistent with this, we observed that the specific inhibition of robust late AKT activity, but not early AKT activity, following high‐dose doxorubicin resulted in a >2.5‐fold increase in apoptotic cell death (Supplementary Figure S6), indicating that late Akt signaling is a pro‐survival signal under these treatment conditions. The loss of this signal in cells co‐treated with TNFα may therefore at least partially explain the dramatically higher levels of apoptotic cell death seen in that treatment condition. Interestingly, the U2OS cell line used in this investigation contains two wild‐type copies of a functional NAD(P)H:quinine oxidoreductase 1 gene (NQO1*1) gene (Fagerholm et al, 2008). Recent findings have implicated NQO1 in regulating cell death in response to both anthracyclines (i.e., epirubicin and doxorubicin) (Fagerholm et al, 2008; Jamshidi et al, 2011), and TNFα (Ahn et al, 2006) through effects on NF‐κB, p53, and cellular anti‐oxidant activity. These findings suggest that NQO1 may be a key integrator of the cellular response to these inputs, potentially playing a role in the observed synergistic response to TNFα and doxorubicin co‐treatment.

In contrast to the early cell death observed following treatment of cells with high levels of doxorubicin, low‐dose doxorubicin treatment resulted in the accumulation and synchronization of a population of cells in late G1 or at the G1/S boundary at early times after damage. This cell‐cycle arrest in G1 or at G1/S, which could also include some cells that were damaged during M when checkpoints are relatively inactive (Giunta et al, 2010; van Vugt et al, 2010), serves as a cell‐fate decision point, with cells either dying actively by apoptosis from this arrested state, or re‐entering the cell cycle by proceeding into S‐phase in a more or less synchronous manner. While we cannot completely exclude the possibility that a small amount of DNA replication had occurred prior to initiation of apoptosis in the cells that die from this arrested state, we observed similar amounts of cell death in doxorubicin‐treated G1 or G1/S‐arrested cells when DNA replication was prevented by inhibiting the major replicative DNA polymerases with aphidocolin (Figure 3C). This result indicates that doxorubicin is capable of inducing apoptosis in G1 or G1/S‐arrested cells in the absence of bulk DNA replication, and may explain why even slowly proliferating tumors with a substantial population of G1 cells are sensitive to killing by doxorubicin (Grdina et al, 1980; Colly et al, 1984; Ling et al, 1996; Wartenberg and Acker, 1996). Importantly, the cells in our experiments were freely cycling at the time of DNA damage. Whether doxorubicin induces a similar extent of apoptosis in quiescent G0‐arrested cells, was not investigated, although some data indicate that this may occur (Ritch et al, 1982; Wartenberg and Acker, 1996). Finally, our treatment protocol in which cells were exposed to doxorubicin for 4 h followed by drug removal is likely to at least partially mimic the bolus administration of doxorubicin that is used clinically. Intriguingly in this regard, the early G1/S arrest followed by synchronous S‐phase entry that we observed following low‐dose doxorubicin treatment in our in vitro cell culture model has also been observed in a rat model study of human acute myelogenous leukemia following doxorubicin administration in vivo (Colly et al, 1984).

Some molecules that are directly activated by DNA DSBs displayed the expected response, with high‐dose doxorubicin treatment leading to greater activation than low‐dose treatment. Other molecules, however, displayed a reversal of this relationship over some or the entire portion of the time course. Unexpectedly complex dynamics were observed even for signaling molecules in the same ‘linear’ pathway. For example, the kinetics of ATM activation/autophosphorylation on Ser‐1981 was closely correlated with ATM‐mediated phosphorylation of Nbs1 on Ser‐343 after low‐dose DNA damage, but not after high‐dose damage, where instead we observed enhanced ATM autophosphorylation, but muted Nbs1 phosphorylation. Likewise, both low‐ and high‐dose damage resulted in a biphasic profile of ATM activation, but phosphorylation of p53 on its ATM site, Ser‐15, tracked only with the late phase of ATM activity. These types of observations indicate complex coupling of kinases and substrates with phosphatases and events controlling protein production and turnover. Interestingly, a transcriptional modulator, homeo‐domain interacting protein kinase 2 (HIPK2), has emerged as a player in the core DDR response to DSBs, providing an additional level of regulation between ATM and p53 activities (Rinaldo et al, 2007; Winter et al, 2008; Puca et al, 2010). HIPK2 will be an important molecule to include in future studies investigating the complex regulatory interplay between p53 and ATM and life–death decision following DNA damage.

We made use of two distinct mathematical modeling approaches, PLSR and TI‐SWR, to facilitate the identification of significant relationships between measured signaling events and phenotypic responses. PLSR is an established relational modeling approach that we have used previously to identify signaling pathways and autocrine cascades involved in cytokine‐induced apoptosis in HT‐29 colon carcinoma cells (Janes et al, 2004, 2005). TI‐SWR, a new, complementary data‐driven relational modeling approach explicitly elucidates temporal relationships between particular signaling events and phenotypic responses (see Supplementary information). In these analyses, PLSR showed that the majority of the variance in the signaling data set was captured by two PCs that respectively associated with cell survival versus death or with cell‐cycle arrest versus progression, based on plots of the cellular‐response loadings in this same PC space. This PLSR analysis surprisingly implicated ERK1/2 activity in both the G1/S arrest and apoptotic cell death phenotype following DNA damage. This indication was then buttressed by the TI‐SWR analysis that discerned a strong role for ERK activity at 2, 4, and 8 h following DNA damage in maintenance of the early, transient cell‐cycle arrest in G1/S, and in promotion of apoptotic cell death from G1/S arrest following DNA damage. These predictions from the models were verified experimentally, confirming a role for ERK in cell‐cycle arrest and programmed cell death (Figure 7E). While we cannot rule out the possibility that ERK accomplishes this function indirectly by affecting the activity of other DDR molecules, a computational analysis does not suggest this to be the case (Supplementary Table S2).

These results are unanticipated in light of ERK1/2's well‐established role in promoting survival and progression into the cell cycle in the absence of DNA damage (Ballif and Blenis, 2001; Chambard et al, 2007; Meloche and Pouyssegur, 2007; Junttila et al, 2008). Furthermore, several studies have implicated ERK activity in these same pro‐survival and cell‐cycle progression responses in cells under DNA damaging conditions. Tsakiridis et al (2008), for example, reported that in non‐small cell lung tumors from patients undergoing high‐dose radiotherapy with or without chemotherapy, ERK activation was directly correlated with a poor response to treatment, while Kumar et al (2007) found that inhibition of ERK resulted in greater apoptotic cell death after γ‐irradiation in human endothelial cells that over‐expressed Bcl‐2. Work from Nishioka et al (2009) showed that inhibition of ERK caused greater growth arrest and apoptosis in NB4, HL60 and freshly isolated AML cells following administration of the DNA‐damaging chemotherapeutic agent cytarabine, while Hayakawa et al (1999) found that Erk inhibition enhanced cell death following cis‐platinum treatment of human ovarian carcinoma cell lines.

Other studies, however, have provided findings consistent with at least part of our results. A small but growing number of studies have indicated roles for Erk in promotion of cell death under a variety of conditions including following DNA damage (Zhuang and Schnellmann, 2006; Cagnol and Chambard, 2010). Inhibition of ERK, for example, has been reported to improve survival in A172 human glioma cells exposed to cisplatin (Choi et al, 2004), and in human hepatocarcinoma cell lines exposed to doxorubicin (Alexia et al, 2004). Similarly, Liu et al (2008) reported that Erk inhibition reduced the extent of apoptosis in rat cardiomyocytes exposed to doxorubicin, while Woessmann et al (2002) found that Erk inhibition reduced apoptosis in osteosarcoma and neuroblastoma cell lines treated with cis‐platinum. Furthermore, a few reports have indicated that Erk may be involved in G2/M control following DNA damage. Tang et al (2002) have suggested that ERK mediates G2/M arrest following low‐dose etoposide treatment and promotes apoptosis following high‐dose etoposide. Similarly, Yan et al (2007) concluded that ERK is necessary for the induction of a G2/M checkpoint following irradiation in MCF‐7 cells, where it is apparently required for activation of ATR, Chk1, and Wee1 in response to ionizing radiation.

Our study extends these findings by demonstrating that ERK1/2 is required for the maintenance of a G1/S arrest—a point of regulation that has not been previously reported—which is induced following exposure to low doses of doxorubicin, and that apoptotic cell death from this induced arrest is ERK dependent. These unexpected observations regarding the role of ERK in cell‐cycle arrest and death following genotoxic injury, as revealed through PLSR and TI‐SWR analysis of a systems‐level, time‐dependent signal–response data set, could have important implications for the application of MEK inhibitors in cancer therapy. It is important to note that these findings were made and validated in the U2OS cell line. While U2OS contain two wild‐type copies of the p53 gene, this cell line is known to be hypermethlyated at the INK4a/ARF locus, leading to lack of expression of both p16ARF and p14ARF (Park et al, 2002). G1/S cell‐cycle arrest in response to DNA damage may be expected to be less robust in the absence of p16/p14 molecules to reinforce this response. It is possible that a role for ERK in the maintenance of the G1/S arrest response, and death from this arrest, represents a ‘fail‐safe’ mechanism, and that a role for ERK in this response may be expected to be less pronounced in cell types that express p16/p14. Furthermore, the absence of a robust G1/S checkpoint in these U2OS cells may partially account for the unexpected doxorubicin‐induced death in the G1 state. Importantly, hypermethylation and mutations at the INK4a/ARF locus are found in many human cancers, including a large fraction of melanomas and carcinomas (Tannapfel et al, 2001; Sharpless and Chin, 2003; Cheung et al, 2009). Thus, Erk may fulfill this role in a variety of other tumor contexts. Our results imply that ERK pathway inhibitors should be used cautiously in combination chemotherapeutic treatments.

In future work, it will be crucial to determine under what genetic background and DNA‐damaging conditions ERK1/2 is pro‐survival and under what conditions it can promote apoptosis in order to effectively implement chemotherapeutic regimens that combine DNA damaging chemotherapeutic agents with MEK inhibitors. Computational approaches such as those used and developed in this work will aid in the unraveling of these complex dependencies. Our complementary use of PLSR and TI‐SWR for the analysis of a quantitative signal–response data set and identification of novel relationships between molecular signaling and cellular response represents a useful extension of established computational analysis methods. TI‐SWR is not limited to exploring signal–response relationships, but could potentially also be used to explore direct relationships between different molecular signals at the same or different points in time (c.f. Supplementary Table S2). Furthermore, this method seems to offer generalizability to any quantitative, time‐dependent signal–response data set and will be useful in the continuing quest to elucidate connections between signaling pathways and cell fate such as life–death decisions following exposure to DNA damaging chemotherapeutic agents.

Materials and methods

Reagents

Antibodies against the phosphorylated forms of histone H3 (Ser‐10) ATM (pSer‐1981) were obtained from Upstate Biotechnology. Antibodies detecting phosphorylated forms of H2AX (Ser‐139), p38 MAPK (Thr‐180/Tyr‐182), JNK (Thr‐183/Tyr‐185), ERK1/2 (Thr‐202/Tyr‐204), and Akt (Ser‐473) were obtained from Cell Signaling Technologies. Antibodies against total ATM and Nbs1 were obtained from Abcam, total H2AX and phospho‐p53 (Ser‐15) from R&D Systems, and phospho‐Nbs1 (Ser‐343) from Novus Biologicals. Antibodies against total p38 MAPK, p53, Akt, Cyclin A, and Cyclin B were obtained from Santa Cruz Biotechnology, against total JNK, total ERK1, and cleaved caspase‐3 from BD Biosciences, and against cleaved caspase‐3 and cleaved PARP from Molecular Probes. Antibodies against β‐actin were obtained from Sigma, IR‐conjugated anti‐mouse and anti‐rabbit secondary antibodies from Rockland Immunochemicals, and Alexa Fluor‐conjugated anti‐rabbit secondary antibodies were from molecular probes. TransAM transcription factor ELISA kits for detection of NF‐κB activity were obtained from ActivMotif. Doxorubicin, PD98059, aphidicolin, propidium iodide (PI), formaldehyde, methanol, and thymidine were obtained from Sigma, TNFα from Peprotech, protease inhibitors from Roche, and nitrocellulose from Millipore.

Cell culture

U2OS cells were obtained from ATCC were maintained at 21% oxygen and 5% CO2 in Dulbecco's modified Eagle medium supplemented with 10% fetal bovine serum, penicillin, streptomycin, 2 mM l‐glutamine, and used within 15–20 passages. For cell synchronization experiments, cells were seeded at 5 × 105 on 10 cm2 tissue‐culture dishes, and 24 h later supplemented with thymidine (2.5 μM final concentration) for 18 h. The first thymidine block was released by washing the plates three times with PBS, and incubating them in fresh thymidine‐free media for 12 h. A second thymidine block was then performed by re‐addition of thymidine to 2.5 μM followed by incubation for an additional 18 h. Media was aspirated, plates were washed 3 × with PBS, and replaced with fresh media in the presence or absence of 10 μM aphidicolin.

Sample collection and immunoblotting

To retain both adherent (live and dying) and floating (dying) cells, plates were removed from the incubator, placed on ice at the indicated time points, and media removed and placed in a pre‐chilled 15 ml Falcon tube. Plates were then rinsed with 4 ml of ice‐cold PBS, which was pooled with the reserved media. Cells were scraped in residual PBS with non‐lysing scrapers, and transferred to the same tube containing the media and PBS rinse. Following centrifugation for 5 min at 1000 g at 4°C, the supernatant was aspirated and the pellet resuspended in 1 ml ice‐cold PBS and transferred to 1.5 ml pre‐chilled Eppendorf tubes and microfuged for 5 min at 4000 r.p.m. at 4°C. The supernatant was discarded and the cell pellets snap‐frozen with liquid nitrogen and stored at −80°C.

Cell pellets were thawed on ice and lysed in 350 μl of fresh, ice‐cold lysis buffer containing 20 mM Tris–HCl, pH 7.4, 50 mM NaCl, 1% Triton X‐100, 50 mM NaF, 15 mM Na4P2O7, 2 mM Na3VO4, and protease inhibitors for 20 min with 15 rounds of sonication on ice. Aliquots totaling 75 μl were removed for protein quantification (micro‐BCA assay, Pierce) and ELISA assays. Sample buffer (55 μl of 6 × stock) was added to the remaining lysate, and 10–15 μg of each lysate, along with 10 μg of positive and negative control lysates for normalization, then analyzed by SDS–PAGE on 6–12% gels SDS–PAGE and transferred to nitrocellulose. Membranes were blocked with Odyssey Blocking Buffer (OBB) and stained with primary antibodies (in OBB+0.05% Tween) overnight at 4°C. All primary antibodies were used at a 1:250 to 1:1000 dilution. Membranes were then washed in PBS+0.1% Tween and probed with an IR‐conjugated anti‐mouse or anti‐rabbit secondary antibodies at a 1:10 000 dilution and visualized using an Odyssey Imaging System with IR laser excitation (Licor). Signal intensity was quantified using the Odyssey v2.1 software program and normalized to an anti‐β‐actin signal as a loading control.

Quantitative immunoblotting assay validation

All antibodies were analyzed for linearity of response in quantitative immunoblotting assays using positive control lysates in which the signaling pathway of interest was strongly activated by specific treatments. ATM, H2AX, and Nbs1 were maximally activated by incubation of U2OS cells with 45 nM neocarzinostatin for 1 h, p38 MAPK and JNK by treatment of cells with 50 J/m2 of UV irradiation, p53 by incubation with 10 μM doxorubicin, ERK and Akt by serum starvation followed by incubation with 100% FBS for 5 min, and NF‐κB by incubation with 100 ng/ml TNFα for 30 min. Lysates from untreated cells served as negative controls. Cell pellets were collected and lysed as described above, and samples containing 2.5–30 μg of each positive control lysate or 10 μg of negative control lysate analyzed by SDS–PAGE/immunoblotting as described. Signal intensity was quantified using Odyssey v2.1 software and the degree of linearity in the relationship between amount of activated/phosphorylated species loaded on the gel and intensity of resulting signal from the assay was evaluated (R2≥0.92 was used as a cutoff for antibody validation). Signal to noise of the assay was evaluated by comparing intensity of signal from negative and positive control lysate (S/N⩾2 was used as a cutoff).

Flow cytometry

Cells (both floating and attached) were collected at the indicated time points for analysis of cell‐cycle progression and apoptosis by pooling the media together with a single 2 ml PBS rinse as described above, followed by release of the attached cells by trypsinization. Samples were then centrifuged for 5 min at 1000 g and fixed according to an ethanol‐ (cell‐cycle outcomes) or formaldehyde‐ (cell‐death outcomes) fixation protocol. For ethanol fixation, the cells were prepared as a single‐cell suspension in 200 μl of ice‐cold PBS, and 2 ml of ice‐cold 70% ethanol was added dropwise with gentle vortexing. Samples were placed at 4°C for at least 2 h, centrifuged, resuspended in 200 μl ice‐cold PBS, and transferred to wells of a V‐bottom 96‐well plate. Plates were centrifuged for 4 min at 1800 r.p.m. at 4°C, washed once with PBS, and cells permeabilized in 200 μl of ice‐cold PBS containing 0.25% Triton X on ice for 15 min. The plates were re‐centrifuged as above, rinsed once with PBS+1% BSA, and stained with anti‐pHH3 antibody at a 1:100 dilution in 50 μl PBS+1% BSA overnight at 4°C. Following two washes with PBS+1% BSA, cells were incubated with an Alexa Fluor 488‐conjugated anti‐rabbit secondary antibody at 1:100 dilution for 1 h at room temperature. Cells were washed once with PBS+1% BSA, once with PBS, and counterstained for DNA content with 0.1 mg/ml PI in PBS containing 1 mg/ml boiled RNase A. Samples were stored at 4°C, protected from light, for at least 1 h and not more than overnight prior to analysis on a FACScalibur.

Formaldehyde fixation was performed in a similar manner by resuspending the cell/media/PBS pellet in 400 μl of 4% formaldehyde to a single‐cell suspension with incubation for 10–15 min. Cells were centrifuged for 5 min at 1000 g, washed with 800 μl ice‐cold PBS, and resuspended in 175 μl ice‐cold 100% methanol for at least 2 h for permeabilization. Cells were transferred to wells of a V‐bottom 96‐well plate, centrifuged for 4 min at 1800 r.p.m., rinsed in 200 μl PBS+0.1% Tween (PBS‐T), and resuspended in 50 μl PBS‐T containing 1% BSA (PBS‐TB) and anti‐cleaved caspase‐3 antibody at a 1:500 dilution. Following a 1‐h incubation at room temperature, 200 μl PBS‐T was added, the plates were centrifuged for 4 min at 1800 r.p.m., washed once with PBS‐T, and then resuspended in 50 μl PBS‐TB containing both an Alexa Fluor 488‐conjugated anti‐rabbit secondary antibody, and an Alexa Fluor 647‐conjugated anti‐cleaved PARP primary antibody, both at a 1:250 dilution, for 1 h at room temperature. PBS‐T (200 μl) was added to each well, the plates centrifuged for 4 min at 1800 r.p.m., washed once with PBS‐T, and resuspended in 200 μl PBS‐TB. Samples were stored at 4°C, protected from light, for at least 1 h and not more than overnight before analysis on a FACScalibur.

Relational modeling

For PLSR analysis, signaling and response measurements were formatted and pre‐processed using mean centering and variance scaling as described previously (Janes et al, 2005). PLSR analysis was performed using the SIMCA‐P 10.0 (Umetrics) software suite.

TI‐SWR.

TI‐SWR analysis was performed in Matlab using the stepwise function in conjunction with several scripts to automate analysis iterations and post‐processing of results. A dependent variable data vector (6 × 1) was generated to capture the change in a particular cellular response from time, t–1, to time, t, for a subset of the cellular responses measured under each of the six different treatment conditions by calculating

Embedded Image

The subset of cellular responses analyzed included (1) ΔcC3+/cPARP+6−12 h, (2) ΔcC3+/cPARP–6−12 h, (3) ΔcC3−/cPARP−6−12 h, (4) ΔG112−24 h, (5) ΔS12−24 h, (6) ΔG2/M12−24 h based on the large dynamic changes in these response measurements over the indicated time periods.

An independent variable matrix (6 × 14) was constructed for each of the 10 time points at which signaling data were measured, containing measurements of 14 signals at that time point under each of the six treatment conditions investigated. A forward stepwise regression algorithm was used to regress the dependent variable data vector (change in cellular response between two time points) against each of the independent variable matrices (signals) containing signaling data collected at times preceding the observed response measurement. All possible univariate regression models were sampled and the single most strongly correlated signal was used to initiate the model. For inclusion, a variable was required to display an R2⩾0.65, and the slope of the subsequent regression line had to be statistically significantly different from 0 (as confirmed by a P‐value ⩽0.05 using a linear regression t‐test). Next, all possible bi‐variate regression models that added one additional signal to the initial univariate model were sampled, and the signal whose addition to the model most significantly improved the correlation was added. This process of model building was continued until no further addition of signals significantly improved the signal–response correlation. This process was repeated for each of the six cellular responses listed above yielding a set of time point specific multiple linear regression models for each cellular response investigated.

Conflict of Interest

The authors declare that they have no conflict of interest.

Supplementary Information

Supplementary Information

Supplementary Figures 1‐6, Supplementary Tables 1 and 2 [msb20121-sup-0001.pdf]

Acknowledgements

We gratefully thank the Koch Institute core facilities for their technical support, especially Glen Paradis for help with flow cytometry. We thank Duaa Mohammad, Brian Joughin, Suzanne Gaudet, Kevin Janes, Pam Kreeger, and Ericka Noonan for advice and suggestions throughout the course of this project. This work was supported by NIH Grants CA112967 to MBY and DAL, and ES015339 and GM60594 to MBY. MJL received support from DOD fellowship BC097884.

Author contributions: This project was conceived and experiments designed by ART, MJL, GJO, LDS, and MBY. All experiments were performed and analyzed by ART. A subset of secondary assays and analysis was performed by ART and MJL. Computational analysis and interpretation was performed by ART, DAL, and MBY. The manuscript was prepared by ART, MJL, GJO, DAL, and MBY.

References

This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.

View Abstract