Impact of p53 Knockout and Topotecan Treatment on Gene Expression Profiles in Human Colon Carcinoma Cells: A Pharmacogenomic Study

The p53 protein is a central transcription factor activated in response to a variety of cellular stresses, including DNA damage, mitotic spindle damage, heat shock, metabolic changes, hypoxia, viral infection, and oncogene activation (1, 2). p53 can induce growth arrest and apoptosis, events that prevent the survival of damaged cells. p53 can also promote early senescence in response to unregulated mitogenic signaling (3, 4). The transactivation function of p53 is mediated through sequence-specific binding of its central domain to cis-acting elements within the promoters or introns of responsive genes. At present, more than 20 genes are known to be activated by p53, most of them in growth arrest or apoptotic pathways (5, 6). Other promoters (many of which are viral or growth stimulatory) are repressed by p53 (7). Consequently, the downstream effects of activating p53 are complex, and no single pathway mediates the full range of functions of p53.

To analyze the p53 dependence of molecular events after DNA damage, we compared gene expression changes in a p53 wild-type human colon carcinoma cell line, HCT-116 (p53+/+), with those in an isogenic p53 knockout (p53−/−; Ref. 8) after treatment with the topoisomerase I inhibitor topotecan. Topotecan (Hycamtin), a semisynthetic water-soluble derivative of camptothecin, is a clinically useful agent. It is approved for first-line therapy of cisplatin-refractory ovarian cancer and second-line therapy of small cell lung cancer (9). Like other camptothecins, topotecan converts topoisomerase I into a cellular poison by trapping topoisomerase I in a covalent complex with DNA. The cytotoxic lesions result from breaks, generated by collision of the complexes with DNA or RNA polymerase (10). We chose to use the p53 knockout in this study, rather than an overexpressing p53 transfectant, so that the p53 expression would be physiological.

To simultaneously study the effects of p53 status and topotecan treatment at different concentrations and time points, we developed (and present here) a new experimental design for microarray studies. We term it a cross-referenced network (Fig. 1). The most frequently used design for two-color microarray experiments simply compares each sample with a single internal reference sample by cohybridization. The network design uses internal reference samples, but it also provides the additional, global set of comparisons indicated in Fig. 1. Given this design, we were able to analyze the entire network of data by multiple regression in an appropriately weighted fashion to increase the statistical reliability of the results and conclusions in terms of statistical consistency test. The data were then displayed using a novel visualization, GEM,³ to identify transcripts that differed in expression level in relation to p53 status and/or drug treatment.

Isogenic p53-null (p53−/−) and wild-type (p53+/+) human HCT-116 colon carcinoma cells, kindly provided by Dr. Bert Vogelstein (8), were grown as monolayers in DMEM (Life Technologies, Inc., Gaithersburg, MD) supplemented with 10% fetal bovine serum (Hyclone, Logan, UT). Exponentially growing cultures at 80% confluence were treated for 1 h with 0.1 μm (LD) or 1 μm (HD) topotecan obtained from the Developmental Therapeutics Program, National Cancer Institute, NIH (Bethesda, MD). Cells were then washed twice with PBS (pH 7.4) resuspended in drug-free medium and harvested 0, 1.5, 3, or 6 h later. Total RNA was isolated from the cells using the RNeasy Midi Kit (Qiagen, Valencia, CA) and diluted to a concentration of 5 μg/μl in diethyl pyrocarbonate-treated water (Invitrogen, Carlsbad, CA). The RNA samples were then aliquoted and stored at −80°C until used for cDNA hybridization.

cDNA microarray hybridization was carried out based on a protocol developed in the Laboratory of Cancer Genetics, NHGRI,⁴ with minor modifications. In cDNA labeling reactions, 90 and 78 μg of total RNA were used for Cy5 and Cy3 labeling, respectively. Seventy-five μmol of Cy3-dUTP or Cy5-dUTP were used in each cDNA labeling reaction. Cy3- and Cy5-labeled cDNA samples (one experimental, the other a reference sample) were mixed and hybridized to pin-spotted cDNA microarray chips (Human Oncochips) developed in the Microarray Facility of the National Cancer Institute’s Advanced Technology Center. The microarrays contained 6720 cDNA clones representing approximately 6500 individual cancer-related genes. The complete gene list is posted at online.⁵ The hybridization data were acquired using a GenePix 4000 fluorescence scanner (Axon Instruments, Inc., Union City, CA).

