Targeted Sequencing Reveals Low-Frequency Variants in EPHA Genes as Markers of Paclitaxel-Induced Peripheral Neuropathy

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

The anticancer agent paclitaxel is a microtubule inhibitor widely used in the treatment of many solid tumors (1). Peripheral neuropathy is its dose-limiting toxicity (2), and severe neuropathy cases with an important reduction in the quality of life of the patients are not rare (3, 4). The lack of effective treatments for the neuropathy creates an urgent need to identify markers that can help to personalize treatment and avoid severe neuropathy events. The patient genetic background has been proposed to play a relevant role in the susceptibility for suffering neuropathy (5). In this regard, genes in paclitaxel pharmacokinetic (6, 7) and pharmacodynamic (8, 9) pathways have been included in studies of candidate genes and, more recently genome-wide association studies (GWAS) have been performed (10, 11).

Candidate gene studies, by us and other groups, have demonstrated that common variants in paclitaxel metabolizing enzymes and paclitaxel target [i.e., CYP2C8*3 (12–15), CYP3A4*22 (7), TUBB2A rs909964, and rs909965 (8, 9)] influence neuropathy risk, whereas genome-wide genotyping has uncovered novel genes (10, 11). A GWAS by our group (11) suggested that the EPHA gene family, which plays a key role in the development of nervous system and in nerve injury repair (16–18), was a key player for paclitaxel neuropathy susceptibility. Meta-analysis of GWAS top hits showed that EPHA5 rs7349683 reached genome-wide significance (11), and follow-up studies further supported that this variant (19), EPHA6 rs301927 (9, 19), and EPHA8 rs209709 (19) moderately increased paclitaxel-induced neuropathy risk. However, large part of the variation in paclitaxel-induced neuropathy remains unexplained.

Low-frequency variants with strong effects may contribute to the neuropathy variability observed in patients. To investigate this hypothesis, sequencing technologies are required and, so far, only two exploratory studies following different strategies have been performed. In one, we applied whole-exome sequencing to few extreme neuropathy patients, and identified defective CYP3A4 variants associated with the neuropathy (20). The second study sequenced genes causative of familial polyneuropathies (Charcot–Marie–Tooth, CMT), and suggested ARHGEF10 and PRX as chemotherapy-induced neuropathy markers (21). These initial studies are promising; however, the statistical power for a whole-exome sequencing study is low, and in the CMT analysis, key genes were not included.

Here, we performed targeted exome sequencing of genes with common variants associated with paclitaxel-induced neuropathy (EPHA4, EPHA5, EPHA6, and EPHA8) plus genes involved in paclitaxel pharmacokinetics and in CMT. In total, we sequenced 39 genes in 228 selected patients with high or no/low paclitaxel-induced neuropathy. The strongest association corresponded to EPHA6, and the relevance of low-frequency EPHA5/6/8 nonsynonymous coding variants was validated in an independent cohort of 202 paclitaxel-treated patients. These results reveal EPHA genes as key players in chemotherapy-induced neuropathy and stress the importance of gene sequencing for identifying genetic risk factors of neuropathy.

The discovery cohort was derived from a series of 449 breast or ovarian cancer patients treated with paclitaxel (97% in first line), with DNA available, with clinical data and with a homogenous neuropathy assessment (20). Part of the patients has been described in previous studies (15, 19, 20). From these, 228 were selected for whole or targeted exon deep-sequencing, based on extreme-neuropathy phenotype. Among them, 131 were high-neuropathy patients that fulfilled the following criteria: grade 3 or 2 neuropathy (NCI-CTC v4) during paclitaxel treatment, no neuropathy risk factors (diabetes, alcoholism, AIDS, or previous neuropathies), and treatment modifications due to neuropathy (dose reduction or treatment suspension) or neuropathy that lasted >6 months after paclitaxel treatment finished. The remaining 97 patients were no/low-neuropathy patients with no neuropathy signs or grade 1 neuropathy after receiving paclitaxel (Table 1).

The validation of results was performed in an independent series of 202 paclitaxel-treated patients with neuropathy data recorded cycle by cycle. Most patients had breast or ovarian tumors, 109 were Spanish (54%) and 93 Swedish (46%). One hundred and twenty-nine samples corresponded to a previous GWAS study (11), 37 to Spanish patients already described (19), and 36 samples were new cases collected in Spain. From all patients, cumulative paclitaxel dose up to grade 2 (NCI_CTC v2/4) neuropathy was available (Table 1).

