Using transcriptome sequencing to identify mechanisms of drug action and resistance.

Determining mechanisms of drug action in human cells remains a major challenge. Here we describe an approach in which multiple-drug-resistant clones are isolated and transcriptome sequencing is used to find mutations in each clone. Further analysis of mutations common to more than one clone can identify a drug's physiological target and indirect resistance mechanisms, as indicated by our proof-of-concept studies of the cytotoxic anticancer drugs BI 2536 and bortezomib.

Currently, approaches to analyze how a drug works fall into two broad categories. First, several strategies rely on model organisms that are compatible with genetic manipulations[1-3]. However, many drugs are inactive in these organisms, possibly due to multi-drug resistance mechanisms and target divergence[4]. Second, affinity-based methods are used to identify proteins that bind the drug[5, 6]. These approaches are effective when the drug is potent and the targets are reasonably abundant in vivo. However, establishing that a drug-binding protein is the physiological target typically depends on correlations between binding/activity assays in vitro and protein knockdown phenotypes. These correlations can be misleading due to differences between chemical inhibition, which can be acute, and the protein knockdown phenotypes, which can be indirect or from cumulative effects[7]. Currently, we lack an approach that addresses these limitations and can determine mechanisms of drug action in human cells.

The ‘gold standard’ in identifying a drug’s target is achieved when two criteria are met. First, resistance to a drug in a physiological context occurs through mutation in the target protein. Second, this mutation should suppress inhibition of the target’s activity by the drug. This standard is met for a few drugs (e.g. Gleevec), for which focused analyses of the expected target led to such mutations[8, 9]. However, the human genome’s large size and complexity has limited unbiased analyses of mutations conferring drug resistance.

We sought to develop an approach that could achieve ‘gold standard’ validation of a drug’s target and be applied to relevant human cell types (e.g. cancer cells). We reasoned that deep sequencing and bioinformatics could be used to screen the genome of drug resistant clones isolated from cultured human cells for resistance mutations. Moreover, we hypothesized that analyzing mutations common to multiple clones would focus our analysis on a handful of potential drug targets. When combined with biochemical analyses, these data could lead to ‘gold standard’ validation of the drug’s physiological target. As a proof-of-concept we analyzed the mechanisms that confer resistance to drugs whose targets are known.

BI 2536 (Fig. 1a) is currently in clinical trials and inhibits Polo-like kinase 1 (PLK1), a major cell cycle regulator[10]. To isolate resistant clones we used a human colon cancer cell line, HCT-116, which is deficient in mismatch repair and has low expression of multi-drug resistant (MDR) pumps[11]. Therefore, HCT-116 behaves like a mutagenized cell line, facilitating the rapid identification of mutations that confer drug resistance[12, 13]. From fifteen BI 2536-resistant clones (isolated from selections with 10 nM BI 2536, LD50: 3.9 ± 2.8 nM), we selected six for transcriptome sequencing (also called RNA-seq). The LD50s for BI 2536 were 3-9 fold higher in these clones than in the parental cell line (Fig. 1b, Supplementary Results, Supplementary Table 1). In parallel, the parental cell line was also processed for transcriptome sequencing.

Figure 1

Characterization of BI 2536-resistant clones

(a) Structure of PLK1 inhibitor BI 2536. (b) Median lethal dose (LD50s) measured for the parental cell line and two drug-resistant clones (LD50s: 3.9 ± 2.8 nM (parental), 33.5 ± 2.6 nM (Clone B), 14.2 ± 7.8 nM (Clone D); n = 3, mean ± sem). (c) Graph-based analysis of similarities (0 = low similarity, 2 = high similarity) between BI 2536-resistant clones allows clones D, E, and F to be grouped together (group 4). (d) Proliferation assay showing the effects of BI 2536 exposure on hTERT-RPE1 cells, transfected with GFP-PLK1 wt, GFP-PLK1 G63S, or GFP-PLK1 R136G. Median lethal dose (LD50s) measured for each transfected cell line (LD50s: 44 ± 5 nM (GFP-PLK1 wt), 83 ± 9 nM (GFP-PLK1 G63S), 76 ± 8 nM (GFP-PLK1 R136G); n = 6, mean ± sem, p < 0.01 for both, two-tailed paired t-test). (e) Proliferation assay showing the effect of 20 nM taxol on HCT-116 parental cells and clones A, B, and C, normalized to untreated cells (n = 3, mean ± s.d.).

