Integration of Cell Line and Clinical Trial Genome-Wide Analyses Supports a Polygenic Architecture of Paclitaxel-Induced Sensory Peripheral Neuropathy

This article is featured in Highlights of This Issue, p. 309

Paclitaxel is a tubulin-targeting agent, widely used in the treatment of malignant disease, including ovarian, breast, lung, and head and neck cancers. Its long-term use is often limited by sensory peripheral neuropathy, although the mechanism of this toxicity is poorly understood. In one recent large study of more than 1,500 patients with breast cancer, severe (grade 3) sensory peripheral neuropathy occurred in 4% of patients treated with 4 cycles and 10% of patients treated with 6 cycles of single-agent paclitaxel (1). Currently, genetic prediction of which patients with cancer may experience severe side effects induced by paclitaxel treatment is not possible (2), but several preliminary genetic associations have been made (3–8). If patients likely to experience such toxicities could be identified before beginning a paclitaxel regimen, patient care might be improved by implementing a reduced dose or an alternative treatment.

Because accruing large patient cohorts receiving the same drug regimen for discovery genome-wide and replication studies in oncology is challenging, several groups of investigators have used lymphoblastoid cell lines (LCL) as a discovery tool and for follow-up functional studies (9–13). LCLs are easy to experimentally manipulate and the genetic background and expression environment is known. The LCL model also permits functional validation studies of candidate markers and genes discovered in both preclinical and clinical studies (14, 15). However, a critical question is how well this cell-based model generates clinically relevant markers and genes associated with patient response to drug. Recently, a few chemotherapeutic response single-nucleotide polymorphisms (SNP) discovered in LCLs have been replicated in patient populations by associating with phenotypes such as tumor response and overall survival in patients receiving the same drug (16–19); however, these studies tested the individual variants most associated with the LCL phenotypes. We sought to understand to what extent the overall genetic architecture of patient response to chemotherapy can be captured by LCLs by investigating beyond just the top few signals. In contrast to previous studies that tested single SNPs, we use an enrichment method (20) to determine in a systematic manner whether top genome-wide association study (GWAS) SNPs for paclitaxel-induced sensory peripheral neuropathy in patients with breast cancer (3) are more likely to also be paclitaxel-induced cytotoxicity SNPs identified in LCLs than expected by chance.

In this study, we found that SNPs associated with patient paclitaxel-induced neuropathy are enriched for SNPs associated with paclitaxel-induced cytotoxicity in HapMap LCLs. This significant enrichment confirms that LCLs are a useful model in the study of a subset of shared genes involved in patient toxicity. The overlap SNPs are predominantly expression quantitative trait loci (eQTL) as defined previously (21), therefore supporting an enriched functional role for these significant SNPs. We show a functional role for one eQTL host gene (RFX2) in paclitaxel toxicity, using a cellular model of peripheral neuropathy. These results are consistent with the hypothesis that the cell-based models capture components of the underlying genetic architecture for paclitaxel-induced sensory peripheral neuropathy.

HapMap LCLs from a population with Northern and Western European ancestry from Utah (HAPMAPPT01, CEU, n = 77), a Yoruba population in Ibadan, Nigeria (HAPMAPPT03, YRI, n = 87), and an African-American population from the Southwest of the United States (HAPMAPPT07, ASW, n = 83) were treated with 12.5 nmol/L paclitaxel and cytotoxicity was determined using an AlamarBlue (Invitrogen) cellular growth inhibition assay as described (22). The cytotoxicity phenotype used in the LCL GWAS was mean percentage survival at 12.5 nmol/L paclitaxel determined from 6 replicates from 2 independent experiments. Percentage survival values for each cell line were log₂-transformed before statistical analysis to form an approximately normal distribution in each population.

A GWAS of paclitaxel-induced cytotoxicity was conducted on each of the 3 populations separately. Greater than 2 million SNPs from HapMap r27 [minor allele frequency (MAF) > 0.05 within the panel, no Mendelian errors and in Hardy–Weinberg equilibrium (P > 0.001)] were tested for association with paclitaxel cytotoxicity in each population, using the quantitative trait disequilibrium test total association model (23). To control for population structure in the admixed ASW population, local ancestry at each genotyped SNP locus was estimated using HAPMIX (24) and to increase genome coverage of the ASW, ungenotyped markers were imputed using BEAGLE (25) as previously described (26). Genomic control lambda (λ_GC) values (27) were calculated for the GWAS of each population. Studies with λ_GC values greater than 1 were corrected for residual inflation of the test statistic by dividing the observed test statistic at each SNP by the λ_GC (27), and then the corresponding P values were carried through the meta-analysis.