All individuals participating in the study were over 18 years of age, had been diagnosed of cancer with histologic confirmation, a life expectancy of ≥12 weeks and Eastern Cooperative Oncology Group performance status ≤2, adequate bone marrow, and renal and hepatic function. The recruitment of patients and collection of samples was approved by local internal ethical review committees, and all patients gave written informed consent to participate in the study.

From the 228 patients used in the discovery series, 196 samples were processed using the TruSeq Custom Amplicon Kit (Illumina) covering the coding plus 25 bp intronic flanking region of 39 genes that included: EPHA4, EPHA5, EPHA6, and EPHA8 (10, 11) plus additional genes involved in paclitaxel metabolism and transport (ABCB1, CYP2C8, CYP3A4, SLCO1B1, and SLCO1B3) and a selection of 30 genes associated with CMT hereditary peripheral neuropathies (Fig. 1). Very conserved CMT genes with no/very few variants reported were not selected for sequencing (e.g., ATL1, EGR2, GDAP1, GJB1, LMNA, PRPS1, RAB7A, and YARS). In brief, 150 ng of DNA extracted from peripheral blood (FlexiGene DNA Kit; Qiagen) was used to construct libraries and sequenced in a MiSeq sequencer (Illumina) with a paired-end mode using MiSeq Reagent Kit V3 (Illumina) and 600 cycles. In addition, whole-exome sequencing was performed on the remaining 32 patients [16 with high neuropathy (8 have been reported; ref. 20) and 16 patients with no neuropathy], as previously described (20). For the validation of the results, a TruSeq Custom Amplicon kit (Illumina) including the coding and intronic flanking region of EPHA5, EPHA6, and EPHA8 was used.

