Capturing drug responses by quantitative promoter activity profiling.

Quantitative analysis of cellular responses to drugs is of major interest in pharmaceutical research. Microarray technologies have been widely used for monitoring genome-wide expression changes. However, this approach has several limitations in terms of coverage of targeted RNAs, sensitivity, and quantitativeness, which are crucial for accurate monitoring of cellular responses. In this article, we report an application of genome-wide and quantitative profiling of cellular responses to drugs. We monitored promoter activities in MCF-7 cells by Cap Analysis of Gene Expression using a single-molecule sequencer. We identified a distinct set of promoters affected even by subtle inhibition of the Ras-ERK and phosphatidylinositol-3-kinase-Akt signal-transduction pathways. Furthermore, we succeeded in explaining the majority of promoter responses to inhibition of the upstream epidermal growth factor receptor kinase quantitatively based on the promoter profiles upon inhibition of the two individual downstream signaling pathways. Our results demonstrate unexplored utility of highly quantitative promoter activity profiling in drug research.CPT: Pharmacometrics & Systems Pharmacology (2013) 2, e77; doi:10.1038/psp.2013.53; published online 25 September 2013.

In the development of new drugs, monitoring cellular responses to drug candidates is a fundamental approach for assessing the efficacy and safety of using these drugs. A range of molecular profiling approaches, such as mass spectrometry, antibody-based proteomics, quantitative reverse transcriptase–polymerase chain reaction, and DNA microarray, has been used for monitoring responses at the molecular (protein and RNA) level.[1,2,3,4] For instance, quantification of mRNA abundance is an effective way for monitoring gene expression changes in response to the drugs.[5,6,7,8] Microarray-based technologies have been widely used for monitoring such changes on a genome-wide scale[9] also in a context of drug effects on particular signaling pathways.[10,11,12] However, microarrays have several limitations in terms of coverage of targeted RNAs, sensitivity, and dynamic quantitative range because they rely on predesigned oligonucleotide probes and hybridization-based detection.[13] The quantitative limitation has forced researchers to use semi-quantitative interpretation of gene expression changes, for example, by rank-based analysis.[14,15,16,17,18] RNA-seq is one of the latest techniques for profiling the transcriptome[19,20] by sequencing random fragments of RNA; here, most of the protocols rely on second-generation sequencers and polymerase chain reaction–based amplification. Cap Analysis of Gene Expression (CAGE) is an alternative method for quantifying the transcriptome by sequencing the 5′-end of RNAs,[21] and transcription start site profiles based on a polymerase chain reaction–dependent CAGE protocol are used as a reference of promoter activities in quantitative modeling based on multiple epigenetic marks in the ENCODE consortium.[22] Recently, we improved on the CAGE method by adapting it to a third-generation (single-molecule) sequencer, which allowed us to avoid any amplification steps from the library preparation to the sequencing reaction, suggesting that the resulting read counts represent the absolute number of observations of RNA presence.[23,24] In this article, we ask whether cellular responses can be modeled quantitatively from the aspect of the transcriptome, in particular, promoter activities. We demonstrate quantitative modeling based on accurate quantification of subtle cellular responses induced by low-dosage drug treatment.

Results

Promoter activity profiling of cellular responses to drugs