Images were analyzed using a software package (F-SCAN) written by some of the authors [P. J. M., L. Y., and V. V. P. (11); Center for Information Technology, NIH]. The F-SCAN software is freely available online.⁶ The output of F-SCAN includes signal and background intensities for each spot in each channel, as well as several quality control measures. The JMP 4.0 statistical package (SAS Inc., Cary, NC) was used to inspect the data for artifacts, apply gene selection criteria, and calculate the hierarchical clustering.

Measured raw signal intensities were log-transformed: I₁, I₂ = log₁₀ (raw intensity) for channels 1 and 2, respectively. Log-intensity ratios related to spot colors were calculated as L = I₁ − I₂. Systematic drift in color was occasionally seen across the 16 blocks of the array, corresponding to the 16 different pins used in the printing process. To ameliorate this effect, the median log-intensity within each block was subtracted and a pin-effect corrected log-intensity, L′, was computed. For a particular spot, S, printed by pin p = pin (S), L^′_S = L_S – median{L_T:pin (T) = pin (S)}. Because L′ sometimes varied as a function of spot intensity, it was necessary to reduce this effect by the following procedure. In the plot of L′ versus average log-intensity, I = (I₁ + I₂)/2, a smooth, running median curve modeled the color trend as a function of intensity. For a particular spot S, the corrected log-ratio value was computed as L^″_S = L^′_S – running_median (I_S), where running_median (I) is the value on the smooth curve corresponding to average log-intensity I.

The 12 distinct cDNA preparations were as shown in Table 1.

Each chip was hybridized using two cDNA preparations, each labeled with either Cy3 or Cy5. Every such comparison of preparation i with preparation j was also performed using reciprocal labeling, and the pair was replicated n (i, j) times. To compensate for any possible systematic bias caused by the chemical labeling difference between the two dyes, the experimental log-ratio, L″(i, j, k), and that from the reciprocal experiment, − L″(j, i, k), for each replicate k were averaged to give the final estimate of the log-expression ratio y(i, j) for that comparison.

When an observation fell below quality control cutoffs (based on the number of pixels in a spot and the intensity level above background), it was dropped from the average. When the reciprocally labeled data were compared with replicate data for the p6LD versus p0 comparison, we noted that agreement between reciprocally labeled data were better, in general, than would be based expected from the agreement of similarly labeled replicates. This fact supports the effectiveness of the block normalization/color correction procedure used in this study.

Initially, each sample was hybridized against its corresponding zero-time sample, producing 10 pairwise comparisons. Subsequently, the p53 wild-type cells (p) and p53 null cells (m) were compared with each other directly at t = 0, 3, and 6 h at HD, bringing the total to 13 pair-wise comparisons. The resulting network of comparisons (Fig. 1) can be used to deduce the relative expression difference between any two treatment conditions in the study, not just those cohybridized on the same array slide.

A gene was judged to differ significantly for a particular comparison of preparation i with preparation j if the minimum log-ratio over all replicates, k, min{L″(i, j, k), − L″(j, i, k)} exceeded a fixed cutoff, c, or the maximum over k, max{L″(i, j, k), − L″(j, i, k)} < −c. For this study, we established a log-ratio cutoff of 0.2, corresponding to about a 1.6-fold (i.e., 60%) difference. Using a newly devised consistency test, it became possible to calculate the probability of exceeding that cutoff. In our study, about 2.5% of genes exceeded the selected cutoff (in absolute value) in a single comparison. If there were indeed no true changes of expression for any genes, the variation of the observed expression would arise only from random experimental noise (the null hypothesis). Assuming that the noise acted independently in duplicate experiment, one would expect (under the null hypothesis) less than (2.5%)² = 0.06% of the genes to exceed the cutoff in duplicate comparisons, (2.5%)³ = 0.000015 in triplicates, and so forth. For duplicate and quadruplicate experiments with 6700 genes, one therefore expects fewer than 4.0 and 0.003 false positives, respectively. In most comparisons reported here, the false discovery rate (the number of false positives expected/number detected) was less than 20% by this approach. Altogether, 809 of 6720 genes were seen to change by 60% or more in all replicates of at least 1 of the 13 comparisons. One should interpret these false discovery rates with caution because the replicate hybridization sometimes used cell lysates from the same flasks, although with different labeling, and thus might not have been fully statistically independent.