Targeted next-generation sequencing (NGS) data were demultiplexed with MiSeq Reporter (Illumina). Alignment was performed using Smith–Waterman algorithm (22) using GRCh37/hg19 assembly as reference, and Genome Analysis Toolkit v2 (GATK; ref. 23) was used for raw variant calling. For the 32 samples with whole-exome sequencing data, alignment and variant calling were performed by RUbioSeq software v3.7 (24). In this software, the alignment was performed using Burrows–Wheeler alignment (25), unmapped reads are realigned using BFAST (26), and for variant calling, GATK v2 was used (23). Variants were annotated with Snp Eff (http://snpeff.sourceforge.net/) and Variant Effect Predictor (http://www.ensembl.org/info/docs/tools/vep/index.html), and only nonsynonymous coding variants and those altering canonical splice sites, with P > 0.001 for Hardy–Weinberg equilibrium, were considered in subsequent steps. Supplementary Table S1 indicates gene and transcript references.

Variants included in the analysis were: (i) those previously described in public databases [dbSNP, http://www.ncbi.nlm.nih.gov/SNP/; Exome Aggregation Consortium (ExAC), http://exac.broadinstitute.org], and (ii) variants not previously described with: high variant call quality (Q>30), read depth >10X, and alternative variant frequency higher than 0.3 in at least one individual. Sequencing artifacts, defined as nucleotide changes detected in >20 samples in the sequencing panel but not present in ExAC, were omitted from the analysis. We defined loss of function (LOF) variants as those introducing stop codons (nonsense), variants disrupting canonical splice sites, and indels disrupting the reading frame. Template and configuration files for alignment and scripts are available at https://github.com/htejero/PaclitaxelNeuropathy.

Validation of variants was performed by Sanger sequencing with an ABI PRISM 3700 DNA Analyzer capillary sequencer (Applied Biosystems) on 3% of the LOF and missense variants included in the analysis.

Variants were classified as “common variants” if they had a minor allele frequency (MAF) ≥0.5% in the more than 30,000 sequenced non-Finnish Europeans from ExAC. Variants were classified as “low frequency variants” if they had an MAF <0.5% in the non-Finnish Europeans from ExAC and MAF <1% in 578 Spanish exomes from the CIBERER Spanish Variant Server (http://csvs.babelomics.org/). The purpose of including the Spanish data was to detect population-specific variants, and because of the small sample size (n < 600), the MAF threshold in this population was less stringent. For common variants, the frequency of each variant in the high versus no/low neuropathy group was compared with a χ² or Fisher test. For low-frequency variants, the association with paclitaxel-induced neuropathy was assessed with the gene-based Burden test (27) using the SKAT package and R statistical software (http://www.R-project.org/). Scripts are available at https://github.com/htejero/PaclitaxelNeuropathy. Based on statistical power calculations, only genes with ≥ four rare variants were included in the analysis.

The study followed a two-step design in which the best candidates from the discovery phase were selected for validation in an independent cohort of paclitaxel-treated patients (Table 1 shows discovery and validation series). No correction for multiple testing was performed. For samples with cycle-by-cycle neuropathy data, the association between EPHA variants and paclitaxel neuropathy risk was tested using the Kaplan–Meier analysis, modeling the cumulative dose of paclitaxel up to the development of neurotoxicity grade ≥2. Patients with no or low neuropathy (grade 0/1) were censored at total administered cumulative dose. We also evaluated the association using univariate and multivariable Cox regression analysis (14). Country of origin and treatment schedule (1 hour vs. 3 hour infusion) were included as covariates in the multivariate analyses. SPSS software package v.19 was used for these analyses. P values less than 0.05 were considered statistically significant.

NGS was performed on selected cases: 131 patients with high neuropathy (grades 2/3 that lasted a mean of 55 months) despite low accumulated paclitaxel dose (median, 1,295 mg) and 97 patients with no/low neuropathy (grades 0/1) despite high accumulated paclitaxel dose (median, 1,485 mg; Table 1). In addition, 33% of patients in the high neuropathy group had paclitaxel dose reductions or treatment suspensions caused by the neuropathy.

Sequencing of 39 candidate genes in the 228 patients identified 277 coding nonsynonymous or canonical splice site variants (266 missense, 3 in-frame deletions, 8 LOF; Supplementary Table S1). From these, 86 were common variants and 191 low-frequency variants.

At least one common variant was identified in each sequenced gene, except for CYP3A4, EPHA4, HSPB1, HSPB8, NEFL, NDRG1, and SPTLC2. When the presence of these common variants was compared among the neuropathy groups, association with paclitaxel neuropathy was found for only 2 SNPs located in CYP2C8 and PRX (P < 0.05; Supplementary Table S2).

The 191 low-frequency variants were distributed among all sequenced genes, except for NEFL and NGF. Of these 191 variants, 8 were LOF (3 altered canonical splice sites, 2 were nonsense variants, and 3 were indels causing frameshifts leading to premature stop codons; Table 2).

Analysis of the low-frequency variants identified EPHA6 as the gene most significantly associated with paclitaxel-induced neuropathy (Table 3). The five carriers of these variants were all high neuropathy patients with an amino acid change in the ephrin receptor ligand-binding domain of the protein. Remarkably, no EPHA6 variant carriers were present in the no/low-neuropathy group, suggesting a strong effect on neuropathy. One additional gene had this characteristic (SEPT9), but results did not reach statistical significance level. The other two EPHA genes analyzed, EPHA5 and EPHA8, have a similar biological function as EPHA6 (16–18), and also belonged to the high-neuropathy risk group of genes (Table 3). In EPHA5, five carriers had high neuropathy versus one with low neuropathy; and in EPHA8, nine carriers were in the high neuropathy and six in the no/low neuropathy group (Fig. 2; Supplementary Table S1). The highly conserved EPHA4, with only two variant carriers, one in each group, could not be analyzed.

Some of the discovery series patients had cycle-by-cycle neuropathy data available, and among these, three were carriers of low-frequency variants in EPHA5/6/8 genes (one variant in each gene). Accumulated paclitaxel dose analysis revealed that these patients had a significantly higher risk to suffer from neuropathy than patients without EPHA low-frequency variants [HR, 14.60; 95% confidence interval (CI), 2.33–91.62; P = 0.0042; Fig. 3A].

Sequencing EPHA5/6/8 in an independent cohort of 202 patients treated with paclitaxel and detailed cycle-by-cycle neuropathy data (Table 1) revealed 15 carriers of low-frequency missense variants in these genes (1 corresponded to EPHA6, 1 to EPHA5, and 13 to EPHA8). These variants were combined, and an accumulated paclitaxel dose analysis revealed that low-frequency EPHA5/6/8 variants conferred increased risk of neuropathy (HR, 2.07; 95% CI, 1.14–3.77; P = 0.017; Fig. 3B).

Combining discovery and validation series resulted in an HR of 2.50 (95% CI, 1.46–4.31) with a P value of 9.1 × 10⁻⁴ (Fig. 3C).

Paclitaxel-induced neuropathy is a clinically relevant toxicity affecting large number of cancer patients. Genetic variation has been shown to influence susceptibility to paclitaxel-induced neuropathy; however, a large part of the variation remains unexplained. Low-frequency variants with strong effects may explain part of the variability. To investigate this hypothesis, we performed massive sequencing of candidate genes in patients selected based on extreme-neuropathy phenotype. Gene-based analysis identified, for the first time, low-frequency genetic variants in EPHA5/6/8 as risk factors of chemotherapy-induced neuropathy. These results may provide a basis for personalizing paclitaxel treatment and decreasing the incidence of severe chemotherapy-induced neuropathies.

GWAS have identified common variants in EPHA genes with moderate effects on paclitaxel-induced neuropathy (EPHA5-rs7349683, EPHA6-rs301927, EPHA8-rs209709, and EPHA4-rs17348202; refs. 10, 11), and subsequent studies further supported the association of EPHA5, EPHA6, and EPHA8 polymorphisms (9, 19). Nonsynonymous coding variants, potentially affecting protein function, are expected to have stronger effects on neuropathy than common regulatory variants (28). Following this idea, we performed an NGS study in EPHA genes, together with paclitaxel pharmacokinetics and hereditary peripheral neuropathy-related genes. Gene-based analysis of our data revealed that low-frequency missense variants in EPHA6 increased paclitaxel-induced neuropathy risk. All these variants were located in the ephrin receptor ligand-binding domain, suggesting an alteration of the protein function and further supporting the association. EPHA5 and EPHA8 followed a similar trend (Fig. 2). In total, 15% (19 of 131) of patients in the high neuropathy group carried low-frequency nonsynonymous coding variants in EPHA5/6/8 genes. In the 202 patients of the validation series, 13 EPHA8 variant carriers were identified, but only 1 EPHA6 and 1 EPHA5 carriers were detected, suggesting that EPHA6 and EPHA5 variants (present in 5 of the 131 patients with high-neuropathy of the discovery) are less frequent in an unselected patient population, including many moderate-neuropathy patients (not represented in the discovery set). Thus, EPHA6 and EPHA5 variant carriers were scarce in the validation series, and the calculated EPHA effect mainly derived from EPHA8. Despite this, the accumulated dose analysis is a sensitive approach (15, 19) and was able to detect a statistically significant association. Altogether, these data suggest a relevant role for EPHA5/6/8 genes in paclitaxel-induced neuropathy and indicate a high impact of low-frequency variants missed in GWAS.

Eph receptors are tyrosine kinases involved in neural development (16) and nerve regeneration after damage (18, 29) among other functions: EphA4 controls axon sprouting/nerve regeneration after spinal cord injury (30–32); EphA5 plays an important role in the initiation of the early phases of synaptogenesis (33), and it has been found upregulated in mice with injured sciatic nerve (34); EphA6 is involved in neural circuits underlying aspects of learning and memory (35); and EphA8 induces neurite outgrowth through induction of sustained MAPK activity (36) and lack of this gene produces aberrant axonal projections (37). Knocking out EphA4, EphA5, EphA6, and EphA8 genes in mice results in viable and fertile animals with different neurological phenotypes. EphA4 knockout mice have gross motor dysfunction (38–40) and altered axonal regeneration and functional recovery following spinal cord injury (41). Knocking out the tyrosine kinase domain of EphA5 results in axon aberrations in topographic mapping and altered behavioral patterns (42, 43). EphA8 knockout mice have abnormal axonal projections in the spinal cord (37), and EphA6 knockout mice experienced behavioral deficits in learning and memory tests (35). Thus, these are crucial genes for neural development and nerve regeneration with a plausible link for the association found with paclitaxel-induced neuropathy.

In ExAC database, 0.1% of the European non-Finish population are carriers of LOF variants in either EPHA5, EPHA6, or EPHA8, and on >100,000 Islandic individuals, two complete human knockouts for EPHA5 and one for EPHA6 were identified (44). So far, no phenotype has been assigned to these individuals who are apparently healthy subjects. However, based on the literature and on our results, a high susceptibility to drug-induced neuropathy would be expected.

Concerning other genes potentially associated with the neuropathy, in line with Beutler and colleagues (20), we postulated that variants moderately affecting the function of CMT genes, while not being pathogenic, may increase the susceptibility to drug-induced neuropathy. We did not find low-frequency variants in PRX and common variants in ARHGEF10 associated with paclitaxel-induced neuropathy, although the second and third top protective genes were these two. For the ARHGEF10 common variant rs9657362, we also found a trend toward protection (21, 45). We also observed a trend toward increased neuropathy risk for other CMT genes (SEPT9 and SH3TC2). Variability in results among studies may be related to differences in neuropathy definitions/assessments, in tumor types and patient treatments, or in the distribution of low-frequency variants, which have shown to be population-specific. Thus, results need to be further explored and validated in large independent series.

With regard to the LOF variants detected in this study, three occurred in CMT genes (ARHGEF10, IKBKAP, and DHTKD1). The patients with variants in ARHGEF10 and IKBKAP belonged to the no/low neuropathy group, in agreement with the fact that activating rather than LOF mutations in ARHGEF10 cause CMT (46) and that no phenotype is observed for IKBKAP heterozygous individuals (47). The variant in DHTKD1 was present in 2 patients with different neuropathy, but recent data question the role of this gene in CMT disease (48, 49). Two LOF variants affected EPHA genes (EPHA5 and EPHA8) and were found in high-neuropathy patients. One LOF variant occurred in the paclitaxel uptake transporter SLCO1B1, in a high neuropathy patient. Two occurred in CYP3A4, a gene in which we have demonstrated that defective variants increased neuropathy risk (20). Two patients were carriers of the CYP3A4*20 frameshift allele and belonged to the high-neuropathy group, but 1 patient with a splicing defect affecting the last exon belonged to the no/low neuropathy group. The effect of this latter variant on the splicing of the gene and how it affects function remains to be studied.

Although the main goal of this study was to identify neuropathy-associated low-frequency coding variants, we also found two common polymorphisms associated with the neuropathy: CYP2C8 rs1058930 (CYP2C8*4), for which previous studies have found contradictory results (9, 14), and PRX rs268674, which was associated with neuropathy risk here for the first time. Further studies should evaluate the relevance of these results.

Limitations of this study include gene selection, because relevant genes not yet connected with neuropathy susceptibility may have not been studied. There are also differences in the selection of patients in the discovery and validation series. In the discovery series, patients were mainly treated with paclitaxel as single agent, whereas in the validation cohort, the majority of the patients were treated with paclitaxel in combination with carboplatin. However, no major differences in neuropathy development between paclitaxel/carboplatin therapy versus paclitaxel as single agent exist (50, 51). In addition, we adjusted the analysis using treatment schedule as covariate. Nevertheless, using a more homogenous series may have resulted in stronger association results. Detection of low/moderate effects on neuropathy may require even larger samples sets, although the number of patients in this study is substantial and the neuropathy assessment was homogenously performed to reduce subjectivity (11, 20). On the whole, additional studies validating the results in extensive and well-characterized series of patients, the development of a model integrating all different risk markers identified, and providing with a standardized methodology to perform the genetic testing would be required to implement these risk factors into the clinics.

In conclusion, this study proves a relevant role of EPHA5, EPHA6, and EPHA8 genes in paclitaxel-induced neuropathy susceptibility and suggests that sequencing studies, rather than genotyping, would be adequate approaches to study genetic markers of neuropathy. Moreover, taking into account the role of these proteins in neural development and injury repair, EPHA variants may also confer increased neuropathy risk to many additional neurotoxic drugs. The final goal is to identify genetic risk factors that can help to personalize neurotoxic drug treatments and avoid severe chemotherapy-induced neuropathies that can seriously affect patients' quality of life.

No potential conflicts of interest were disclosed.

Conception and design: M. Apellániz-Ruiz, L. Inglada-Pérez, J. García-Donás, F. Al-Shahrour, A. Cascón, M. Robledo, C. Rodríguez-Antona

Development of methodology: M. Apellániz-Ruiz

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): L. Sánchez-Barroso, G. Gutiérrez-Gutiérrez, I. Calvo, A. Redondo, J. García-Donás, N. Romero-Laorden, M. Sereno, M. Merino, M. Currás-Freixes, V. Mancikova, E. Åvall-Lundqvist, H. Green

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M. Apellániz-Ruiz, H. Tejero, L. Inglada-Pérez, H. Green, F. Al-Shahrour, M. Robledo, C. Rodríguez-Antona

Writing, review, and/or revision of the manuscript: M. Apellániz-Ruiz, H. Tejero, L. Inglada-Pérez, L. Sánchez-Barroso, G. Gutiérrez-Gutiérrez, I. Calvo, B. Castelo, A. Redondo, J. García-Donás, M. Merino, M. Currás-Freixes, C. Montero-Conde, V. Mancikova, E. Åvall-Lundqvist, H. Green, F. Al-Shahrour, A. Cascón, M. Robledo, C. Rodríguez-Antona

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): L. Sánchez-Barroso

Study supervision: B. Castelo, J. García-Donás, C. Rodríguez-Antona

This work was supported by projects from the Spanish Ministry of Economy and Competiveness (grant number SAF2015-64850-R). M. Apellániz-Ruiz and Veronika Mancikova are predoctoral fellows of "la Caixa"/CNIO international PhD programme. Maria Currás is a predoctoral fellow supported by the Severo Ochoa Excellence Programme (project SEV-2011-0191). Cristina Montero-Conde is supported by a postdoctoral fellowship from the Fundación AECC. Part of the work was financially supported by grants from the Swedish Cancer Society, the Swedish Research Council, and LiU Cancer.

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.