Genetic factors influencing the risk of multiple myeloma bone disease

A major complication of multiple myeloma (MM) is the development of osteolytic lesions, fractures and bone pain. To identify genetic variants influencing the development of MM bone disease (MBD), we analyzed MM patients of European ancestry (totaling 3774), which had been radiologically surveyed for MBD. Each patient had been genotyped for ~6 00 000 single-nucleotide polymorphisms with genotypes for six million common variants imputed using 1000 Genomes Project and UK10K as reference. We identified a locus at 8q24.12 for MBD (rs4407910, OPG/TNFRSF11B, odds ratio=1.38, P=4.09 × 10–9) and a promising association at 19q13.43 (rs74676832, odds ratio=1.97, P=9.33 × 10–7). Our findings demonstrate that germline variation influences MBD and highlights the importance of RANK/RANKL/OPG pathway in MBD development. These findings will contribute to the development of future strategies for prevention of MBD in the early precancerous phases of MM.


INTRODUCTION
Multiple myeloma (MM) is a B-cell malignancy characterized by the expansion of clonal plasma cells in the bone marrow. 1,2 The disease is typified by varying numbers of osteolytic lesions that are the result of reduced osteoblastic and increased osteoclastic activity. 3,4 The presence of such bone lesions is a major criterion used to distinguish MM that requires treatment from precursor entities such as monoclonal gammopathy of undetermined significance and smoldering myeloma, which can be managed expectantly, 5 as such, it is important to understand factors impacting on the development of MM bone disease (MBD). Although most patients develop osteoblastic lesions, a subset is unaffected by MBD, the reasons for which are not understood. Understanding the genetic mechanisms that are responsible for these differences in MBD is a pressing clinical issue, which has important implications for the development of novel treatments and in determining which patients might benefit from alternate bone therapies. In this respect, it has been shown in the setting of osteoporosis that heritable factors account for 50-85% of the variation in bone mineral density (BMD). 6,7 As the interaction of bone-forming cells with MM cells is central to the development of an osteolytic lesion and MM clonal growth, we hypothesized that germline variation could also contribute to the development of MBD. We have recently conducted genome-wide association studies (GWAS) of MM searching for susceptibility alleles. [8][9][10] Linking these genetic data to the extent of MDB at baseline has allowed us to search for genetic variants influencing MBD risk.

METHODS Patients
We studied four independent cohorts of MM patients that had been the subject of previous GWAS [8][9][10][11]  Radiological assessment of bone lesions MBD was detected using axial survey in My9 and My11, axial skeletal survey (2001-2010) and whole body computed tomography (2011 onwards) in HdB and combined skeletal survey and skeletal computed tomography in ArK. The frequency of MBD was marginally higher in the HdB cohort (P = 0.01). Owing to the differences in sensitivity of the radiological methods used to detect MBD, patients were classified as either affected (MBD) or unaffected (no MBD). Age was not significantly associated with MBD in any of the four cohorts (that is, P40.05).

Tumor karyotyping
Conventional cytogenetics of tumors was conducted using standard karotyping methodologies, and standard criteria for the definition of a clone were applied. 13 Fluorescent in situ hybridization and ploidy classification of My9 and My11 samples were conducted using the methods described by Chiecchio et al. 14 Fluorescent in situ hybridization analysis and ploidy classification of HdB samples were performed as previously described. 15 The XL IGH Break Apart probe (MetaSystems, Altlussheim, Germany) was used to detect immunoglobulin H translocations in HdB samples.
Genotyping and quality control All cases were genotyped using Illumina Human OmniExpress arrays adhering to the manufacturer's protocols (Illumina, San Diego, CA, USA). Standard quality control was performed on all scans, excluding individuals with low call rate (o 90%) and extremely high or low heterozygosity

Imputation
Genotypes for common variants across the genome were imputed using data from 1000 Genomes Project (phase 1 integrated release 3, March 2012) and UK10K as reference in conjunction with IMPUTE2 v2.1.1 16 after pre-phasing with SHAPEIT software; 17 poorly imputed SNPs defined by an information measure o0.90 were excluded. All genomic locations are given in National Center for Biotechnology Information Build 37/UCSC hg19 coordinates. All SNPs having a minor allele frequencyo1% were excluded.