The sequencing data were analyzed to identify single nucleotide variations and short insertions/deletions (indels) in each clone, and only those present in coding sequences were further considered. We found 6-14 single nucleotide variations that were significantly increased (with a 0.5% false discovery rate) in the BI 2536-resistant clones, compared to the parental cell population (Supplementary Tables 2-7; no indels were found). Groups of similar clones were identified by analyzing single nucleotide variations using a clustering approach (see Supplementary Methods). Among the six clones, A, B, and C were independent (groups 1, 2, and 3, respectively), while D, E, and F formed a single group (group 4, Fig. 1c). We further focused our analysis on genes that were mutated in more than one BI 2536-resistant group. While there were no genes common to all groups, PLK1 was the only gene mutated in more than one group (Supplementary Table 8).

Two distinct PLK1 mutations were identified, G63S and R136G, both of which map to the binding site of BI 2536 in the crystal structure of PLK1 (Supplementary Fig. 1)[14]. To examine whether these mutations are sufficient for conferring resistance to BI 2536, we generated stable cell lines expressing wildtype PLK1 and each mutation as a GFP-tagged full-length PLK1 construct. Independent cell lines were used as secondary mutations were less likely to be common between these lines and HCT-116. The R136G mutation suppressed BI 2536 toxicity in hTERT-RPE1 (Fig. 1d), and the G63S mutation suppressed toxicity in hTERT-RPE1 (Fig. 1d) and in HeLa cells (Supplementary Fig. 2). These data suggest that PLK1 is the major physiological target of BI 2536. Published data show that the R136G mutation has reduced sensitivity to BI 2536 in vitro[15]. Therefore, these data provide the ‘gold standard’ proof that PLK1 is the target of BI 2536.

BI 2536-resistant group 4 (representing clones D, E, and F) did not have a mutation in the PLK1 gene, suggesting a different resistance mechanism. An unbiased survey of all transcript levels (Supplementary Data 1 and 2), which can be determined from the transcriptome sequencing data, revealed that ABCB1 (P-gp, a drug efflux transporter) was among the most highly over-expressed mRNAs in group 4 clones (Supplementary Fig. 3). Consistent with the hypothesis that the intracellular concentration of BI 2536 was lowered by increased drug efflux, group 4 clones also had reduced sensitivity to the chemically-unrelated drug taxol, a compound known to be transported by the ABCB1 pump[16] (Fig. 1e). While further experiments are needed to confirm this resistance mechanism, our data suggest that our approach can lead to testable hypotheses for indirect mechanisms of drug resistance.

We next examined if our approach could be applied to another drug. Bortezomib (Fig. 2a) inhibits the proteasome by targeting its subunit PSMB5 and is used clinically to treat multiple myeloma and mantle cell lymphoma[17]. Nineteen clones were isolated from HCT-116 cells grown in the presence of bortezomib (8 -12 nM, LD50: 6.3 ± 0.9 nM). Five clones, with reduced bortezomib sensitivity (LD50s 2.4 – 6.5 fold higher, Fig. 2b and Supplementary Table 9), were processed for transcriptome sequencing. 15-28 single nucleotide variations in each clone were identified (Supplementary Tables 10-14). Clustering analysis (Fig. 2c) grouped clones A, B, C, and D together (group 1) while clone E was independent (group 2). Five genes were mutated in both bortezomib-resistant groups, and the only gene with two distinct mutations (M104V and A108T) was PSMB5, the known target of bortezomib (Supplementary Table 15). If the drug target was unknown, all five genes would have to be examined as potential targets. However, the existence of two distinct resistance mutations would make PSMB5 the highest priority for further analysis.

Figure 2

Characterization of bortezomib-resistant clones