Using the software METAL, we combined SNP P values across the 3 population studies, taking into account a study-specific weight (sample size) and direction of effect (positive or negative β; ref. 28). This approach converted the direction of effect and P value observed in each study into a signed Z-score, such that very negative Z-scores indicate a small P value and an allele associated with higher drug sensitivity, whereas large positive Z-scores indicate a small P value and an allele associated with higher drug resistance. Z-scores for each SNP were combined across studies in a weighted sum, with weights proportional to the square-root of the sample size for each study (28).

Cancer and Leukemia Group B (CALGB) 40101 is a phase III trial comparing the efficacy of standard therapy cyclophosphamide and doxorubicin with single-agent paclitaxel as adjuvant therapy for breast cancer in women with 0 to 3 positive axillary lymph nodes. All study participants were enrolled in CALGB 40101 and gave their additional consent to participate in the pharmacogenetic companion study (CALGB 60202), which has been published (3). All patient research met state, federal, and Institutional Review Board guidelines. Germline DNA was isolated from 1,040 patients on the paclitaxel arm of CALGB 40101 and genotyped using the Illumina 610-Quad platform as described previously (3). Following quality control analysis, genotypes were available for 520,679 SNPs. Principal component (PC) analysis identified 855 genetic Europeans that were used in a GWAS of sensory peripheral neuropathy (3). A dose-to-event analysis was conducted, with an event defined as grade 2 or greater sensory peripheral neuropathy. The Cox score test, powered for additive genetic effects, was used to test these marginal associations. Only SNPs with MAFs more than 0.05 in the patient population and in Hardy–Weinberg equilibrium in the CEU (P > 0.001) were used in the LCL GWAS comparisons.

We conducted a permutation resampling analysis (29) to test for an enrichment of cytotoxicity-associated SNPs (LCLs) among the paclitaxel-induced sensory peripheral neuropathy-associated SNPs (patients). To this end, the patient outcomes (cumulative dose and event indicator vectors) were randomly shuffled while keeping the genotype data fixed to preserve linkage disequilibrium. On the basis of this permutation replicate, the standardized Cox score statistics were recalculated for all the SNPs. This process was conducted 1,000 times. For each of the 1,000 permutation replicates, the number of SNPs that had P < 0.05 in the patient data, P < 0.001 in the LCL data, and the same direction of effect (the same allele associated with increased neuropathy and increased cytotoxicity) was calculated. The overlap distribution from the permutations was compared with the observed SNP overlap to generate an empirical P value, calculated as the proportion of permutations in which the number of LCL/patient overlap SNPs exceeds the observed number. To test the robustness of our findings, we calculated an empirical P value across a range of inclusion thresholds from P < 0.001 to P < 0.1. We also tested for enrichment of patient SNPs among the LCL SNPs by generating 1,000 randomized SNP sets the same size and MAF distribution as the observed LCL data at a range of P value thresholds to calculate empirical P values. In addition to the paclitaxel LCL cytotoxicity data, we compared the patient sensory peripheral neuropathy data with LCL cytotoxicity GWAS data from capecitabine (30) and carboplatin (13) as negative controls.

To test for eQTL enrichment in the LCL, patient, and LCL/patient overlap SNPs, we generated 10,000 randomized SNP sets each of the same size as the observed set of LCL cytotoxicity (P < 0.001), patient neuropathy (P < 0.05), or LCL/patient overlap SNPs. The randomized SNP sets were matched on MAF distribution of the observed list and sampled (without replacement) from the set of SNPs on the Illumina 610-Quad platform, similar to the method of Gamazon and colleagues (31). We grouped the platform SNPs into discrete MAF bins of a width of 5%, from which the SNPs used in the simulations were selected. For each of the 10,000 sets, we determined the number of eQTLs (P < 10⁻⁴) and calculated an empirical P value for enrichment. The eQTLs were defined previously and are available in the SCAN database (21, 31).