Statistical analysis
We compared the relationship between genotype and presence or absence of MBD by logistic regression including covariates found in univariate analysis to influence MBD. We adjusted for the method of radiological assessment in each cohort and used eigenvalues in the analysis of the HdB series, to adjust for population substructure. P-values presented correspond to the significance of a test difference among all three of the genotype groups (common allele homozygote, heterozygote and rare allele homozygote). We confined our analysis to SNPs with a minor allele frequency 41% because of extreme value of the test statistics. Overall statistical significance tests for each SNP were performed by combining the results for each cohort by a fixed-effects meta-analysis. All statistical tests were two-sided. Inflation of the test statistics, λ, was estimated by dividing the 45th percentile of the test statistic by 0.357-the 45th percentile for a χ 2 distribution on 1 degree of freedom. Betweenstudy heterogeneity was quantified using the I 2 statistic. Associations were regarded as statistically significant at a P-value ⩽ 5.0 × 10 -8 (that is, genome-wide significance). All statistical analyses were performed using PLINK v1.07 18 and R (v3.1.3) software. 19 Functional prediction To explore the epigenetic profile of genomic location associated with MBD, we used ENCODE histone modification data and HaploReg and RegulomeDB 20,21 to examine whether any of the SNPs or their proxies (that is, r 2 40.8 in the 1000 Genomes European reference panel) annotate transcription factor binding or enhancer elements. We assessed sequence conservation using Genomic Evolutionary Rate Profiling scores. 22