Using the full network of treatment comparisons (Fig. 1), we determined log₁₀ relative expression levels relative to a single preparation, m0, the p53-null cells at t = 0. The experimental design contained multiple, redundant comparisons that are reflected as cycles in the graph. For example, p0 can be compared directly with m0 or, alternatively, can be compared with m0 indirectly through p3HD and m3HD. This redundancy provides an opportunity to refine or “average” the results further. To reduce the 13 averaged comparison measurements to estimates of relative log-ratios, we used least-squares estimation with weighting by the number of available replicates.

Let β(i) be the concentration of a particular mRNA species in the i^th cDNA preparation, i.e., β(1) = m0 log-concentration, β(2) = m1.5LD log-concentration, and so forth. We can write a system of equations y(i, j) = β(j) − β(i) for each comparison of preparation i with preparation j. The least-squares estimation procedure can best be presented in the following vector notation. We define Y to be the vector of the 13 log ratios for each comparison, indexed by preparation numbers,

and β, the vector of 12 unknown log-concentrations in each cDNA preparation,

The system of equations relating Y to β can be written compactly as follows after defining the design matrix X with 13 rows and 12 columns (showing only the non-zero entries):

The system of equations is simply Y = Xβ. Because this system is not of full rank, we must add a constraint on β to obtain a unique solution. This can be accomplished by setting β(1) = m0 = 0, which is equivalent to treating all other β(i) values relative to that for m0. Operationally, we augment X with a single row having a 1 in the first column only and augment the vector Y with a single element, zero. Finally, to account for the different numbers of replicates, n(i, j), contributing to each Y value, we define the diagonal weight matrix

The factor of 2 reflects the presence of two reciprocally labeled chips in each replicate. The least-squares solution to the system of equations is given by:

a 12-component least-squares estimate of the desired set of relative log-concentrations. The least-squares computation is carried out for each gene represented on the chip. The RMS error from the least squares regression is a measure of the “lack of fit” of the model, i.e., how much inconsistency one finds in the network of treatment comparisons. The RMS error is computed as shown below.

No single criterion has been found that can effectively remove all outlying or bad data in the experiment. Therefore, each gene in the set of 809 that showed significant differences in at least one pairwise comparison was further required to meet three conditions: (a) a sufficient number of the comparisons schematized in Fig. 1 were present for that gene to determine a solution β̂ to the least-squares equations (1); (b) the RMS error for the gene fell below a specified level, in other words, the gene must give self-consistent results across all available comparisons; and (c) the responsiveness (SD of the expression pattern) of the gene exceeded a specified level, i.e., the gene must clearly be among those that responded to the treatment stimuli. The RMS cutoff was set so that genes with significantly more than 30% error would be rejected. The upper 95% quantile for a χ² distribution with 2 degrees of freedom is 5.99. The corresponding cutoff for the RMS error is therefore 5.99*log₁₀(1.3) = 0.68, which was rounded to 0.7. The responsiveness was measured as the SD of the relative concentration estimates, as shown below.

The responsiveness cutoff was empirically set at 0.12, corresponding to about a 30% variation across the cDNA preparations. Although somewhat arbitrary in nature, these filters effectively removed many obvious outliers and much of the lower quality data, leading to a set of 167 consistently and significantly changing, highly responsive, high quality spots, which formed the basis of all further analyses.

The experimental design (Fig. 1) included three factors: p53 presence or absence; drug treatment at low concentration (LD) or high concentration (HD); and duration of treatment (0, 1.5, 3, or 6 h). The full pattern of response is given by the vector β calculated previously. Complex patterns of response could be expected (e.g., overexpression at 1.5 h followed by underexpression at 3 h, and so forth). Various linear combinations of the estimated β parameters are likely to be more interpretable than the original vector. In general, the vector space of possible patterns now spans an 11-dimensional subspace of a 12-dimensional space. We defined an orthogonal rotation of this space to that after rotation; the second dimension represents the p53 effect, the next five dimensions include all treatment (dose and time) effects, and the remaining five dimensions define interactions between p53 and treatment effects. The selected orthogonal rotation matrix E is given by:

