Novel method development and application
In this study, we have built upon a previous methodology using the NCI60 panel to identify genetic markers associated with drug response: namely gemcitabine and selenium by our group and perifosine by Zhang et al [9, 10, 21]. Instead of investigating specific candidate biological processes/pathways, here we used a genome wide approach with an aim of identifying genomic variants not originally thought to be related to the drug's mechanism. Previous studies on paclitaxel have investigated the pharmacogenetics of specific genes associated with paclitaxel response, genes associated with paclitaxel metabolism, and cDNA, protein activation status and siRNA changes associated with response [22–26]. However, this is the first genome-wide SNP association study done for paclitaxel response. Furthermore, the use of cell lines in these types of GWAS studies allowed access to high throughput genomic and transcriptomic data, along with robust drug screening protocols allowing for more accurate results and also allowed for greater control of non-genetic factors such as treatment compliance and administration time [27, 28].
Relevance to paclitaxel dosing in clinical settings
The GI50 response data was extrapolated from a range of doses which started at a 10-6 M dose and this was used as the basis for our case-control analysis. Although this dose was above therapeutic concentrations, the actual drug concentrations determined to cause 50% growth inhibition (Additional File 1) ranged from 0.104 nM to 7.96 μM. The minimum is below the therapeutic concentration range's minimum of 100 nM, while our maximum is slightly below the maximum therapeutic concentration 10 μM . This supports the clinical relevance of the GI50 values utilized in this study.
Impact of significant variants on respective genes and gene expression
Using FastSNP, four of the SNPs were found as intronic enhancers of their respective genes (SGCD, DCT, GRIK1, and CFTR). This suggests that the SNP variant does not alter the protein product significantly, but may alter transcription factor binding leading to altered expression patterns . In all cases, these genes were found to have significantly increased mRNA expression in sensitive cell lines for all probes available, supporting their predicted effect by FastSNP. The differential mRNA expression shown for SNTG1 suggests that its intronic variant may also serve as an intronic enhancer. Since only half of our protein-coding genes showed significant differences in mRNA expression, this emphasizes the importance to investigate genetic changes as the remaining five genes (ROBO1, LPHN2, PTPRD, BTBD12 and ZNF607) would not have been detected solely based on mRNA expression data. As well, these suggests that mechanisms other than changes in that specific gene's expression may help to explain how our identified SNPs can alter paclitaxel response. For the remaining protein-coding genes with intronic SNP variants: ROBO1, LPHN2 and PTPRD, it remains unclear how their variant may alter their gene product. However, SNPs located in intronic and intergenic regions can regulate the expression of many other genes (located on the same chromosome or different chromosomes); serving as master regulators and this may be how our intergenic and remaining intronic variants can be altering paclitaxel response . In the case of BTBD12 and ZNF607, their variants each caused conservative missense mutations with predicted benign effects on protein structure. These variants would not cause major changes in protein or cellular function but may still be significant enough to alter paclitaxel response. As well, the potential alternative splicing of ZNF607 can lead to further changes in the protein product, altering drug response. Alternatively, these SNPs may be in linkage disequilibrium with the real genetic factor that has a biological consequence directly or indirectly contributing to the paclitaxel response. Although FastSNP did accurately predict the potential effect for the 4 SNPs serving as intronic enhancers along with the intronic SNPs in ROBO1 and LPHN2 as not having any effects, these predicted changes are only suggestive at this stage and will need further validation using; for example, in vitro cell based functional assay systems. In addition, since the mRNA expression data from BioGPS served as adjunctive support for our findings, this data will need further validation by Real-time PCR or Northern blot analysis.
Haplotype blocks associated with paclitaxel response
In addition to individual SNP markers, we searched for haplotypes that were more significantly and strongly associated with paclitaxel response and hence can be more informative as response predictors. A total of five protein-coding genes had such haplotypes. In each of these haplotypes, our original SNP marker was found as part of the haplotype. Many of the remaining SNPs in the haplotypes were not found significantly associated with paclitaxel response when working alone and would have remained undetected unless they were working in concert with the original SNP marker. This supports the synergistic effect that multiple SNPs can exert on a gene which has been found to occur in processes such as cancer metastasis . This synergism can occur among SNPs in the same gene or between different genes [32–34]. The location of the SNPs in these haplotypes may also help to identify different gene regions or functional domains that can directly impact drug response and provide more insight into the drug's mechanism. Furthermore, one of our identified SNPs at 9p23 (marker number: 2836880) was found as part of the identified PTPRD haplotype. This suggests that by performing haplotype analyses, SNPs identified in the GWAS that were originally unclear on their role in paclitaxel sensitivity, may be found related to part of the gene through haplotypes.
Impact of genes on drug response
Our GWAS identified 10 genes containing SNP variants associated with paclitaxel response. Due to the in vitro nature of our study, it identifies only genes playing a direct role in the cellular response to paclitaxel, but do not reveal any variation potentially associated with paclitaxel's pharmacokinetics. Hence, the SNPs such as those previously identified in the drug metabolizing enzymes: CYP2C8, 3A4/5 and 2C19 and drug transport proteins: p-glycoprotein (MDR1) and SLCO1B3 would not be identified here [3, 22, 25]. Nevertheless, our analysis suggested that the genes identified here fall into two groups: those altering paclitaxel response through the β-catenin and p53 transcriptional regulatory axis, and those altering paclitaxel response through microtubule interactions.
1. β-catenin and p53 axis and paclitaxel response
Mutations of p53 and β-catenin have commonly been found associated with cancer risk and progression and in the case of p53, also response to chemotherapy [35, 36]. β-catenin plays a role in cell growth/morphogenesis, cell polarity and cell adhesion and induces expression of p53 and related candidates, which in turn regulates β-catenin's transcriptional activity [37–39]. p53 detects DNA damage and its activation can lead to cell cycle arrest and/or apoptosis [38, 39]. p53 is directly associated with paclitaxel chemosensitivity and paclitaxel increases p53 expression, phosphorylation status and its nuclear accumulation leading to the arrest of cell growth and apoptosis [4, 40–42].
Two of our identified genes PTPRD and BTBD12 interact with these two genes or their related pathways and may impact paclitaxel response through them. p53 decreases expression of protein tyrosine phosphatase receptor type A (PTPRA) which interacts with our candidate, PTPRD [43, 44]. PTPRA and PTPRD are receptor tyrosine phosphatases, and PTPRD can dephosphorylate v-src sarcoma viral oncogene homolog (SRC), altering v-raf-1 murine leukemia viral oncogene homolog 1 (RAF-1) activation, a major pathway for paclitaxel induced apoptosis [2, 6, 45, 46]. BTBD12 is a scaffold protein that helps to regulate various DNA repair enzymes including those used for double stranded breaks and cross-links and also interacts with enzymes for cell-cycle control, homologous recombination and replication forks [47–50]. In addition to arresting cell division, paclitaxel can also induce DNA damage, while repressing various DNA repair genes leading to cytotoxicity . Thus, the missense substitution caused by BTBD12's variant could reduce its DNA repair ability leading to increased paclitaxel susceptibility or alter paclitaxel's ability to arrest mitosis.
2. Microtubule interactions and paclitaxel response
Microtubules play an important role in paclitaxel response. Paclitaxel binds to β-tubulin when it is present as assembled tubulin and causes cell-cycle arrest, and eventually apoptosis . Previous studies have found that point mutations in tubulin, which weaken drug-tubulin interactions and altered expression of a microtubule stabilizing protein, microtubule-associated protein 4 (MAP4), both alter paclitaxel response in cells [2, 4].
Several of our identified genes may also affect paclitaxel response through changes in microtubules dynamics or stability. ROBO1 is a neuronal receptor used in neurodevelopment and when activated by its ligand Slit, directly interacts with the beta 4 tubulin subunit (TUBB4) of microtubules to regulate microtubule dynamics [52–55]. CFTR is a chloride transporter channel and is required for proper microtubule distribution, as studies have shown that cells having mutant CFTR (DeltaF508) show microtubular disorganization . ZNF607 is a zinc finger protein involved in transcriptional regulation and it interacts with Ataxin 1 (ATXN1), which interacts with microtubules to alter cell morphology during neuronal development [54, 57, 58].
As well, many of our identified genes use microtubules for their own transport or transport of other proteins and therefore these genes can also cause changes in microtubule dynamics, potentially altering paclitaxel response. GRIK1 is a kainite receptor used for interneuronal communication and uses microtubules to get transported to the surface of dendrites [59, 60]. Additionally, GRIK1 can indirectly affect cellular binding to paclitaxel through co-regulation of expression by Atrophin 1 (ATN1), which also regulates Stathmin 1 (STMN1) expression and STMN1 reduces cellular binding of paclitaxel [61–63]. LPHN2 interacts with microtubules when it binds to a ligand, α-Latrotoxin in order to cause neurotransmitter release . DCT belongs to a group of enzymes used in the biosynthesis of melanin and this process uses microtubules to transport melanosomes containing melanin, to the cell surface [65, 66]. However, DCT is also a marker for melanoma and the observed results for this gene may be caused by the melanoma cell lines all belonging to the sensitive cell line category .
The last two of our genes, SGCD and SNTG1 which represent members of their respective families, sarcoglycans and syntrophins form a part of a larger complex called the dystrophin associated glycoprotein complex (DAGC). The DAGC connects the cytoskeleton to the extracellular matrix in muscle . The beta isoform of sarcoglycan, SGCB along with SGCD have been found to co-localize with tubulin during the M phase of the cell cycle, and this corresponds to the cell cycle stage when paclitaxel interacts with tubulin [4, 25, 69]. Members of the syntrophin family, including β2-syntrophin have been found to associate with purified microtubules through a PDZ domain which is shared by SNTG1 [70, 71]. The association of these genes with microtubules suggests that they may play a role in regulating microtubule dynamics; thereby potentially altering paclitaxel response.
The NCI60 panel comprises of a limited number of human cancer-derived cell lines, and as such, there may be confounding by either race or tumor type in our identification of pharmacogenomic markers. For instance, DCT is also a marker for melanoma and is melanocyte specific [65, 67]. The melanoma cell lines in this case were all found in the sensitive group and this may serve as a confounding variable. Validation of our results in non-cell line data and in other tumor types are therefore necessary prior to clinical translation.
Working through the p53/β-catenin axis and microtubule interactions, our identified SNP variants can alter their respective gene products and in turn modulate the cellular response to paclitaxel causing the observed heterogeneity in both cellular and clinical responses.
Requirement for an integrative approach in pharmacogenomics
Traditionally, the discovery of markers associated with drug response or prognosis of complex disease risk has depended upon investigating changes in the gene expression, in particular mRNA expression levels [72, 73]. However, if this was applied here, only five of our genes would have been identified. Thus, it is important to use both genomic and expression data sets to complement each other as genomic variants which ultimately alter gene functions other than expression may not be detected . Furthermore, it has been suggested that the discovery of such markers should use an integrated functional genomics analysis where information from genomic (i.e., mutation, copy number variation, SNPs, DNA methylation), expression and functional experiments should be combined to find genetic marker profiles . This analysis could also include haplotype analysis which can reveal markers that are more meaningfully linked with traits, as done here, but also include analysis such as proteomics, epistatic changes and reverse-phase protein arrays (RPPA) [72, 73, 75].
A similar type of integrated analysis has been performed for paclitaxel by Park et al, where RPPA profiling along with mutation analysis and siRNA screening analysis was used to find and assess markers associated with paclitaxel sensitivity . However, our results differed from those by Park et al, in part due to the different approaches used to identify markers associated with paclitaxel sensitivity. Our study focused on investigating a genetic association through categorization (sensitive or resistant) while Park et al, correlated sensitivity with protein expression level. As seen in our mRNA expression analysis, not all identified SNPs belonging to protein-coding genes showed changes in their respective gene's mRNA expression and this can hold true for protein expression as well. Also, although not investigated here, it is likely that other genes, not identified in our GWAS, show significant differences in expression levels between sensitive and resistant cell lines, despite not having any SNPs associated with drug response. As discussed previously, one possible mechanism for this could be through SNPs in other genes or intergenic SNPs altering the expression of genes on the same or different chromosomes and reinforces the need to complement expression data with genomic analysis . Hence more integrative analysis approaches, integrating data from genetic (SNPs and copy number variations), mRNA expression, protein expression studies should be utilized to better understand the role of genetic variants in drug response.