First, we determined which of the LCL/patient overlap SNPs from the enrichment analysis were located in or near (within 2 kb) gene transcripts (dbSNP build 129, human genome assembly build 36). Eleven of 24 overlap SNPs were in or near genes and genotyping intensity plots for these SNPs in the patient data are available in Supplementary Fig. S1. Second, we determined which SNPs within genes were also eQTLs (31) and prioritized by which had the most target genes (P < 10⁻⁴). We also tested whether the expression of the eQTL target genes associated with paclitaxel-induced cytotoxicity (P < 0.05) using previously published exon array data (32). A general linear model was constructed between gene expression and paclitaxel-induced cytotoxicity with growth rate (33) and population as covariates. A Toeplitz covariance structure with 2 diagonal bands was used to allow for familial dependencies in the data as previously described (9).

Neuroscreen-1 (NS-1) rat pheochromacytoma cells (Cellomics Inc.) were maintained in NS-1 media (RPMI supplemented with 10% horse serum, 5% fetal calf serum and 1% l-glutamine). Cells were seeded at a density of 1 × 10⁵ cells/mL on collagen I–coated plates and induced to differentiate by adding 20 ng/mL nerve growth factor (NGF, BD Biosciences) to the media 24 hours before transfection. Cells for cytotoxicity assays were plated in 96-well collagen I–coated plates, whereas cells for expression quantification and neurite-outgrowth assays were plated in 6-well collagen I–coated plates. Pooled Rfx2 siRNA (25 nmol/L; Qiagen; S101639659, S101639666, S101639673, and S101639680) or nontargeting control siRNA (Qiagen; 1027292) was transiently transfected into the NS-1 cells using DharmaFECT Reagent #1 (Dharmacon). Quantitative reverse transcription PCR (qRT-PCR) was conducted for Rfx2 (Rn00501380_m1) and control gene Gapdh (4352338E) using TaqMan Gene Expression Assays (Applied Biosystems) 24 hours posttransfection in the neurite-outgrowth assays and 24, 48, 72, and 96 hours posttransfection in the cytotoxicity assays to assess Rfx2 knockdown in NS-1 cells. Expression of the potential Rfx2 target genes Cyp51 (Rn01526553_m1), Bach1 (Rn01477344_m1), and Cbara1 (Rn01644475_m1) was also measured by qRT-PCR at 24 hours post-siRfx2 transfection. Each qRT-PCR was run in duplicate and individual samples were run in triplicate on each plate. Percentage knockdown was calculated by dividing the relative Rfx2 expression levels in the siRfx2 sample by those in the nontargeting control sample.

Twenty-four hours following siRNA transfection, transfection media was removed from the NS-1 cells and 0, 12.5, or 100 nmol/L paclitaxel in NS-1 media (supplemented with 20 ng/mL NGF) was added to either the siRfx2 or nontargeting control cells. After 24 hours in the presence of paclitaxel, phase-contrast images (×10) of the cells were taken using an Axiovert 200M inverted widefield fluorescence microscope (Zeiss). At least 500 cells per treatment in 6 randomly chosen fields were imaged and the longest neurite per cell was measured using ImageJ (34) software. The entire experiment was carried out in duplicate and mean neurite lengths were normalized relative to the 0 nmol/L drug treatment for each siRNA. Because tracing neurite lengths is somewhat qualitative, 2 scientists independently measured neurite lengths and the second scientist was blinded to siRNA/drug treatment. The mean of each set of measurements between the 2 scientists was assessed for significance by 2-way ANOVA (factors: siRNA treatment and drug treatment) to determine if the siRfx2 affected neurite length upon paclitaxel treatment.

Twenty-four hours after siRNA transfection, transfection media was removed from the NS-1 cells and 0, 6.25, 12.5, 25, 50, or 100 nmol/L paclitaxel in NS-1 media (supplemented with 20 ng/mL NGF) in triplicate was added to either the siRfx2 or nontargeting control cells. After 72 hours of paclitaxel treatment, ATP levels were measured using the CellTiter-Glo assay (Promega) and percentage survival curves were generated. The entire experiment was done in duplicate and 2-way ANOVA was used to determine if the siRfx2 significantly affected overall cytotoxicity upon paclitaxel treatment.