where D = sqrt(inv(diag([12, 12, 60, 60, 12, 4, 4, 60, 60, 12, 4, 4]))), a diagonal normalization matrix.

The new vector of effects is computed as an orthogonal rotation of the original vector β.

Names may be given to these effects as follows: Mean; p53; treatment; HD; early time at LD; late time at LD; early time at HD; p53-treatment interaction; p53-HD interaction; p53-early time at LD interaction; p53-late time at LD interaction; and p53-early time at HD interaction; respectively. The p53 effect is given by Effect(2), or equivalent, the difference in expression between p53+/+ and p53−/− cells given the same treatment:

The root-12 appears here to normalize the result and make it comparable in scale with the other effects.

Because the treatment effects are spread over the next five dimensions [Effect(3) through Effect(7)], we chose to define a composite treatment effect whose magnitude was the RMS of the five treatment effects, and its sign was determined by Effect(3).

The complex interaction components were neglected because most of the variability of the genes in this study could be interpreted as a function of the two variable treatment effect and p53 effect.

Depicting the results of the study in only two variables [p53 effect and treatment (Rx) effect] allows one to view the overall landscape of gene expression for the entire set of genes being studied. We call this depiction (Fig. 5) a GEM. By analogy with cartographic depictions of the world, the roles of longitude and latitude are played by p53 effect and treatment effect, respectively. Genes found in a particular region of the map have similar expression profiles and may also be associated in a more fundamental, physiological way, perhaps in some cases occurring in a common pathway. This map helps us to identify the most prominent landmark genes, which are found at the periphery of the map. Genes with little or no response to the treatment and no p53 dependence are found close to the center of the map and are omitted to simplify the picture.

The genes can be classified according to whether their expression shows a positive, negligible, or negative p53 effect (i.e., in terms of p53+/+ minus p53−/− levels). Likewise, a positive, negligible, or negative treatment effect indicates the behavior of gene expression as a result of topotecan treatment (at various doses and times studied here). The territories, named for convenience according to the points of the compass, are: (a) N, p53 independent, positive treatment effect; (b) NE, positive p53 effect, positive treatment effect; (c) E, positive p53 effect, treatment independent; (d) SE, positive p53 effect, negative treatment effect; (e) S, p53 independent, negative treatment effect; (f) SW, negative p53 effect, negative treatment effect; (g) W, negative p53 effect, treatment independent; and (h) NW, negative p53 effect, positive treatment effect.

Gene expression levels found in N and S are essentially independent of p53 status but responsive to treatment, whereas gene expression levels found in E and W are essentially independent of treatment but differ between the p53+/+ and p53−/− cells.

Results for the 167 selected genes, rendered as a 167 × 12 matrix of log-ratio comparisons with the m0 condition, were hierarchically row-clustered using Ward’s method implemented in JMP (SAS Inc.). Clustering the rows of the matrix (corresponding to genes) had the effect of bringing similar expression patterns close together in the figure. The reordered matrix was then color-coded, with red representing overexpression and green representing underexpression in comparison with m0. The present application of the CIM visualization (12, 13) is novel in that all expression ratios were first averaged over reciprocally labeled, replicate comparisons and then reconciled, and each treatment condition was compared with a single control condition. The underlying data were obtained from 35 separate pairwise comparisons but were reconciled into just 12 expression ratios/gene, thus reducing the noise inherent in separate measurements. Finally, each gene name was annotated with the region of the GEM in which the gene fell. There was good correspondence between these two representations of the expression profiles for the entire study.

The expression levels of selected transcripts were verified by real-time RT-PCR using an ABI Prism7700 Sequence Detection System (PE Applied Biosystems, Foster City, CA). TaqMan primers and FAM-TAMRA-labeled probes were designed for each gene of interest using Primer Express 1.0 software (PE Applied Biosystems) with the manufacturer-specified parameters.

Cells were incubated with 0, 0.1, or 1 μm topotecan and then washed twice with PBS (pH 7.4) and harvested. Cell lysates were prepared as described previously (14). One hundred-μg protein samples were then separated by SDS-PAGE (12% polyacrylamide gel) and transferred into Immobilon membranes (Millipore, Bradford, MA). p53 and PCNA proteins were identified using anti-p53 and anti-PCNA primary antibodies (Dako, Carpinteria, CA). The reactive bands were then visualized using an enhanced chemiluminescence detection system (New England Nuclear Life Products, Boston, MA). S100A4 was detected using anti-S100A4 primary antibody (Dako).

