A Genome-Wide Association Study Identifies a Novel Locus for Bortezomib-Induced Peripheral Neuropathy in European Patients with Multiple Myeloma

Some patients with multiple myeloma have subclinical or even clinical peripheral neuropathy at diagnosis. This peripheral neuropathy can be related to comorbidities, such as diabetes mellitus, or associated with the M-protein itself. In the course of the disease, peripheral neuropathy is mostly induced by therapies, especially thalidomide (thalidomide-induced peripheral neuropathy, TiPN) and bortezomib (bortezomib-induced peripheral neuropathy, BiPN), which may be considered as distinct clinical entities (1). TiPN may arise after prolonged administration of thalidomide (in 30%–55% of patients treated for 12 months, including 15%–25% with grade 2 or higher peripheral neuropathy) and appears to be due to a cumulative effect. Initial symptoms include sensory changes, such as paresthesia and hyperesthesia, later followed by motor symptoms and autonomic dysfunction. BiPN is characterized by neuropathic pain and a length-dependent distal sensory neuropathy with suppression of reflexes. Motor neuropathy may follow and infrequently results in mild to severe distal weakness in the lower limbs. There may also be a significant autonomic component, which manifests as dizziness, hypotension, diarrhea or constipation, and/or extreme fatigue. BiPN is thought to occur at a certain threshold of treatment (within five cycles but rarely beyond) in 40% to 60% of the patients, including 15% to 40% who will develop severe peripheral neuropathy (grade 2 or higher). This drug-induced toxicity is well known by physicians and nurses, and patients are now systematically informed about these potential side effects. The use of subcutaneous bortezomib reduces the incidence of BiPN but does not abrogate this toxicity (2). As no effective prophylactic treatment is available, prompt action in case of symptoms, including dose reduction and weekly administration of bortezomib, is crucial to manage this severe toxicity, which may dramatically affect the quality of life (3–5). Therefore, the identification of patients at risk of developing BiPN or TiPN is an important issue. This is especially true because the triplet combination of bortezomib–thalidomide–dexamethasone is considered one of the best induction regimens prior to high-dose therapy and autologous stem cell transplantation for the treatment of younger patients with de novo multiple myeloma (6).

The interindividual differences in the onset of BiPN or TiPN is in agreement with an underlying genetic susceptibility to this toxicity. Rare variants in bortezomib or thalidomide target proteins could affect the patient's sensitivity to these drugs. Among the pharmacogenomic methods to discover genetic loci associated with drug-induced toxicities, the candidate gene approach has shown a significant genetic contribution to the risk of developing TiPN or BiPN (7–10). However, a genome-wide association study (GWAS) has the capacity to identify new genetic variants that will have a direct or indirect effect on drug sensitivity. Here we report the results of a GWAS of 583 patients with multiple myeloma treated with bortezomib to discover genetic variants associated with severe BiPN. This is the first GWA pharmacogenomic study of bortezomib treatment toxicity and provides novel insights into bortezomib-related pathways.