We conducted a genome-wide meta-analysis (see Materials and Methods) to test common SNPs for association with paclitaxel-induced cytotoxicity in LCLs. We compared the results from this analysis with those from clinical trial CALGB 40101, a GWAS of paclitaxel-induced sensory peripheral neuropathy in patients with breast cancer (3). Neither study produced genome-wide significant results (α < 0.05) nor did the very top SNPs match between the 2 studies (Fig. 1). However, through a permutation resampling analysis of the CALGB patient data, we found that the top sensory peripheral neuropathy-associated SNPs (P < 0.05) are significantly enriched for SNPs associated with paclitaxel-induced cytotoxicity in LCLs (P < 0.001) with consistent allelic directions of effect (Fig. 2; empirical P = 0.007). The observed enrichment of 24 SNPs between the LCL and patient studies is likely paclitaxel-specific, due to the sensory peripheral neuropathy SNPs not being enriched for either capecitabine- or carboplatin-induced cytotoxicity SNPs, which were tested as negative controls (Fig. 2). Positional information and effect sizes of all 24 overlap SNPs in the LCL and patient data can be found in Supplementary Table S1. When the inclusion thresholds for overlap SNPs were relaxed and when the LCL SNPs were tested for enrichment of patient SNPs, the significant overlap was present at a range of P value thresholds from 0.001 to 0.1, showing the robustness of our findings (Supplementary Table S2).

We tested the top paclitaxel-induced LCL cytotoxicity SNPs (P < 0.001) and the top paclitaxel-induced patient sensory peripheral neuropathy SNPs (P < 0.05) for eQTL enrichment because these were the thresholds used in the primary overlap analysis. We compared the observed number of eQTLs at these thresholds to the number of eQTLs in 10,000 randomly selected MAF-matched SNP sets (for details, see Materials and Methods). Neither cytotoxicity-associated SNPs nor neuropathy-associated SNPs alone were enriched for eQTLs (Fig. 3). However, we found that the 24 paclitaxel LCL/patient overlap SNPs at these thresholds are enriched for eQTLs when compared with MAF-matched SNP sets (empirical P = 0.0447), potentially revealing an important role for this functional class in paclitaxel toxicity.

First, we determined that 11 of 24 overlap SNPs from the enrichment analysis were located in or near (within 2 kb) gene transcripts (Table 1). The relationship of these 11 SNPs with paclitaxel-induced sensory peripheral neuropathy in patients and LCL cytotoxicity is shown in Supplementary Fig. S2. Next, we determined which of these 11 SNPs within genes were also eQTLs (31). Of the 8 eQTLs, we determined which had the most potential target genes at an arbitrary threshold of P < 10⁻⁴. The SNP in RFX2 had 18 target genes, more than any other of the 8 eQTLs. In addition, we tested the expression of the target genes for association with paclitaxel-induced cytotoxicity adjusted for growth rate (see Materials and Methods). We found that expression of 3 of the RFX2 target genes associated with paclitaxel-induced cytotoxicity (Table 1 and Supplementary Table S3); therefore, we pursued evaluating RFX2 in a model of neuropathy.

We used neuroscreen (NS-1) cells, a subclone of the rat pheochromocytoma cell line PC-12 that has previously been used as a research model for chemotherapy-induced neuropathy (35, 36), to test Rfx2, the rat ortholog of RFX2, for functional involvement in paclitaxel response. Using siRNA, we decreased expression of Rfx2 resulting in increased sensitivity of NS-1 cells to paclitaxel, as measured by reduced neurite outgrowth and increased cytotoxicity (Fig. 4). The 3 RFX2 SNP target genes whose expression associated with paclitaxel-induced cytotoxicity in LCLs are CYP51A1, BACH1, and CBARA1 (Table 1; Fig. 5A–C; P < 0.05). We measured the expression of these 3 potential Rfx2 target genes upon knockdown of Rfx2 in NS-1 cells and found that the expression of 1 of 3 genes, Cyp51 (rat ortholog of CYP51A1), significantly decreased 24 hours posttransfection (P < 0.05), which is the expected direction of effect based on the LCL expression versus cytotoxicity data (Fig. 5D).