The cDNAs obtained from all cell samples, HCT-116 (p53+/+) and HCT-116 (p53−/−), treated and untreated, were hybridized to microarrays in reciprocal pairs as indicated by the arrows in Fig. 1. We analyzed the data using two approaches. The first was a pairwise analysis in which genes were selected based on a consistent difference (>1.6-fold) between samples for at least one pair of conditions (as explained in “Materials and Methods”). This approach led to the selection of 809 genes. The second approach was a multivariate analysis in which the 809 genes selected were further analyzed based on consistency of results, as described in “Materials and Methods.” This approach led to the selection of 167 genes, which were then classified based on both p53 and treatment dependence.

The data summarized in Table 2 show the number of genes that were expressed differentially as a function of drug concentration and time of exposure. When p53+/+ cells were treated with LD topotecan, a total of 24 genes were induced, and 16 genes were repressed within the 6 h of drug treatment. In the p53−/− cells, 11 genes were induced, and only 4 genes were repressed. This difference was even more pronounced when the cells were treated with HD topotecan; 178 genes were induced and 483 were repressed in p53+/+ cells, but only 20 were induced and 49 were repressed in p53−/− cells. Thus, the number of differentially expressed genes was much greater in the p53+/+ cells than in the p53−/− cells, indicating that a large number of genes are controlled at the transcriptional level directly or indirectly by p53. It is also noteworthy that the number of genes induced by HD treatment in p53+/+ cells relative to p53−/− cells was much greater at 6 h than at 3 h. The same was true for the repressed genes, indicating that p53-dependent expression/repression (as reflected at the RNA level) takes some time (>3 h) to develop fully.

We next asked whether the observed changes in gene expression were, in fact, correlated with differences in p53 protein expression. The Western blots in Fig. 2 show that p53 protein levels increased by 1.5 h with LD topotecan without further detectable increase at 3 or 6 h. The level of p53 protein was much higher after treatment with HD, and the maximum level of p53 protein was detected at 6 h, indicating that the change in p53 protein expression after topotecan treatment is correlated with that of RNA expression.

The 809 selected genes were filtered further to yield 167 genes that showed substantial changes under some conditions of drug concentration, time of exposure, or p53 status and also showed overall consistency of responses. The results are displayed in the form of the CIMs (12, 13) shown in Figs. 3 and 4. These 167 differentially expressed genes are also listed in Table 3 in the order in which they appear in the CIM (Fig. 3). The most striking feature of the CIM is the difference in gene expression patterns between the two cell types, clearly demonstrating that the presence of functional p53 profoundly affected the stress response profiles of a large number of genes.

Fig. 4 focuses on the gene clusters marked α, β, and γ in Fig. 3. The expression patterns of these genes clearly distinguish the responses of p53+/+ and p53−/− cells. The transcripts in cluster α were more strongly expressed in p53+/+cells, and in these cells, some transcripts tended to decrease slightly with topotecan treatment, e.g., HNRPK (spot 25), SLC7A5 (spot 32), and ENO1 (spot 37).

Cluster β includes p53-dependent transcripts whose expression increased with treatment. Genes in this cluster include the known p53-responsive transcript CDKN1A/p21^Waf1/Cip1 (spot 58), as well as some members of the TGF-β superfamily [PLAB (three spots; spots 59, 60, and 62) and MADH6/SMAD6 (spot 61)] and the proapoptotic gene SIVA (spot 63), whose p53-dependence has not been reported previously.

Cluster γ consists of transcripts expressed at higher levels in the p53−/− cells than in the p53+/+ cells after LD topotecan exposure. After HD exposure, the expression levels of these genes were uniformly low. Expression in the p53−/− cells after LD exposure peaked at 1.5 or 3 h and then declined (Fig. 4). Cluster γ contains genes related to the HSP70 family [HSPA1A (spots 96 and 97), HSPA10 (spot 100), and HSPA1L (spots 101 and 102)].

To test the microarray results, we performed real-time RT-PCR for a cluster of three selected genes (indicated by a blue arrow in Fig. 3) whose expression was drug dependent and greater in the p53−/− cells than in the p53+/+ cells. The same RNA batches were used for both microarrays and real-time RT-PCR. As indicated in Table 4, the two methods for expression analysis were in good agreement.