Peripheral blood DNA samples were collected from 598 patients with newly diagnosed multiple myeloma who received bortezomib–dexamethasone (VD) induction therapy. Patients were treated in randomized clinical trials of the Intergroupe Francophone du Myelome (IFM; IFM 2005-01, IFM2007-02) or routine practice in France (n = 482) and in a randomized clinical trial of the Dutch/Belgian Haemato-Oncology Foundation for Adults in the Netherlands (HOVON) and the German-Speaking Myeloma Multicenter Group (GMMG; HOVON-65/GMMG-HD4; n = 116). The IFM VD treatment consisted of four 3-week cycles of bortezomib 1.3 mg/m² administered intravenously on days 1, 4, 8, and 11 plus dexamethasone 40 mg days 1–4 (all cycles) and days 9–12 (cycles 1 and 2). The HOVON 65/GMMG-HD4 VD treatment consisted of three cycles of bortezomib 1.3 mg/m² administered intravenously on days 1, 4, 8, and 11 plus dexamethasone 40 mg on days 1–4, 9–12, and 17–20 (patients enrolled in the HOVON-65/GMMG-HD4 trial received doxorubicin 9 mg/m²/day on days 1 through 4, in addition to VD according to the bortezomib, doxorubicin, and dexamethasone (PAD) regimen. Adverse events including peripheral neuropathy were graded by NCI Common Toxicity Criteria Version 3.0. All patients provided written informed consent for both the treatment and companion protocols.

Data quality assessment and control steps carried out during GWAS are summarized in Supplementary Fig. S1. A total of 482 multiple myeloma samples in the exploratory IFM cohort and 116 multiple myeloma samples in the Dutch and German replication cohort were genotyped using Affymetrix SNP6.0 Human DNA chips. Affymetrix CEL files were analyzed either by using Affymetrix Genotyping Console software v4.0 (GTC 4.0), followed by application of the Affymetrix Birdseed algorithm v2.0 to generate SNP genotype calls for the IFM exploratory cohort (GEO accession GSE65777) or by application CRLMM v2 algorithm to generate SNP genotype calls for the replication cohort (GEO accession GSE66903).

Stringent quality control (QC) thresholds were applied to filter out poorly genotyped subjects: if contrast QC <0.4, call rate <97% and outlying heterozygosity rate (het_rate > mean het_rate +3SD), the individual was removed. Principal component analysis (PCA) was performed to visualize the genetic ancestry of the IFM samples that passed the QC and assess whether population adjustment should be made (Supplementary Fig. S2). Random 60,000 genotypes of IFM subjects (IFM, n = 469) and unrelated individuals from three HapMap phase III populations representing Northwest European (CEU, n = 162), African (YRI, n = 163), and Chinese (CHB, n = 82) ancestries were combined to calculate the PCA. This method identified samples not clustered with Northwest European group (IFM outliers, n = 34), given that these patients were equally distributed between the cases and controls groups (Fisher exact test P = 0.36), no adjustment was needed, and therefore they were kept for the GWAS. Inspection of the observed and expected distribution of the neuropathy association statistic showed absence of hidden population substructure (Cochran–Armitage test of association; genomic inflation factor λ = 1.05).

SNP QC was conducted in four steps to remove suboptimal markers of the GWA data (Supplementary Fig. S1). (i) unannotated SNPs according to hg19 na32 SNP6 Affymetrix annotations (n = 130) along with SNPs from mitochondrial and sex chromosomes (n = 37,326) were not considered in the study, (ii) SNPs with missing genotype in more than 5% of the subjects (n = 16,743), (iii) SNPs of low minor allele frequency (MAF) less than 5% (n = 483,984), (iv) SNPs showing extensive deviation from Hardy–Weinberg equilibrium (HWE) with an HWE P < 1 × 10⁻⁵ (n = 834).

Statistical analyses were performed using SNPTEST v2.5 (11).First, we compared 370,605 genotypes from 155 grade ≥ 2 BiPN IFM patients to 314 control IFM patients defined as grade 1 BiPN or no BiPN patients. Second, we performed a validation using the HOVON-65/GMMG-HD4 cohort for the highest associated SNPs (P_trend< 10⁻⁵) as identified in the exploratory cohort. We compared 41 bortezomib-treated grade ≥ 2 BiPN patients with 75 bortezomib-treated control patients. We applied a one-sided logistic regression with 10,000 label-swapping permutations to correct for multiple testing to confirm BiPN association in this independent cohort. The predictive value of the SNP validated in the external series was assessed on the overall population (n = 583, i.e., 195 cases and 388 controls) with 10⁷ label-swapping permutations.

We conducted a pharmacogenomic GWA study to identify genetic variants associated with bortezomib toxicity in newly diagnosed patients with multiple myeloma who received VD induction therapy. Using SNP 6.0 Affymetrix arrays, we genotyped 909,622 tagging SNPs in 482 multiple myeloma cases. Of the 482 DNA samples genotyped, 469 cases passed strict QC criteria (Supplementary Fig. S1). We considered only the 370,605 autosomal SNPs with homozygosity in at least 5% of patients, a genotype call in at least 95% of patients and with an HWE P > 10⁻⁵. We compared the genetic contribution of patients who developed BiPN of grade ≥ 2 (n = 155) with that of patients who did not develop severe BiPN or without BiPN (n = 314). We separated grade 0, 1 versus 2 or more based on the clinical impact of such a toxicity. Grade 1 neuropathy requires a careful follow-up, but doses of bortezomib are not modified. Doses of bortezomib in the routine clinical practice must be adapted (from 1.3 mg/m² to 1 mg/m², or from the biweekly to the weekly schedule administration) according to the onset of grade 2 peripheral neuropathy, or stopped in case of grade 3 or more, and resumed in case of recovery. The GWA study showed association for six SNPs with OR > 1.8 and P_trend< 10⁻⁵ (Table 1 and Supplementary Table S1; Supplementary Fig. S3 and S4), although none reached the actual significance in a GWA study (0.05/370,065 = 1.35 × 10⁻⁷). To replicate these findings, a validation was performed using SNP 6.0 Affymetrix arrays in 114 newly diagnosed patients with multiple myeloma enrolled in the HOVON-65/GMMG-HD4 clinical trial who received VD induction therapy. A significant association was seen for rs2839629 (OR, 2.04; 95% CI, 1.11–3.33; P = 8.3 × 10⁻³; Table 2) which maps within the 3′ UTR of PKNOX1 (transcription factor PBX/knotted 1 homeobox 1). The overall estimate of BiPN associated with rs2839629 was an OR of 1.89 (95% CI, 1.45–2.44; P = 5 × 10⁻⁷). Moreover, rs2839629 is in strong linkage disequilibrium (LD) with rs76516641 (r² = 0.94) and rs915854 (r² = 0.86) which map within the intergenic region of 19.5 kb between PKNOX1 and cystathionine-ß-synthetase (CBS; Fig. 1).

Both PKNOX1 and CBS appear to be strong candidates for BiPN susceptibility genes. PKNOX1 is known to modulate transcriptional activity of chemokine monocyte chemoattractant protein-1 (MCP-1) gene (12–15). Through interaction with its cognate receptor CCR2, MCP-1 contributes to paclitaxel CIPN through changes in dorsal root ganglion neurons (16). MCP-1 is universally increased in different models of neuropathic pain and may be considered as a biomarker of chronic pain (17). MCP-1 is an important mediator of macrophage-related neural damage in different animal models of inherited neuropathies and acute inflammatory demyelinating neuropathy (18, 19). CBS encodes the endogenous H₂S-producing enzyme CBS. CBS–H₂S signaling pathway is implicated in the pathogenesis of a variety of neurodegenerative and inflammatory disorders, diabetic gastric hypersensitivity and plays a crucial role in inflammatory pain in temporomandibular joint (20–23).

To explore the possibility that this association might be mediated through differential expression of PKNOX1 or CBS or both, we examined the correlations between rs2839629, rs76516641, and rs915854 genotypes and tissue-specific gene expression levels by using the expression quantitative trait locus analysis available on the SNiPA portal (http://www.snipa.org; ref. 24) that used GTEX Portal v6 (25) and MuTHER consortium (26) as primary sources. PKNOX1 expression was significantly associated with rs2839629 and rs915854 genotypes in tibial nerve tissue (P = 5.6 × 10⁻⁸ and P = 1.9 × 10⁻⁷, respectively; Supplementary Table S2 and S3) with higher expression associated with rs2839629 risk alleles (Fig. 2), whereas CBS expression was significantly associated only with rs2839629 in skin (P = 2.6 × 10⁻¹⁵) and in blood (P = 3.1 × 10⁻⁸; Supplementary Table S2). rs915854 is annotated with a regulatory feature cluster characterized by histone marks H3K27ac and H3K4me1 enrichment in blood and cervix cells (Supplementary Table S3). In addition, rs915854 minor allele is predicted to disrupt binding site for the general pioneer factor FOXA1, this could render the enhancer less active for target gene expression such as PKNOX1 (Fig. 3). Conversely, rs76516641 genotype has no significant effect on both PKNOX1 and CBS gene expression (data not shown).

As we show that the risk allele A for rs2839629 is associated with higher levels of PKNOX1 expression and previous report have demonstrated that PKNOX1 binds preferentially to the -2578G (rs1024611G) polymorphism leading to increase MCP-1 levels (14), we analyzed the relationship between rs2839629A and rs1024611G in the IFM exploratory cohort. We found a significant association between rs2839629 A/A homozygous genotype and rs1024611G-bearing allele (Fisher exact test P = 0.01) suggesting a possible epistatic interaction between rs2839629 and rs1024611 to regulate MCP-1 expression.

The current pharmacogenomic GWA study also confirmed the modest association of the rs619824 genotype with BiPN (P =0.043) previously identified by Corthals and colleagues (8). Although there has been no overlap with a previous study on late-onset of BiPN-associated variants reported by Broyl and colleagues (10) as shown in Supplementary Table S4. This lack of overlap could reflect the potential complexity of predisposition to BiPN. More importantly, the design of the custom SNP chip used previously only contained 3404 SNPs in 983 hypothesis-driven genes which were thought to be functionally relevant in abnormal cellular functions, inflammation and immunity, as well as drug responses rather than adverse drug reactions which are less obvious candidates (27).

To date there are no established predictors of BiPN, it is impossible to predict which patient will develop neuropathy. Previous studies performed by our group and others using candidate gene approach have revealed significant association between SNP and BiPN; however, the clinical relevance of these findings is not clear (8–10). To increase our chance to discover variants that might provide new insights in the mechanisms underlying gene phenotype, we used a hypothesis-free approach. GWAS in cancer pharmacogenomics is challenging and few reports have been published to date. This is mainly due to insufficient statistical power in studies (28). To partially overcome these limitations, we designed our analysis to identify high-effect SNP (OR > 1.8) with MAFs greater than 0.05, in a large cohort of IFM patients with multiple myeloma (n = 469) uniformly treated to achieve convincing statistical power (29) and we verified our findings in an independent cohort of Dutch and German patients with multiple myeloma. Furthermore, our GWAS approach eliminates selection case–control bias as both case–control studies included patients in cohort studies, that is, IFM or HOVONGMMG clinical trials cohorts. When evaluating toxicity, it is sometimes difficult to distinguish between BiPN and neuropathic pain in general. It is also recognized that the sole use of the NCI CTC for assessment of sensory peripheral neuropathy is suboptimal. It is also recognized that detailed patient-reported symptom data and a quality-of-life assessment more accurately describes this toxicity and that physician-reported NCI-CTC grading underreports peripheral neuropathy. These systematic evaluations are difficult to apply in a multicenter study in the context of pharmacogenomics analyses. Of note, our study has enrolled patients without peripheral neuropathy at baseline, and patients were treated with the doublet combination of bortezomib and dexamethasone, and did not receive other neurotoxic agents.

Our analysis revealed an SNP associated with BiPN (rs2839629; OR, 1.89; P = 7.6 × 10⁻⁶) that was replicated in an independent cohort (OR, 2.04; P = 8.3 × 10⁻³) in high LD with SNP rs915854. Both variants are in noncoding regions; rs2839629 is located in the 3′UTR of PKNOX1, and rs915854 is in the intergenic region between PKNOX1 and CBS. Expression quantitative trait loci showed that these variants alter PKNOX1 and CBS expression presumably via cis-regulatory elements in the case of rs915854 as it falls within a regulatory region (Fig. 3). Finally, we found a significant association between rs2839629A and rs1024611G that could have an impact on MCP-1 expression levels. Given that these genes encode proteins, directly or indirectly, involved in neuropathic (16) and inflammatory pain (23), the functional significance of these predictive SNPs is established. This discovery opens the way to investigate novel pathways linked to PKNOX1 and CBS activities for a better understanding of mechanisms underlying this neurotoxicity.

This work generated a new hypothesis hypothesis regarding neurotoxicity mechanisms and provides new targets for neuroprotective strategies; however, additional international collaborative efforts including non-European countries are warranted to confirm or refute these findings and examine the impact of differential expression of both PKNOX1 and CBS effects on bortezomib exposure in cell model. Our results are preliminary and cannot be proposed yet for a systematic use in the routine clinical practice. Nevertheless, our findings are one of the first steps that may allow for the identification of patients at increased risk of severe BiPN, and these patients may benefit from the use of alternative drugs, such as carfilzomib, and/or a more focused clinical management of this toxicity.

R. Kuiper is an employee of SkylineDx. N.C. Munshi is a consultant/advisory board member for Takeda. P. Sonneveld reports receiving commercial research grant from Celgene and Janssen; and is a consultant/advisory board member for Celgene and Janssen. P. Moreau reports receiving speakers bureau honoraria from Amgen, Bristol-Myers Squibb, Celgene, Janssen, Novartis, and Takeda. No potential conflicts of interest were disclosed by the other authors.

Conception and design: F. Magrangeas, R. Kuiper, H. Avet-Loiseau, P. Sonneveld, C. Dumontet, P. Moreau, M. van Duin, S. Minvielle

Development of methodology: R. Kuiper, C. Dumontet, L. Campion

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): P. Sonneveld, C. Dumontet, P. Moreau, M. van Duin

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): F. Magrangeas, R. Kuiper, H. Avet-Loiseau, W. Gouraud, C. Guérin-Charbonnel, A. Aussem, H. Elghazel, N.C. Munshi, P. Sonneveld, M. van Duin, L. Campion, S. Minvielle

Writing, review, and/or revision of the manuscript: F. Magrangeas, R. Kuiper, H. Avet-Loiseau, C. Guérin-Charbonnel, M. Attal, N.C. Munshi, C. Dumontet, P. Moreau, M. van Duin, L. Campion, S. Minvielle

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): L. Ferrer, J. Suhard, H. Der Sakissian, M. van Duin

Study supervision: F. Magrangeas, C. Dumontet, S. Minvielle

Other (technical and computational analysis support): H. Der Sakissian

The authors thank Elise Douillard, Magali Devic, Emilie Maurenton, and Nathalie Roi for excellent technical expertise.

This study was supported by Intergroupe Francophone du Myélome; the French Institute National du Cancer Grant POBNIMM PROG/08/10 (to H. Avet-Loiseau, S. Minvielle); and NIH grant PO1CA155258-01 (to S. Minvielle, H. Avet-Loiseau, N.C. Munshi).

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.