We conducted a GWAS of paclitaxel-induced cytotoxicity in LCLs and showed significant enrichment of the top cytotoxicity-associated SNPs in a clinical GWAS of paclitaxel-induced sensory peripheral neuropathy in patients with breast cancer. This robust enrichment shows that susceptibilities to increased cytotoxicity in LCLs and sensory peripheral neuropathy in patients with breast cancer likely have some genetic mechanisms in common and supports the role of LCLs as a preclinical model for paclitaxel toxicity studies. Furthermore, the top SNPs that overlap between the 2 studies were enriched for eQTLs. This eQTL enrichment indicates that SNPs associated with paclitaxel-induced toxicity phenotypes may be functioning through gene regulatory mechanisms. Interestingly, neither GWAS alone was enriched for eQTLs. Thus, our integration method may be reducing noise and revealing important functional SNPs. An enrichment of eQTLs has previously been shown in SNPs associated with 6 other chemotherapeutic drugs, which indicates that susceptibility to these drugs may depend on subtle gene expression differences across individuals (31).

The enrichment analyses were likely affected by the different linkage disequilibrium patterns among the populations studied. The LCL GWAS was a meta-analysis of African, African American, and European populations, whereas the patient GWAS was conducted in Europeans. In the meta-analysis, SNPs that are associated with cytotoxicity in all populations are prioritized over those associated in only one of the populations. We may have missed identifying European-specific overlap alleles. However, because the population linkage disequilibrium patterns differ and because African populations have shorter linkage disequilibrium blocks, overlap SNPs are more likely to be functional SNPs rather than SNPs that simply tag a functional locus (37).

We functionally assessed the involvement of one overlap eQTL, RFX2, in the NS-1 neuropathy cell model. Paclitaxel has previously been shown to decrease neurite outgrowth in the parent clone of the NS-1 cell line (36). Here, we showed that decreased expression of Rfx2 sensitizes NS-1 cells to paclitaxel by reducing neurite outgrowth and survival. This result validates our approach by showing that patient neuropathy and LCL cytotoxicity overlap analyses can reveal genes mechanistically involved in paclitaxel response. Although most previous work on RFX2 in mammalian cells describes its role in spermatogenesis (38, 39), several studies point to a potential role for the protein in sensory neurons. RFX2 and the related protein RFX1 have been shown to directly bind and regulate the transcription of ALMS1 (40). Mutations in ALMS1 cause the rare genetic disorder Alström syndrome, which is characterized by neurosensory degeneration, metabolic defects, and cardiomyopathy (40). In addition, the regulatory factor X transcription factors present in Caenorhabditis elegans and Drosophila, which are called DAF-19 and RFX, respectively, regulate ciliated sensory neuron differentiation (41, 42).

Upon knockdown of Rfx2 in NS-1 cells, the potential target gene Cyp51 also decreased expression, which was the expected direction of effect based on the preliminary gene expression analysis in LCLs. However, CYP51A1 does not contain an X-box RFX-binding domain (43) in the promoter region (2 kb upstream of the transcription start site), which means it is unlikely a direct target of RFX2 and may instead be further downstream in the pathway. Alternatively, RFX2 could be regulating an enhancer of CYP51A1 that is further outside the gene region. CYP51A1 is a member of the cytochrome P450 superfamily of enzymes, which catalyze many reactions involved in the metabolism of drugs and endogenous compounds. Specifically, CYP51A1 is known to participate in the synthesis of cholesterol (44). CYP51A1 has not been previously implicated in paclitaxel metabolism (45).

In the CALGB GWAS, one of the top SNPs that associated with patient paclitaxel-induced sensory peripheral neuropathy (rs10771973, P = 2.6 × 10⁻⁶) is located in FGD4 (3). Mutations in FGD4 can cause the congenital peripheral neuropathy Charcot–Marie–Tooth disease type 4H, and thus the gene is a plausible candidate for involvement in variation in peripheral neuropathy induced by paclitaxel. This SNP association was replicated in a second cohort of self-reported White patients with breast cancer (n = 154; P = 0.013) and in a cohort of self-reported African American patients with breast cancer (n = 117; P = 6.7 × 10⁻³; ref. 3). However, this SNP was not associated with paclitaxel-induced cytotoxicity in LCLs (P = 0.65). FGD4 is not expressed in LCLs (21), and thus the SNP is not expected to function in this model system. While our integrative approach can reveal variants and genes acting in paclitaxel response in both patients and LCLs, it does not identify genes potentially acting in patients that are not expressed in LCLs.