For concurrent examination of the two experimental variables, namely, drug treatment and p53 status, we devised a novel approach for analysis and visualization. The results are displayed as a two-dimensional GEM (see “Materials and Methods”). The GEM (Fig. 5) shows how the 167 differentially expressed genes can be grouped on the basis of treatment effect and p53 effect. The GEM regions where the individual genes fall are listed in the second column in Table 3. The last two columns in Table 3 contain the coordinates of gene locations in the GEM of Fig. 5.

The Northeast (NE) region contains genes that were expressed at higher levels in p53+/+ cells than in p53−/− cells or that tended to be up-regulated by topotecan treatment. As expected, CDKN1A/p21^Waf1/Cip1 (spot 58) fell into this category, as did other known p53 target genes [PPM1D/WIP1 (spot 13), TNFRSF6/Fas (spot 14), and ATF3 (spots 53 and 54); see Fig. 5 and Table 3]. Other genes not previously linked to p53 transcriptional activity appear in the NE region. Among them are a set of genes related to the TGF-β growth regulatory pathway. As noted earlier, these genes [PLAB (2 spots; spots 59 and 62) and MADH6/SMAD6 (spot 61)] and the proapoptotic SIVA gene (spot 63) belong to the β cluster. Also in the NE region is KIAA0838 (spot 60), a chimeric clone that includes sequences for both glutaminase C and PLAB. The similar expression patterns of PLAB (two separate spots; spots 59 and 62), KIAA0838 (spot 60), and MADH6/SMAD6 (spot 61) transcripts are shown quantitatively as functions of time and drug concentration in Fig. 6. Another gene in the NE region is the proapoptotic SIVA gene (spot 63). The time course and dose response of SIVA showed behavior parallel to that of the TGF-β pathway genes (Fig. 6).

The opposite of the NE region is the Southwest (SW) region, which contains transcripts expressed at lower levels in p53+/+ cells and down-regulated by drug treatment. This region contains ENG/CD105 (endoglin; spot 148), a TGF-β type III receptor (15, 16). The expression behavior of this transcript is inverse to that of PLAB and SIVA, as shown in Fig. 6. Another transcript in this region is HS3ST1 (spot 135 in Table 3). It codes for a heparan sulfate glucosamine sulfotransferase, which may be involved in transmission of growth factors from the extracellular matrix to the cell surface.

The East (E) and West (W) regions contain p53-dependent transcripts that are not substantially responsive to drug treatment. In the E region, the S100A4/mst1 transcript (spot 56) has the most intensely p53-dependent expression (>25-fold). S1004A4/mst1 is a M_r 11,000 protein that belongs to the S100 family of Ca²⁺-binding proteins, different members of which have diverse cellular functions (17). Its involvement in tumor metastasis has generated recent interest (18).

Differential expression of S100A4 was confirmed by immunoblots (Fig. 7), which demonstrated a strong signal in the p53+/+ cells and no detectable signal in p53−/− cells. This observation suggests a close association between S100A4 expression and p53 function.

The W region, which corresponds to drug-independent transcripts expressed more in the p53−/− cells than in the p53+/+ cells, contains the M_r 70,000 HSPs [HSPA1A (spots 96 and 97), HSPA10 (spot 100), and HSPA1L (spots 101 and 102)] already mentioned above in connection with cluster γ in Figs. 3 and 4. Consistent with this observation, Kannan et al. (6) have reported down-regulation of HSP70 in cells that overexpress a temperature-sensitive p53. In that study, the down-regulation might have been associated with a nonphysiological level of expression of p53; however, our results were obtained with physiological expression of wild-type p53.

The North (N) and South (S) regions include transcripts that were enhanced or suppressed, respectively, by topotecan exposure, independent of p53 status. The N region includes the FOS/JUN genes (spots 47 and 49–51), consistent with the role of activator protein 1 in stress responses (19). The N region (Fig. 5) also includes CFLAR/cFLIP (spot 106) and casein kinase I (CSNK1G2; spot 108), both of which can modulate apoptosis. cFLIP (CFLAR), a potent inhibitor of the death receptor pathway (type I apoptosis), binds to Fas-associated death domain and blocks activation of FLICE (caspase 8; Ref. 20). Our finding that cFLIP expression can be induced by topotecan independent of p53 status (Fig. 5; see also Fig. 3, spot 106; and Table 3) is novel (21). Casein kinase I was recently shown to inhibit apoptosis by preventing the activation of Bid and therefore inhibiting mitochondrion-dependent apoptosis (type II apoptosis; Ref. 21).