We monitored the effects of U0126, wortmannin, and gefitinib on human breast cancer MCF-7 cells using the quantitative and genome-wide promoter profiling method. U0126 and wortmannin are specific inhibitors of the Ras-ERK and phosphatidylinositol-3-kinase (PI3K)-Akt pathways, respectively (). Gefitinib is a potent inhibitor of the epidermal growth factor receptor (EGFR) kinase and mainly inhibits the Ras-ERK and PI3K-Akt pathways located downstream of this receptor. After determination of dosage of these drugs that show significant but not saturating effects on the cells (see Supplementary Figure S1 online), we prepared three replicate samples followed by CAGE profiling. On average, we obtained ~14 million reads mapped on the reference genome per sample. By aggregation of neighboring transcription start sites (see Supplementary Methods online for detailed parameters and thresholds), we defined 10,298 promoters with characteristics consistent with those in a previous research[23] (see Supplementary Figures S2–S4 and Table S1 online). Of note, even when we treated with drugs at low concentrations, promoter activities across triplicate samples were highly reproducible (average of three drug samples and standard deviation of Pearson's correlation coefficient = 0.9984 ± 0.0016; a scatter plot of the biological replicates is shown in and Supplementary Figure S5 online). By differential analysis comparing with no drug treatment, we identified 139, 168, and 157 promoters significantly affected by U0126, wortmannin, and gefitinib treatment, respectively (false-discovery rate <2%; see Supplementary Table S2 online). Although drug treatment of the cells induces changes in gene expression, the magnitude of such changes is, in general, not very extraordinary because fundamentally they are of the same cell type.[5,6,9] U0126 significantly suppressed EGR1 and FOS, target genes of the Ras-ERK pathway,[25] wortmannin upregulated ErbB3, which was previously reported to be upregulated in a PI3K-Akt and FoxO3a-dependent manner,[26] and gefitinib suppressed CCND1 promoter activity, one of the target genes of the EGFR pathway.[27] Notably, FOS was not detected as a differentially expressed gene following U0126 treatment as by microarray profiling of the same RNA samples (). Taken together, the promoter activity profiles that we generated in this study are highly reproducible across biological replicates, consistent with previous studies, and could quantify modest changes in the transcriptome at submaximal drug dosage.

Quantitative comparison of drug responses

We compared U0126- and wortmannin-treated profiles. Since they inhibit parallel but distinct pathways, one would expect distinct cellular responses. indicates that the promoter activities emphasize their differences (Pearson's correlation coefficient ~0.60), which are more evident than in gene expression profiles determined by microarray (Pearson's correlation coefficient ~ 0.92). Among the promoters significantly affected by the drugs, ~15% of them are commonly altered (i.e., suppressed or activated), whereas ~84% are affected by only a single drug (). Interestingly, a promoter, designated as chr7_-_74867620 (unannotated promoter), was affected in opposite ways (see Supplementary Table S2 online), even though microarray analysis did not detect significant expression changes at our working drug dosage. Overall, the promoter activity–based cellular response profiling approach successfully identified a distinct set of promoters downstream of individual pathways (~80% of the promoters), as well as a set of commonly affected promoters.

Quantitative explanation of the EGFR inhibitory effect by downstream inhibitors

We asked whether the profile for gefitinib could be explained by using the profiles for U0126 and wortmannin. As EGFR (a gefitinib target) primarily controls both Ras-ERK (target of U0126) and PI3K-Akt (target of wortmannin) pathways (), we hypothesized that EGFR inhibition can be explained by individual inhibition of these two pathways at the transcriptome level. To test this hypothesis, we used multiple linear regression analysis by taking gefitinib profiles as response variables and U0126 and wortmannin profiles as explanatory variables, where only promoters significantly affected by gefitinib are considered. Pearson's correlation coefficient of the regression model to the gefitinib profile, 0.749, is higher than the correlations of either U0126 or wortmannin profiles to the gefitinib profile (0.668 and 0.659, respectively; and Supplementary Figure S6 online). One could speculate that the correlation 0.749 can be achieved by combinatorial regression of any profiles, and thus we tested this possibility by performing the same regression analysis 100,000 times with random permutations of explanatory variables, that is, either U0126 or wortmannin profiles (see Supplementary Method 5.2 online). The correlation coefficients of the result of permutations (green and blue histograms, respectively; see Supplementary Figure S6 online) did not become higher than the one of the regression model based on U0126 and wortmannin with the gefitinib profiles, which indicates P value for the occurrence of such high correlation coefficient is <1 × 10−5. This result cannot be obtained with microarray analysis because the number of genes whose expression changes at equivalent thresholds are reliably measured is too small at an equivalent threshold (see Supplementary Figure S7 online). Even if we take a different threshold, by considering the result that the microarray profiles of U0126 and wortmannin responses are very close each other as described earlier (correlation coefficients 0.92; ), their combination will remain very close to their original profiles rather than to the gefitinib response.

We focused on the promoters that are roughly modeled by using both U0126 and wortmannin profiles. Of the 157 promoters significantly affected by gefitinib treatment, the level of expression changes in 133 (84.7%) promoters is consistent with the regression model (the difference of log2 fold changes in the gefitinib profile to control and the regression model are within 0.2-fold as shown in , and Supplementary Table S3 online for individual classifications). The remaining 24 promoters (15.3%) can be interpreted by postulating the existence of minor pathways affected by gefitinib only. For example, one of them may be a STAT pathway. It has been reported that two transcription factors, STAT1 and STAT2, are perturbed by gefitinib in a manner independent of the Ras-ERK and PI3K-Akt pathways.[28] Two of their target genes identified by ChIP experiments,[29] TNFR and KRT23, are included in the 24 promoters (). When focusing on the 133 promoters in which the regression model is consistent with the gefitinib profile, 50 (37.6%) are altered either by U0126 or wortmannin, and 72 (54.1%) are unaffected by either treatment (see Supplementary Table S3 online and ). This finding indicates that both U0126 and wortmannin profiles contribute independently to the regression model for majority of the promoters and that even subtle changes identified here reflect the actual perturbation at the transcriptome level.

Discussion

This is the first report on promoter activity profiling of drug responses using nonamplified CAGE libraries followed by single-molecule sequencing. Promoter activity profiles in response to the inhibition of ERK and PI3K kinases by U0126 and wortmannin, respectively, in human breast cancer MCF-7 cells are highly reproducible and consistent with previous studies. The sensitive profiling achieved by single-molecule sequencing allowed us to distinguish modest quantitative changes in the transcriptome even at submaximal drug dosage, which is very difficult to achieve with microarray-based profiling. Furthermore, the quantitative profiling allowed us to propose that both the inhibitory profiles of U0126 and wortmannin were constitutive components of the transcriptome profiles obtained by inhibition of the EGFR kinase, which is located upstream of both ERK and PI3K pathways. Although there are still some unexplained effects when using those two profiles as explanatory variables, our regression model successfully explained the changes in majority promoters significantly altered by gefitinib (, and Supplementary Table S3 online), and its statistical significance was confirmed by permutation test (see Supplementary Figure S6 online). Lange et al. explained protein phosphorylation in the same pathway by nonlinear model, which could explain the unexplained effect by our linear model presented here complementarily.[30] Nevertheless, the results here demonstrated the utility of sensitive and quantitative promoter activity profiling in elucidation of cellular responses to the drug.

Quantitative transcriptome analysis is potentially widely applicable in determining target proteins and action mechanisms of uncharacterized compounds. To further verify the applicability of this approach, additional analysis should be performed for different kinds of drugs with known action mechanisms, for example, herceptin (anti-ErbB2 receptor antibody), tamoxifen (estrogen receptor antagonist), cisplatin (DNA synthesis inhibitor), 5-fluorouracil, methotorexate (metabolome inhibitor), or docetaxicel (microtubule assembly inhibitor), which have been widely used for cancer therapies. Nonetheless, this study paves the way for quantitative analysis of drug responses at the promoter level, and moreover, it is potentially applicable for the evaluation of combinatorial or serial drug treatment in a clinical setting.

Methods

Cell culture and RNA preparation. The MCF-7 human breast cancer cell line was obtained from the American Type Culture Collection and maintained in Dulbecco's modified Eagle's medium (Invitrogen, Carlsbad, CA) supplemented with 10% fetal bovine serum. Before drug treatment, the medium was switched to Dulbecco's modified Eagle's medium supplemented with 2% fetal bovine serum and the cultures were maintained overnight. Gefitinib (ZD1839)[31] (1 μM, a generous gift from AstraZeneca, London, UK), U0126[32] (500 nM, Calbiochem, San Diego, CA), or wortmannin[33] (10 nM, Nacalai Tesque, Kyoto, Japan) was added and incubated with the cells for 6 h. Cells treated with dimethyl sulfoxide were used as the no drug treatment control. After the treatment, the cells were washed with phosphate-buffered saline twice, and RNA was isolated using the miRNeasy Mini kit (QIAGEN, Hilden, Germany) and analyzed by Bioanalyzer (Agilent technology, Santa Clara, CA).

Western blot. We confirmed the drug effect on the target proteins at the phosphoprotein level and defined the optimal drug concentrations so as to avoid a saturating effect. Cell lysates clarified by centrifugation were used in western blotting as described previously.[34] ERK1/2, phospho-ERK1/2(T202/Y204), Akt, phospho-Akt(S473), phospho-EGFR(Y1068) antibodies were purchased from Cell Signaling Technology (Danvers, MA), and anti-EGFR antibody was purchased from Fitzgerald Industries International (Acton, MA).

Production and analysis of single-molecule CAGE data. CAGE libraries for single-molecule sequencing were constructed as described previously[23] and sequenced on HeliScope. Sequenced reads were filtered and mapped to the February 2009 human genome assembly in UCSC (hg19 assembly). We grouped mapped reads overlapping within a few base pairs on the same strand into a single cluster as a promoter, and quantified promoter activities by counting the number of integrated reads in each cluster. Differential expression analysis between the control and the drug-treated samples was performed by using the edgeR package of Bioconductor in the R statistical language.[35] False-discovery rate <2% was used as a criterion for the differential expression against control experiments (see Supplementary Methods online).

Microarray analysis. We performed microarray analysis for establishing a baseline of the drug effects caused by our working doses, which was used for comparison with the drug effects of single-molecule CAGE profiles. Five hundred nanograms of total RNA was subjected to the Sentrix Human-6 Expression BeadChip (Illumina) according to the standard Illumina protocols. Variance-stabilizing normalization of Illumina data and B-statistics calculations were carried out using the lumi and limma packages of Bioconductor in the R statistical language.[36,37,38] B-statistics >4 was used as criterion for the differential expression against control (see Supplementary Methods online).

Linear regression model of the effect of EGFR inhibition. As a result of the differential analysis of gefitinib treatment against the control condition, we identified 157 promoters that were significantly altered. We then modeled the expression changes of these altered promoters by gefitinib treatment based on U0126- and wortmannin-response profiles. We defined the fold change of promoter activities by drug treatment (Fp,s) as follows:

where Ep represents the expression of promoter p on the control and e represents the expression in sample s. We used a linear regression analysis to model the response of promoter activities after gefitinib treatment (Fp,gefitinib). The formula of the linear regression is defined with error variable ep as follows:

and optimized the parameters (β1 and β2) by using the least squares approach, 0.7326 and 0.7391, respectively.

All of these analyses totally depend on how to quantify promoter expressions, while we constructed promoters operationally as described in Supplementary Method 3.1 online. One could hypothesize that the regression produces totally different results when we used different definitions of promoter. To examine this possibility, we estimated parameters (β1 and β2) with several definitions of promoter and found that the parameter values estimated above are representative, that is, within 2.5 percentile (0.5690 and 0.5930 for β1 and β2) and 97.5 percentile (0.7892 and 0.8713, respectively) in the promoter sets (see Supplementary Method 5.1 online for details).

Authors contributions

M.O.-H., H.K., and H.S. designed the research. M.O.-H. and H.S. performed the research. K.K. analyzed the data. K.K., M.O.-H., H.K., and H.S. wrote the manuscript.

Conflict of interest

The authors declared no conflict of interest.

Study Highlights