Effectively, LCLs have been used as an additional cohort to study the pharmacogenomics of various chemotherapeutics (16–19) because limited resources and in vivo confounders make obtaining large, homogeneous patient cohorts difficult. Here, we saw greater SNP overlap than expected by chance between SNPs associated with paclitaxel-induced cytotoxicity in LCLs and SNPs associated with paclitaxel-induced sensory peripheral neuropathy in patients at multiple P value thresholds, which confirms a role for the LCL model in the analysis of at least a subset of genes involved in patient neurotoxicity. This significant enrichment among a relatively large number of top SNPs is consistent with an underlying polygenic architecture for paclitaxel-induced toxicity. Functional siRNA studies in the NS-1 neuropathy model validated the involvement of RFX2 in paclitaxel toxicity, supporting our multi-gene hypothesis. Our novel integrative enrichment approach that combines clinical and LCL GWAS results can be used to expand patient cohort sizes for any drug phenotype of interest, including other toxicities, such as neutropenia, to find genes of potential impact that can be studied in cellular models.

E.P. Winer has a Commercial Research Grant from Roche. No potential conflicts of interest were disclosed by the other authors. The content of this article is solely the responsibility of the authors and does not necessarily represent the official views of the National Cancer Institute (NCI).

Conception and design: H.E. Wheeler, E.R. Gamazon, K. Owzar, M. Kubo, C. Hudis, L.N. Shulman, Y. Nakamura, M.J. Ratain, N.J. Cox, M.E. Dolan

Development of methodology: H.E. Wheeler, E.R. Gamazon, C. Wing, K. Owzar, N.J. Cox

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): H.E. Wheeler, C. Wing, U.O. Njiaju, C. Njoku, R.M. Baldwin, M. Kubo, H. Zembutsu, C. Hudis, L.N. Shulman, Y. Nakamura, D.L. Kroetz, M.E. Dolan

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): H.E. Wheeler, E.R. Gamazon, C. Wing, R.M. Baldwin, K. Owzar, D. Watson, I. Shterev, C. Hudis, L.N. Shulman, M.J. Ratain, D.L. Kroetz, N.J. Cox, M.E. Dolan

Writing, review, and/or revision of the manuscript: H.E. Wheeler, E.R. Gamazon, C. Wing, R.M. Baldwin, K. Owzar, C. Jiang, H. Zembutsu, E. Winer, C. Hudis, L.N. Shulman, M.J. Ratain, D.L. Kroetz, N.J. Cox, M.E. Dolan

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): C. Wing, C. Hudis

Study supervision: H.E. Wheeler, C. Hudis, L.N. Shulman, M.E. Dolan

The authors thank the Pharmacogenomics of Anticancer Agents Research LCL core for maintaining and distributing LCLs.

This study is supported by NIH/NCI Cancer Biology Training grant T32CA009594 and National Research Service Award F32CA165823 (to H.E. Wheeler), NIH/National Institute of General Medical Sciences (NIGMS) Pharmacogenomics of Anticancer Agents grant U01GM61393 (to M.J. Ratain, N.J. Cox, M.E. Dolan, and K. Owzar), NIH/NIGMS Pharmacogenetics of Membrane Transporters grant U01GM61390 (to D.L. Kroetz), University of Chicago SPORE Grant NCI P50CA125183 (to M.E. Dolan), NIH/NCI Medical Oncology Training grant T32CA009566 (to U.O. Njiaju), and the Biobank Japan Project funded by the Japanese Ministry of Education, Culture, Sports, Science and Technology. This work is part of the NIH Pharamcogenomics Research Network-RIKEN Center for Genomic Medicine Global Alliance.

The research for CALGB 60202 and 40101 was supported, in part, by grants from the NCI (CA31946) to the CALGB (Monica M. Bertagnolli) and to the CALGB Statistical Center (Daniel J. Sargent, CA33601).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.