(a) Structure of proteasome inhibitor bortezomib. (b) Median lethal dose (LD50s) measured for the parental cell line and two drug-resistant clones (LD50s: 6.3 ± 0.9 nM (parental), 19.4 ± 6.7 nM (Clone A), 41.5 ± 9.5 nM (Clone E); n = 3, mean ± sem). (c) Graph-based analysis of similarities (0 = low similarity, 2 = high similarity) between bortezomib-resistant clones groups clones A, B, C, and D (group 1). (d) Proliferation assay showing the response of hTERT-RPE1 cells, transfected with GFP-PSMB5 wt, GFP-PSMB5 M104V, or GFP-PSMB5 A108T, to treatment with bortezomib. Median lethal dose (LD50s) measured for each transfected cell line (LD50s: 5.55 ± 0.09 nM (GFP-PSMB5 wt), 28.7 ± 2.8 nM (GFP-PSMB5 M104V, p < 0.05, two-tailed paired t-test), 19.2 ± 2.8 nM (GFP-PSMB5 A108T, p < 0.05, two-tailed paired t-test); n = 3, mean ± sem). (e) Schematic highlights the key steps of the approach: Selecting and expanding drug-resistant clones from a heterogeneous parental population. Massively parallel sequencing of mRNA from multiple drug-resistant clones and parental (untreated) cells. Bioinformatics analyses to identify genes that are mutated with high frequency. A subset of sequencing reads for a BI 2536-resistant clone (clone A) and the HCT-116 parental cells are shown. The ~50% mutation frequency indicates heterozygosity.

We found that expression of GFP-tagged PSMB5, carrying either of these mutations, suppressed bortezomib sensitivity in an independent cell line (hTERT-RPE1, Fig. 2d). This is consistent with previous reports that A108T confers bortezomib-resistance[18, 19]. As both mutations in PSMB5 map to the drug’s binding site, we hypothesize that they directly suppress drug interactions (Supplementary Fig. 4)[20]. While additional biochemical tests are needed to examine this further, our data indicate that our strategy can efficiently lead to resistance mechanisms that include a drug’s direct target.

Our analysis thus far indicates that for our approach to be effective, resistance via mutations in a drug’s direct target must occur at high frequency in drug-resistant clones. To examine this, we sequenced the PLK1 gene in each of the nine BI 2536-resistant clones that we had not processed by RNA-seq. PLK1 was mutated in ~45% of these clones (4 of 9, Supplementary Table 16). We next analyzed two kinesin-5 inhibitors, S-Trityl-L-cysteine (STLC), which is known to be selective[21], and 4-(2-(1-phenylcyclopropyl)thiazol-4-yl)pyridine (PCTP), which has been shown to inhibit other related motor proteins in vitro[21, 22]. Kinesin-5 mutations were found in ~30% of the STLC-resistant clones (4 of 14, Supplementary Fig. 5) and in ~15% of PCTP resistant clones (3 of 22, Supplementary Fig. 6). These data indicate that resistance in a drug’s direct target occurs at high frequency when a drug has one major physiological target (as is the case for STLC, BI 2536 and bortezomib). When a drug has multiple targets (e.g. PCTP) it is likely that resistance in a single target will be less frequent and our approach may be limited as a greater number of clones would have to be sequenced to identify mutations common in more than one clone.

In summary, we have developed a method to identify the target of drugs in human cells by examining resistance mechanisms. Our strategy (Fig. 2e) involves isolating multiple drug resistant clones from genetically heterogeneous human cells. Clones with multi-drug resistance can be excluded by testing for reduced sensitivity to unrelated compounds (i.e. taxol). The remaining clones are processed for transcriptome sequencing, along with the parental (untreated) cell population. Bioinformatics is used to find genes mutated in more than one independent clone. These genes are prioritized for further biochemical and cell biological analyses to identify the drug’s direct target and indirect resistance mechanisms.

Our strategy has advantages over other approaches to identify a drug’s target. First, unlike many other target identification approaches, our method does not rely on chemical modifications of the drug of interest. This can be important when small changes in a drug’s chemical composition can alter its mechanism of action[23]. Second, our approach can enable cell-type specific analyses, which would be particularly useful if a drug is toxic in specific tissues. While our approach is currently limited to analyzing cytotoxic drugs, it should be possible to extend its use to non-toxic drugs by using phenotypic or reporter-based read-outs to select clones with reduced drug sensitivity (e.g. change in fluorescence). Furthermore, as our approach focuses on discovering single nucleotide variations and insertions/deletions, it does not report on all potential mechanisms of resistance. This limitation can be addressed by combining our strategy with other genomic methods, such as exome capture, full genome sequencing, and bisulfate conversion followed by sequencing (to detect DNA methylation). Our approach has the potential to reveal the physiological on-targets of a drug in disease cells, unintended off-targets in healthy cells, and can reveal cellular mechanisms of drug resistance. These findings can impact chemical modifications of drugs to improve efficacy and limit toxicity. When unanticipated drug targets are found new uses of the drugs may also be suggested.