Among the transcripts in the S and SW regions (i.e., down-regulated by topotecan treatment, independent or negatively dependent on p53 status) are ARPC1B (spot 134; a gene related to actin), DGS1 (spot 131; DiGeorge syndrome critical region gene), and GRAVIN (two spots; spots 122 and 123 that clustered next to each other). GRAVIN was first identified as an autoantigen in some myasthenia gravis patients. It functions as a scaffold protein for protein kinases A and C and interacts with actin (22). Perhaps future research will indicate why these S and SW region transcripts were down-regulated, but, at present, we cannot discern a pattern.

p53, a key tumor suppressor gene, is mutated in the majority of human cancers (23). Successful outcome of chemotherapy and radiotherapy, in many cases, depends on functional p53 (1, 2, 24, 25). Therefore, elucidation of the function, regulation, and molecular interactions of p53 is of great importance for cancer therapy. Because p53 is a transcription factor with an expanding repertoire of genes that are known to be directly or indirectly under its control, global analysis of gene expression profiles represents the best approach for studying the p53 response (6, 25, 26).

In previous studies, patterns of gene expression were observed following overexpression of wild-type p53. For example, the Vogelstein group used serial analysis of gene expression (SAGE) in a colorectal carcinoma cell line to analyze inducible overexpression of wild-type or mutant p53 (25). Of the 9,954 transcripts identified, 34 (0.34%) were markedly increased in p53-overexpressing cells. Many of those 34 genes had not previously been shown to be regulated by p53. The Levine group (26) used microarrays to study the response of p53-regulated genes in a human colon carcinoma EB-1 cell line stably transfected with a construct that included a zinc-inducible promoter. Using cluster analysis, they found that response patterns depended on whether the genes were induced by γ-radiation, UV-radiation, or the inducible p53 gene. Recently, Wang et al. (27) showed that 1,501 of 33,615 genes (4.5%) that contained the p53 consensus binding sequence responded positively to p53 expression. Although this study used different criteria to define genes with different experiential levels, this number is close to our finding here that 167 of 6,500 genes (2.5%) respond consistently and significantly to both p53 expression and drug treatment.

In the present study, we compared drug-induced transcriptional responses in wild-type HCT-116 human colon carcinoma cells and an isogenic p53-deleted line (8). Comparison of isogenic cell lines, in combination with pharmacologically induced stress, allowed us to study not simply the effects of drug treatment or the effect of p53 but the combination of both treatment and p53 (see Fig. 1) in a physiologically driven system.

The need to analyze both genotypic/phenotypic and pharmacological variables simultaneously prompted us to develop new types of experimental design and data analysis. To encompass a wide range of potential responses in both the p53 wild-type and knockout cells, it seemed important to study different doses and treatment times. Hence, the design included comparisons among 12 conditions of drug concentration, time, and p53 status. Included were standard, direct pairwise comparisons (i.e., on a single slide) of treated and untreated samples for each cell type. The design also included direct comparisons of the p53 wild-type and knockout states in baseline and treated states. Overall, this cross-referenced network design (Fig. 1) provided built-in redundancy. It also enabled us to obtain quite robust estimates of the relative expression levels by multiple linear least-squares regression calculations for the entire network of experimental conditions. The experimental design and approach to analysis represent advances over other studies that use a single “reference” sample in one channel of a two-label microarray hybridization experiment. Interpretation focused on 167 genes whose expression levels responded most clearly in a consistent and significant way (see Fig. 3 and Table 3).

A novel two-dimensional visualization map, the GEM, identified the components of observed gene expression differences attributable to each of the two major factors in the design: p53 status; and treatment with the topoisomerase I-targeted drug topotecan. For simplicity of description, we can identify eight possible generic outcomes (“territories” in the GEM), each representing a different response pattern, as indicated in Fig. 5. If a set of genes from a CIM cluster (Figs. 3 and 4) falls largely within a single territory in the GEM, we may be led to the hypothesis that some of those genes operate in a common pathway. Also genes up-regulated in one quadrant might operate coherently with genes down-regulated in the complementary quadrant.