Relationship between genotype and bone lesions
After quality control measures were applied, genotype data on 55 31 610 SNPs was available for 3774 MM cases with MBD data. Quantile-Quantile plots of test statistics of the relationship between SNP genotype and MBD for each of the four cohorts is shown in Supplementary Figure 3; inflation factors λ = 1.002-1.01, λ = 1.0008 for the meta-analysis (Supplementary Figure 4).
Nine SNPs showed an association with MBD and reached genome-wide significance (Figures 1 and 2, Supplementary Table 1). All nine SNPs were located in the same region at 8q24.12 and were in strong linkage disequilibrium (LD). The strongest association at 8q24.12 was provided by the common SNP rs4407910 (risk allele frequency = 0.50, odds ratio = 1.38, 95% confidence interval = 1.24-1.54, P = 4.02 × 10 − 9 ). The association was consistent across each of the four patient cohorts (Figure 3 and there was no significant between-study heterogeneity (P het = 0.44, I 2 = 0%). An increased prevalence of MBD has been observed in male patients and hyperdiploid MM. 23,24 Stratifying data by sex or ploidy did not provide evidence for a differential effect of rs4407910 genotypes on the risk of MBD ( Supplementary  Figures 5 and 6). Variation at 8q24.12 has been previously associated with BMD and osteoporosis, but not as strongly with risk of fracture. [25][26][27][28] Meta-analysis of My9 and My11 data showed no evidence that rs4407910 genotype influenced the risk of vertebral body fracture (P = 0.28).
Impact of alleles influencing BMD on MM bone disease Variation at 8q24.12 marked by rs23062375 and rs11995824, which are intronic SNPs in OPG and in LD with rs4407910 (respective LD metrics r 2 and D' -0.78, 0.93 and 0.79, 1.0) have been demonstrated to influence BMD. [25][26][27]31,33 To explore the possibility that other genetic variants influencing BMD also influence MBD, we investigated the association at 77 established risk loci for BMD with MBD [25][26][27][28] (Supplementary Table 4). Aside from the 8q24.12 SNPs, no other BMD locus showed an association with MBD after adjusting for multiple testing (that is, P40.001). Moreover, there was no over-representation of association for MBD across the 77 SNPs (Supplementary Table 5).

DISCUSSION
Our findings support the hypothesis that an individual's risk of developing MBD is influenced by germline variation. Specifically, we identified a locus at 8q24.12 (rs4407910) associated with MBD. rs4407910 maps to a region of LD that only contains OPG, and when taken in conjunction with eQTL data it is likely that reduced OPG expression is the functional basis of the 8q24.12 association. 31 Osteoclasts (OCs) are bone-resorptive cells, which are critical for the integrity of bone. OC-differentiation and activation is dependent on activation of nuclear factor-kB ligand (RANKL) signaling through the p38 MAPK pathway. [34][35][36] OPG is a negative regulator of bone resorption acting as a decoy receptor for RANKL, decreasing bone resorption through inhibiting differentiation of OC precursors, activating mature OCs and stimulating OC apoptosis. Germline inactivating mutations in OPG are responsible for the autosomal dominant diseases: early-onset and familial Paget'disease, familial expansile osteolysis and expansile skeletal hyperphosphatasia, which are characterised by the development of expansile osteolytic bone lesions. 37,38 Myeloma cells express RANKL and treatment of mice models of MM with OPG has been demonstrated to prevent lytic bone lesions, maintaining cancellous bone volume and inhibiting OC formation. 12,32 Current clinical management of MBD involves reducing myeloma cell infiltration of the bone marrow using  bisphosphonates to inhibit OC activity; however, there remains a need to develop more targeted treatments. Established MM therapies such as immunomodulatory drugs and proteasome inhibitors partially exert their activity through the OPG/RANK/ RANKL system. 39,40 Hence, targeting the OPG/RANK/RANKL system through specific agents such as raloxifene may have therapeutic potential. 41,42 Raising OPG levels directly by infusion of recombinant OPG (Fc-OPG) suppresses bone resorption. 43 Potential concerns over the generation of anti-Fc-OPG and binding to TRAIL has however shifted further development away from Fc-OPG as a RANK inhibitor, to Denosumab, a human monoclonal antibody against RANKL. 44,45 In addition to the 8q24.12 locus we identified a promising signal at 19q13.43 proximal to genes ZNF444 and GALP. GALP encodes a member of the galanin family of neuropeptides that are active within the central nervous system. 46 ZNF444, a zinc finger protein activates a scavenger receptor gene, which participates in the degradation of acetylated low-density lipoprotein. 47 Intriguingly, acetylated low-density lipoprotein promotes osteoblastic differentiation, making ZNF444 a credible candidate gene for the basis of the 19q13.43 MBD association. 48 In making these conclusions, we have made use of various differing methods for MDB, which have different sensitivity for its detecting, the differences are between, not within, cohorts and hence systematic basis is unlikely to have impacted on our findings. In addition, the association was seen in each of the four patient cohorts and was not confined to a specific MM subtype. It is important to note, that the frequency of no MBD is equivalent in each of the four case series.
Bone loss in MM is unlikely to be exclusively attributable to RANKL/RANK/OPG providing a rationale for conducting further GWAS-based analyses to identify additional MBD risk variants. It is, however, noteworthy that SNPs, other than rs4407910, which are strongly associated with BMD, were not found by us to be associated with MBD. As our power to detect an allele with an odds ratio 41.5 was high (490% power if the minor allele frequency40.2), it implies that few such alleles are likely to exist. However our power to detect an association for an odds ratio~1.2 is poor (that is, o20%), and such variants could readily exist. It is, therefore, apparent that larger studies will be required to identify additional risk loci for MBD. Although the impact of any individual genetic variant discovered by GWAS influencing MBD may be relatively modest, their identification serves to highlight genes and pathways relevant to developing novel intervention strategies. The genotyping of patients from ongoing and future clinical trials is likely to be especially informative in establishing the relationship between markers of MBD and specific therapies.
In summary, our observations provide the first evidence that germline variation influences a MM patient's risk of developing MBD and importantly impacts on this throughout disease progression. These results also provide further support for the importance of the RANK/RANKL/OPG pathway in the development of MBD. As no other established BMD locus showed an association with MBD, our results suggest the interaction between myeloma cells and bone remodeling is primarily dictated by this pathway. These findings will contribute to the development of future strategies for prevention of MBD by defining those MM patients at high risk of MBD and who may maximally benefit from therapeutic intervention.

CONFLICT OF INTEREST
GHJ has received honoraria from Celgene and Janssen-Cilag for speaking at educational meetings. GJM has participated in advisory boards for, received payment for lectures and development of educational presentations from, and has received travel support from Celgene, Novartis, Merck and Johnson & Johnson. FED has participated in advisory boards and spoken at meetings for Celgene, Ortho Biotech and Novartis, and has received travel support to attend meetings from Celgene and Ortho Biotech. JH has participated in advisory boards for Novartis, and has received travel support to attend meetings from Celgene, Janssen-Cilag and Takeda. HG has participated in advisory boards for, received payment for lectures and development of educational presentations from, and has received travel support from Celgene, Jassen-Cilag, Chugai, Novartis, Bristol-Myers Squibb, Millenium, Onyx, Amgen and Takeda. The remaining authors declare no conflict of interest.