Two examples of coherent gene responses are the TGF-β and Jun/Fos clusters. Four genes in the TGF-β pathway appear in the NE quadrant of the GEM (Fig. 5) and in the β cluster (Figs. 3, 4, and 6): PLAB (spots 59 and 62); KIAA0888 (spot 60); and MAD6 (spot 61). The up-regulation of these genes was p53 dependent. This induction of the TGF-β pathway provides an example of the indirect effects of p53 on gene transcription by activation of other transcription factors. Also up-regulated was the Jun/Fos pathway. Three Jun B spots (spots 49, 50, and 55) and the FosB gene (spot 51) are in the N quadrant of the GEM (Fig. 5) and cluster next to each other in the CIM (spots 49–51 and 55 in Fig. 3). The Jun/Fos up-regulation was dependent on topotecan treatment but independent of p53 status.

Using the GEM analysis, we find that p53 regulates (directly or indirectly) a large number of genes that can induce apoptosis or cell cycle arrest. Some of the p53-dependent responses, to our knowledge, have not been reported before. Among the genes that exhibited p53 dependence in the positive direction, we identified two recently described proapoptotic genes (SIVA and S100A4) and several genes from the TGF-β pathways. SIVA (spot 63) binds to CD27, a member of the TNF receptor family (26, 27) and induces apoptosis after phosphorylation by the Abl-related gene (ARG) tyrosine kinase in response to oxidative stress (28). Thus, SIVA could be a proapoptotic p53 effector. To the best of our knowledge, the p53 dependence of SIVA expression has not been reported previously.

The metastasis-associated protein (S100A4/mst1) belongs to the S100 family of small calcium-binding proteins (for review, see Ref. 29) recently reported to cooperate with p53 to induce apoptosis (30, 31). Calcium binding induces conformational changes in the S100 protein structure, allowing interaction with target proteins including p53 (30, 31), the protein tyrosine phosphatase Liprin bl (32), non-muscle myosin heavy chain, and another S100 family member, S100A1 (see the references in Ref. 32). Binding of S100A4 to the COOH-terminal domain of p53 inhibits p53 phosphorylation by protein kinase C (but not by CKII) and modulates p53 transcriptional activity in a gene-specific manner (32). S100A4 stimulates p53-mediated Bax transcriptional activation and suppresses transcription of the antiapoptotic p21^Waf1/Cip1 gene (31). These gene-specific transcriptional effects provide an explanation for the observed cooperation between wild-type p53 and S100A4 in the induction of an apoptotic response (31). Thus, S100A4 might contribute to the apoptotic response to topotecan in wild-type HCT-116 cells, once p53 levels increase. Our finding that both protein and mRNA levels of S100A4 were markedly higher in wild-type HCT-116 cells than in HCT-116 p53 knockout cells was unexpected, considering that Grigorian et al. (31) observed an inverse relationship between p53 wild-type status and S100A4 expression in 26 tumor-derived cell lines, not including HCT-116 cells.

Among the genes that were down-regulated by p53, we identified the antiapoptotic protein chaperone HSP70 genes. These genes [HSP1A (spots 96 and 97), HSP10 (spot 100), and HSP1L (spots 101 and 102)] fell in the γ cluster (Figs. 3 and 4) and in the W quadrant of the GEM (Fig. 5). HSP70 has recently been shown to suppress apoptosis by blocking the apoptosome (33). Two genes with possible antiapoptotic activity fell in the p53-independent, drug-responsive region (N region): cFLAR/cFLIP (spot 106; Ref. 20); and CSKNK1G2 (spot 108; Ref.21). In conclusion, these experiments have added, at least incrementally, to our understanding of how p53 activation can precipitate cell death in response to pharmacological stress. Perhaps as important, investigators elsewhere who approach these carefully generated data sets with different questions and different perspectives will be able to observe relationships that we cannot now observe. More generally, the cross-referenced network experimental design and methods of analysis introduced here can be applied widely to the simultaneous study of genotypic/phenotypic and pharmacological variables.

We thank the Advanced Technology Center of the National Cancer Institute for providing the Human Oncochips. HCT-116 and an isogenic p53 knockout cell line were kindly provided by Dr. Bert Vogelstein (Johns Hopkins Medical Institute, Baltimore, MD).