Genome-wide Association Study Identified Chromosome 8 Locus Associated with Medication-Related Osteonecrosis of the Jaw.

Medication-related osteonecrosis of the jaw (MRONJ) is a rare but serious drug-related adverse event. To identify pharmacogenomic markers of MRONJ associated with bisphosphonate therapy, we conducted a genomewide association study (GWAS) meta-analysis followed by functional analysis of 5,008 individuals of European ancestry treated with bisphosphonates, which includes the largest number of MRONJ cases to date (444 cases and 4,564 controls). Discovery GWAS was performed in randomly selected 70% of the patients with cancer and replication GWAS was performed in the remaining 30% of the patients with cancer treated with intravenous bisphosphonates followed by meta-analysis of all 3,639 patients with cancer. GWAS was also performed in 1,369 patients with osteoporosis treated with oral bisphosphonates. The lead single-nucleotide polymorphism (SNP), rs2736308 on chromosome 8, was associated with an increased risk of MRONJ with an odds ratio (OR) of 2.71 and 95% confidence interval (CI) of 1.90-3.86 (P = 3.57*10-8 ) in the meta-analysis of patients with cancer. This SNP was validated in the MRONJ GWAS in patients with osteoporosis (OR: 2.82, 95% CI: 1.55-4.09, P = 6.84*10-4 ). The meta-analysis combining patients with cancer and patients with osteoporosis yielded the same lead SNP rs2736308 on chromosome 8 as the top SNP (OR: 2.74, 95% CI: 2.09-3.39, P = 9.65*10-11 ). This locus is associated with regulation of the BLK, CTSB, and FDFT1 genes, which had been associated with bone mineral density. FDFT1 encodes a membrane-associated enzyme, which is implicated in the bisphosphonate pathway. This study provides insights into the potential mechanism of MRONJ.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

☑ Medication‐related osteonecrosis of the jaw (MRONJ) is a drug induced adverse reaction related to antiresorptive drugs, such as bisphosphonates (BPs) and RANKL inhibitor denosumab. The mechanisms of MRONJ is unclear.

WHAT QUESTION DID THIS STUDY ADDRESS?

☑ We conducted an international genome‐wide association study meta‐analysis followed by functional analyses to identify pharmacogenomic markers associated with an increased risk of BP induced MRONJ.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

☑ We identified one locus on chromosome 8 to be significantly associated with increased risk of MRONJ in cancer patients treated with i.v. BPs and this locus was validated in patients with osteoporosis treated with oral BPs. This locus is associated with differential mRNA expression of multiple genes important to the bone‐remodeling pathways.

HOW MIGHT THIS CHANGE CLINICAL PHARMACOLOGY OR TRANSLATIONAL SCIENCE?

☑ Our study provides insights into potential mechanisms of MRONJ. If validated, our findings could provide basis for a Precision Medicine approach to antiresorptive therapy.

INTRODUCTION

Medication‐related osteonecrosis of the jaw (MRONJ) is a major adverse drug reaction related to antiresorptive drugs, such as bisphosphonates (BPs), denosumab, and antiangiogenic medications. The American Association of Oral and Maxillofacial Surgeons (AAOMS) defined MRONJ as exposure of jawbone (mandible, maxilla, or both; no history of radiation therapy) with slow healing for more than 8 weeks or no healing. The exposure of jawbone increases the risk of inflammation/infection and progress to a fracture. Other risk factors of MRONJ include older age, female gender, inflammation/infection, autoimmune disease, and diabetes. BPs are the primary agents for the treatment or prevention of osteoclast‐mediated bone loss in patients with multiple myeloma (MM), solid metastatic cancers, or osteoporosis. , MM and solid metastatic cancers, such as breast cancer, lung cancer, and prostate cancer, increase the activation of osteoclasts’ proliferation and differentiation, which commonly induce bone complications, skeletal‐related events, and bone destruction. , , Osteoporosis is a common age‐related disease, which involves bone loss and increased risk of fracture. The incidence of MRONJ ranges from 0.5–12% in patients with cancer, and 0.001–0.15% in patients with osteoporosis depending on the class of medications and doses. , ,

A previous study showed that the incidences of MRONJ vary in genetically diverse populations. We previously performed a whole‐exome sequencing analysis followed by functional studies and identified a genetic variant (rs932658) in the promoter region of the SIRT1 gene to be associated with the development of MRONJ. , Only two MRONJ genomewide association studies (GWAS) have been published to date, and both studies were conducted in patients with cancer treated with i.v. BPs. , Sarasquete et al. published the first MRONJ GWAS, which included 87 patients with MM (22 cases and 65 controls) treated with i.v. BP (pamidronate or zoledronate). This study identified a CYP2C8 single nucleotide polymorphism (SNP; rs1934951) to be associated with MRONJ (odds ratio (OR): 12.75, 95% confidence interval (CI): 3.7–43.5, Bonferroni corrected P value = 0.02). Nicoletti et al. published another MRONJ GWAS in 2012, which included 47 BP‐treated patients with breast cancer (30 with MRONJ) and an additional 1,726 healthy general population controls. SNP rs17024608 in the RBMS3 gene was identified to be significantly associated with MRONJ (OR: 5.8, 95% CI: 3.0–11.0, P = 7.47*10−8). However, the major limitations of these two studies include small number of cases, lack of medication treated control participants, and lack of replication. Furthermore, multiple subsequent independent studies have failed to replicate the association of these two loci. , There has been no GWAS published yet on MRONJ in patients with osteoporosis treated with oral BPs.

The objective of this study was to identify and replicate common genetic biomarkers associated with MRONJ in BP‐treated individuals through a GWAS approach followed by in silico functional analyses. In this GWAS with the largest number of MRONJ cases to date, discovery and replication analyses were performed in patients with cancer treated with i.v. BPs and in patients with osteoporosis treated with oral BPs, followed by meta‐analyses.

MATERIALS AND METHODS

Data collection

The participants included in this study were from multiple institutions, including the University of Florida in Gainesville, FL, USA; the National Cancer Institute (NCI) in Bethesda, MD, USA; Semmelweis University of Medical School and Dental School in Budapest, Hungary, and the University of Bologna, Italy. Additional data were obtained from the international Serous Adverse Event Consortium (iSAEC) and the Vanderbilt University Medical Center (VUMC) biobank (BioVU), which associates DNA and genomic data with de‐identified electronic health records data. The BioVU portion of the study was reviewed by the VUMC Institutional Review Board and determined to be non‐human subjects research (Table  ). All participants signed informed consent at each study site. This study was approved by the Institutional Review Board (IRB201800934) at the University of Florida.

Case and control definition

The MRONJ cases were patients who were treated with BPs and developed osteonecrosis of the jaw (ONJ) based on the criteria for MRONJ according to the AAOMS definition as exposure of jawbone (mandible, maxilla, or both; no history of radiation therapy) with slow healing for more than 8 weeks or no healing. The non‐MRONJ controls were patients who were treated with BPs for at least 12 months but did not develop ONJ. All study sites defined MRONJ cases and control based on the same AAOMS definition criteria.

Genotyping and imputation

The iSAEC samples were genotyped using Illumina HumanOmniExpressExome‐8 version 1.0 kit (Illumina, San Diego, CA), which contains ~ 1 million SNPs. The BioVU samples were genotyped using the Illumina Multi‐Ethnic Global Array (MEGA; Illumina; ~ 1.7 million SNPs). The other samples were genotyped on Illumina HumanOmni2.5S‐8 version 1 (Illumina; ~ 2.5 million SNPs) or Illumina Infinium Global Screening Array version 2 (Illumina; ~ 0.65 million SNPs; Table  ). Quality controls were performed for each of the four datasets separately. SNPs common to all GWAS panels were extracted and the four datasets were then merged. Genotype imputation was performed on the merged dataset using 1000 Genomes phase III version 5 reference panel. The detailed quality control steps for all raw data and the imputation steps were descripted in Supplementary Methods.

Genomewide association studies analyses

Discovery GWAS was performed in randomly selected 70% of the patients with cancer and replication GWAS was performed in the remaining 30% of the patients with cancer treated with i.v. BPs followed by meta‐analysis of all 3,639 patients with cancer. GWAS was also performed in 1,369 patients with osteoporosis treated with oral BPs. Multivariable logistic regression was performed to estimate the ORs and 95% CIs of each variant for the development of ONJ (yes/no) using EPACTS (http://genome.sph.umich.edu/wiki/EPACTS) adjusting for age, gender, principal components PC1‐PC3, and cancer types (in GWAS of patients with cancer). The results were filtered by effect allele frequency (EAF) > 0.01. Meta‐analysis was then performed using R (sommer package) and/or Meta Analysis Helper (METAL) to estimate the combined ORs for all independent studies. The I 2 statistics and Cochran Q test were used to assess the heterogeneity of the studies in the meta‐analyses. The random‐effects model meta‐analysis was used if the inconsistency index I 2 > 50% and fix‐effect model meta‐analysis was used if the I 2 < 50%. SNPs with P < 5*10−8 were considered as genomewide significant, whereas SNPs with P < 10−5 were considered as suggestive.

Gene‐based association analysis

Gene‐set association analysis was performed using Multi‐marker analysis of genomic annotation (MAGMA) through the Functional Mapping and Annotation of Genome‐Wide Association Studies (FUMA GWAS) platform based on our GWAS data. The MAGMA was to detect the multi‐marker effects of SNPs based on P values and linkage disequilibrium (LD) between markers using a multiple regression approach. Our input SNPs were mapped to 17,529 protein coding genes. Genomewide significance was defined at P = 0.05/17,529 = 2.852*10−6.

Enrichment analysis

Pathway and biological enrichment analyses were performed using the FUMA platform, with independent genomewide significant SNPs identified. Positional gene mapping was used to identify genes up to 10 kb from each independent significant SNP.

Topologically associating domains analysis

The 3D Genome Browser platform was used to identify the topologically associating domains of our target SNPs. The Hi‐C data was performed to identify the interaction between the target region and multiple genes. The human Mesenchymal stem cells (H‐MSCs) was selected for our Hi‐C analysis.

GWAS 4D analysis

The SNPs with P value < 5*10−8 from meta‐analysis among cancer and osteoporosis cohorts were analyzed by GWAS4D (http://mulinlab.tmu.edu.cn/gwas4d). GWAS4D detects human regulatory variants by integrative analysis of genomewide associations, chromosome interactions, and histone modifications. The top SNPs in the significant loci and suggestive loci were entered into the GWAS4D to identify functional SNPs.

In silico analyses

All SNPs of interest and the related gene were also analyzed using the UCSC genome browser, TWAS hub, ENCODE, Ensembl, regulomeDB and Open Target Platform to identify the functional significance of SNPs and genes.

RESULTS

We performed a GWAS of MRONJ in a total of 5,008 individuals of European ancestry (including 444 MRONJ cases and 4,564 controls) who were treated with BPs (oral: alendronate or risedronate; or i.v.: zoledronate or pamidronate) for at least 12 months. GWAS of the patients with cancer treated with i.v. BPs and patients with osteoporosis treated with oral BPs were performed separately. To perform MRONJ GWAS in patients with cancer treated with i.v. BPs, we randomly selected 70% of the cancer patients as discovery (217 cases and 2,208 controls) and the remaining 30% as replication (109 cases and 1,105 controls) followed by meta‐analysis in patients with cancer. We also performed a GWAS in 1,369 patients with osteoporosis treated with oral BPs included 118 MRONJ cases and 1,251 controls. Another meta‐analysis was performed to identify SNPs associated with MRONJ in patients with cancer and patients with osteoporosis combined. The study flow diagram is shown in Figure  . The characteristics of the patients are summarized in the Table  . The Manhattan plots and quantile‐quantile plots for each study are shown in Figure  .

Figure 1

Flowchart of GWAS and meta‐analyses. GWAS, genomewide association studies.

MRONJ in patients with cancer

The discovery GWAS in patients with cancer identified two loci at the suggestive level of significance (P < 10−5). The top SNP rs72817334 in an intergenic region on chromosome 2 was associated with an increased risk of MRONJ (OR: 2.27, 95% CI: 1.64–3.14, P = 7.33*10−7). The top SNP rs11948737 in an intergenic region on chromosome 5 was associated with an increased risk of MRONJ (OR: 2.76, 95% CI: 1.83–4.17, P = 1.51*10–6). However, these two SNPs were not replicated in cancer replication GWAS analysis (P > 0.05). The meta‐analysis of the discovery and replication cohort in patients with cancer identified one genomewide significant (P < 5*10−8) locus at chromosome 8p23.1 to be associated with MRONJ in patients with cancer treated with i.v. BPs. The minor allele C of the lead SNP rs2736308 on chromosome 8 located in the intronic region of C8orf12 (also known as FAM167A‐AS1) was associated with an increased risk of MRONJ with an OR of 2.71, 95% CI of 1.90–3.86 (P = 3.57*10−8; Table  , Figure  ). According to the eQTLGen and GTEx databases, rs2736308 is an expression quantitative trait locus (eQTL) for multiple genes including B lymphocyte kinase (BLK), cathepsin B (CTSB), Family With Sequence Similarity 167 Member A (FAM167A), FAM167A Antisense RNA 1 (FAM167A‐AS1; also named C8orf12), Farnesyl‐Diphosphate Farnesyltransferase 1 (FDFT1), Myotubularin Related Protein 9 (MTMR9), Nei Like DNA Glycosylase 2 (NEIL2), Solute Carrier Family 35 Member G5 (SLC35G5), L‐Threonine Dehydrogenase (TDH), and XK Related 6 (XKR6; Table  ). , , , We identified 26 LD blocks in this locus based on the LD information (r2 > 0.8) and the lead SNPs for each block are listed in Table  .

Table 1

Lead SNPs in loci with genomewide or suggestive level of significance

Chr	Position (hg19)	SNP	Ref /Alt	EA	EAF	1000 G EUR	Meta‐stage	OR	95% CI	P value	Annotated genes	Direction
8	11248500	rs2736308	C/T	C	0.34	0.35	Cancer only	2.71	1.90–3.86	3.57*10−8	BLK/FDFT1/SLC35G5/FAM167A/FAM167A‐AS1/NEIL2/MTMR9/CTSB/TDH/XKR6	**
							Cancer and osteoporosis	2.74	2.09–3.39	9.65*10−11
6	34104970	rs2029462	T/A	T	0.44	0.43	Cancer only	0.62	0.48–0.79	1.72*10−4	GRM4	**
							Cancer and osteoporosis	0.60	0.49–0.75	4.15*10−6
4	182447323	rs974578	T/C	T	0.29	0.25	Cancer only	2.14	1.55–2.96	3.93*10−6	TENM3	*
							Cancer and osteoporosis	1.79	1.35–2.36	5.01*10−5
1	108030093	rs599873	T/C	C	0.43	0.44	Cancer only	1.71	1.35–2.16	8.37*10−6	NTNG1/VAV3	*
							Cancer and osteoporosis	1.48	1.20–1.83	2.73*10−4

1000 G, 1000 Genome Projects; Alt, alternative allele; CI, confidence interval; EA, effect allele; EAF, effect allele frequency; OR, odds ratio; Ref, reference allele.

*The P value of the target SNP was < 0.05 only in cancer study.

**The P value of the target SNP was < 0.05 both in cancer and osteoporosis study.

Figure 2

Manhattan plots of GWAS meta‐analyses. (a) The meta‐analysis in patients with cancer identified one locus on chromosome 8 (C8orf12) to be significantly associated with MRONJ (P = 3.57*10−8). Two other loci (TENM3 on chromosome 4 and VAV3 on chromosome 1) were associated with MRONJ at suggestive level of significance (P < 10−5). (b) The meta‐analysis of patients with cancer and patients with osteoporosis validated the chromosome 8 locus to be significantly associated with MRONJ on a genomewide level (P = 9.65*10−11). One novel locus on chromosome 6 (GRM4) reached suggestive level of significance. GWAS, genomewide association study; MRONJ, medication‐related osteonecrosis of the jaw.

Two loci were associated with MRONJ in patients with cancer at suggestive level of significance (P < 10−5). The minor allele T of the chromosome 4 top SNP rs974578 at 4q34.3 in the intronic region of Teneurin Transmembrane Protein 3 (TENM3) was associated with an increased risk of MRONJ (OR: 2.14, 95% CI: 1.55–2.96, P = 3.93*10−6). The top SNP rs599873 on chromosome 1p13, an intergenic SNP between Netrin G1 (NTNG1) and Vav Guanine Nucleotide Exchange Factor 3 (VAV3), was associated with an increased risk of MRONJ (OR: 1.71, 95% CI: 1.35–2.16, P =8.37*10−6; Figure  , Table , Table  ).

The pharmacogenetic marker on the SIRT1 gene identified in our previous study, rs932658, was also replicated in the meta‐analysis in patients with cancer treated with i.v. BPs, with OR of 0.75, 95% CI: 0.61‐0.92 (P = 0.006).

MRONJ in patients with osteoporosis

The GWAS in the patients with osteoporosis did not yield any loci with genomewide significance, with three loci reaching suggestive level of significance ( Figure  e,f). The minor allele G of the chromosome 3 top SNP rs13072463 in the intergenic region was associated with an increased risk of MRONJ (OR: 2.99, 95% CI: 1.86–4.80, P = 6.01*10−6). The top SNP rs599873 on chromosome 5, an intronic region of RAB3C (RAB3C, Member RAS Oncogene Family), was associated with an increased risk of MRONJ (OR: 2.37, 95% CI: 1.64–3.43, P = 4.30*10−6). The top SNP rs111352217 on chromosome 7, an intergenic region SNP, was associated with an increased risk of MRONJ (OR: 7.23, 95% CI: 3.04–17.19, P = 7.60*10−6). In addition, the top SNP from the MRONJ GWAS meta‐analysis in the patients with cancer treated with i.v. BPs, rs2736308, had a P value of 6.84*10−4 in the patients with osteoporosis treated with oral BPs.

Meta‐analysis in patients with cancer and patients with osteoporosis

The meta‐analysis combining patients with cancer and osteoporosis yielded the same lead SNP rs2736308 on chromosome 8 as the top SNP (OR: 2.74, 95% CI: 2.09–3.39, P = 9.65*10−11). In addition, one new locus was identified at suggestive level of significance on chromosome 6. The lead SNP rs2029462 on chromosome 6p23.21 in the intronic region of GRM4 gene (Glutamate Metabotropic Receptor 4) was associated with a decreased risk of MRONJ (OR: 0.60, 95% CI: 0.49–0.75, P = 4.15*10−6). The two suggestive level loci in cancer meta‐analysis were not replicated in the patients with osteoporosis (P = 0.91 for rs974578 and P = 0.63 for rs599873; Figure  , Figure  , Table  ). The regional plots of the top SNPs of these loci are shown in Figure  .

Figure 3

Forest plots for the top SNPs in the meta‐analysis of patients with cancer and patients with osteoporosis. We performed the multiple logistic regression adjusted for age, gender, PCs to calculate the odds ratio (OR) of each SNP. The I 2 was calculated to estimate the heterogeneity. The random effects model meta‐analysis was used if the I 2 > 50% and fix‐effect model meta‐analysis was used if the I 2 < 50%. CI, confidence interval; MAF, minor allele frequency; OR, odds ratio; SNP, single‐nucleotide polymorphism.

Pathway and biological enrichment analyses

To identify the genes related to the top SNP, we performed gene‐based association analyses using MAGMA by the FUMA GWAS. Genes XKR6, BLK, C8orf12, AF131215.5, MSRA, SLC35G5, FAM167A, and MTMR9 on chromosome 8 were identified as genomewide significant genes associated with MRONJ with P < 2.85*10−6 (Figure  , Table  ). The FAM167A, FDFT1, BLK, MSRA, NEIL2, SLC35G5, and CTSB genes were identified as the strongest genes linked to the locus (8p23.1) by FUMA GWAS platform ( Table  ).

Figure 4

Functional mapping and annotation of genomewide association studies. XKR6, BLK, C8orf12, AF131215.5, MSRA, SLC35G5, FAM167A, and MTMR9 as the genomewide significant genes associated with the significant locus identified in the GWAS meta‐analysis. The red dash line in the plot is the genomewide significance, which was defined at P = 0.05/17529 = 2.852*10−6. GWAS, genomewide association study.

The minor allele C of the lead SNP on chromosome 8 rs2736308 was associated with lower expression of FAM167A, FDFT1, and NEIL2, but higher expression of BLK, MSRA, CTSB, and SLC35G5 based on the eqtlDirection (Table  ). The circos plot showed multiple chromatin interactions between the genomic risk locus (chr8:10803465–11401116) and genes FAM167A, FDFT1, BLK, MSRA, NEIL2, SLC35G5, and CTSB (Figure  ). The gene expression heatmap by GTEx V8 dataset showed that CTSB and FDFT1 had a higher mRNA expression compared with other related genes ( Figure  ). The 6p23.21 locus was identified to have chromatin interactions with multiple genes (Figure  ).

Figure 5

FUMA circos plot for the genomewide level locus. Chromatin interaction between genomic risk locus (blue arc) and genes were showed by orange line. The green color line linked the top SNP and the eQTL genes. (a) Circos plot for chromosome 8 top SNP. (b) Circos plot for chromosome 6 top SNP. (c) Circos plot for chromosome 1 top SNP. (d) Circos plot for chromosome 4 top SNP. eQTL, expression quantitative trait locus; FUMA, Functional Mapping and Annotation; SNP, single‐nucleotide polymorphism.

The lead SNP rs599873 of chromosome 1 locus is an intergenic SNP between NTNG1 and VAV3 genes. No interaction between this locus and genes was identified by chromatin interactions analysis (Figure  ). The lead SNP rs974578 of chromosome 4 locus interacted with the DCTD gene (Figure  ).

Topologically associated domains analysis

The 3D Genome browser , predicted that there are seven topologically associating domains surrounding the chromosome 8 locus region (chr8:10803465–11401116; Figure  ). This locus included the peak of H3k4me3 (a marker associated with the activation of transcription of nearby genes), the peak of H3K4me1 (an enhancer mark), and the peak of H3K9me3 and H3K27me3 (two marks associated with heterochromatin). The circular chromosomal conformation capture (4C) plot using 3D Genome browser and the UCSC genome browser further corroborates that there are long‐range promoter enhancer interactions between the region around the lead SNP and other genes, such as MTMR9, BLK, and FAM167A in H‐MSC (Figure  ). The Hi‐C data also showed that the chromosome 6 locus (chr6:34092806–34109404) was associated with the peaks of H3k9me3, and H3k27me3 in H‐MSC. The 4C plot data using 3D Genome browser and the UCSC genome browser showed SNP rs2029462 was interacted with multiple genes, including the GRM4 gene (Figure  ).

Figure 6

Summary of Hi‐C data of the top locus on chromosome 8. Seven Topologically associated domains (TADs) surrounding the chromosome 8 locus region (chr8:10803465–11401116) (a). Circular chromosomal conformation capture (4C) plot measures the interaction frequencies between rs2736308 and other loci (peak signal). The plot showed that there is a chromatin interaction event between rs2736308 and MTMR9, FAM167A, and BLK in human mesenchymal stem cell (H‐MSC). The anchoring point is the interesting SNP position (b). SNP, single‐nucleotide polymorphism.

The Hi‐C data of the chromosome 4 locus (chr4:182367808–182450209) reported a CpG island, the peak of H3M27me3, and H3K4me3 ( Figure  c,d). There was no evidence that this locus had any significant chromatin interaction with protein‐coding genes. The top SNP of chromosome 1 locus had interactions with the NTNG1 and VAV3 gene (Figure  e,f). The 4C data further corroborates that there are long‐range promoter enhancer interactions between the region around the lead SNP and genes NTNG1, and VAV3 in H‐MSC (Figure  ).

GWAS4D

GWAS4D (http://mulinlab.tmu.edu.cn/gwas4d) is an in silico tool that evaluates GWAS results and identifies potential regulatory variants by integrating the latest multidimensional functional genomics resources and the published algorithms by Mulinlab. The lead SNP rs34128921 of block 3 (OR: 2.33, 95% CI: 2.02–2.70, P = 3.60*10−9) on chromosome 8 has a significant Hi‐C interaction with the BLK gene (Figure  ). The GWAS4D predicted that this SNP was associated with decreasing transcription factor binding affinity (THAP1, MLLT1, and CUX1; Figure  , Table  ). SNP rs76154097 in block 15 (the lead SNP is rs77072957) on chromosome 8 also has a significant Hi‐C interaction with the BLK gene (Figure  c). SNPs rs1600252 (Figure  d) in block 8 showed a significant Hi‐C interaction with the CTSB gene. SNPs rs2248699 (Figure  e) in block 20. GWAS4D analysis also identified the SNP rs11753307 to be associated with GRM4 gene regulation (Figure  f). The GWAS4D of chromosome 6 locus showed this locus to have Hi‐C interaction with multiple genes.

Phenome Wide Association Studies

The phenome wide association studies (PheWAS) of the lead SNPs using the UK Biobank dataset based on the data of Neale lab (http://www.nealelab.is) showed alternative allele T of the top SNP rs2736308 was associated with higher bone mineral density (BMD; Figure  ).

Transcriptome‐wide association study analysis

The transcriptome‐wide association study (TWAS) identified the BLK gene to be negatively associated with heel BMD in blood tissue (Table  ). The FDFT1 gene was positively correlated with heel BMD in blood tissue (Table  ). The CTSB gene was negatively associated with heel BMD (Table  ).

GeneMANIA pathway analysis

All the above in silico analyses identified BLK, CTSB, FDFT1, GRM4, NTNG1, and VAV3 genes to be important genes associated with MRONJ. The GeneMANIA pathway analysis showed that the VAV3, BLK, FDFT1, and CTSB genes interact through epidermal growth factor receptor (EGFR) and SLC25A24. These genes are associated with small GTPases or initial immune system pathway (Figure  ).

DISCUSSION

In this GWAS meta‐analysis that included the largest number of MRONJ cases to date, we identified and validated pharmacogenomic markers of MRONJ in patients with cancer treated with i.v. BPs. We also performed the first MRONJ GWAS in patients with osteoporosis treated with oral BPs. One locus on chromosome 8 was identified to be significantly associated with MRONJ at the genomewide level in patients with cancer treated with i.v. BPs and was validated in patients with osteoporosis treated with oral BPs. The minor allele C of the top SNP rs2736308 was associated with an increased risk of MRONJ with an OR of 2.74 (CI, 2.09–3.39, P = 9.65*10−11) in the combined meta‐analysis of patients with cancer and patients with osteoporosis. The in silico analyses and pathway analyses provided evidence that the minor allele C of rs2736308 was associated with lower expression of FDFT1 and CTSB and higher expression of BLK gene.

The FDFT1 gene encodes squalene synthase which is a membrane‐associated enzyme in the mevalonate pathway. This pathway generates cholesterol in cells, and side pathways produce farnesyl and geranylgeranyl groups that are attached to small GTPase, crucial regulators of membrane trafficking, and the cytoskeleton. , BPs are incorporated into bone, and released as bone is resorbed by osteoclasts and are internalized into osteoclasts. The nitrogen‐containing BPs, such as zoledronate and pamidronate, inhibit farnesyl diphosphate synthase, a key enzyme in the mevalonate pathway. This prevents the synthesis of geranylgeranyl diphosphate in the downstream of the mevalonate pathway and blocks the prenylation of small GTPases. This disrupts normal osteoclast function and induces apoptosis of osteoclast cells. , Lower expression of squalene synthase might decrease the competitive inhibition of the activation of small GTPases and thus modify the effects of nitrogen‐containing BPs on osteoclasts. This could alter the propensity for the development of MRONJ. Our proposed mechanism for MRONJ based on the observed association is presented in Figure  .

The BLK gene encodes a nonreceptor tyrosine‐kinase of the src family of proto‐oncogenes. This kinase stimulates the insulin synthesis and secretion in response to glucose and also increases the expression of several pancreatic beta‐cell transcription factors. BLK deficiency was shown to be associated with diabetes, which is a risk factor of MRONJ. The BLK gene was also associated with inflammation/infection responses by stimulating the secretion of metalloproteinase‐9 (MMP‐9), which was reported to be associated with rheumatoid arthritis, atherosclerosis, cancer, diabetes, obesity, and osteoporosis. AAOMS have reported the inflammation/infection as a risk factor of MRONJ.

The CTSB gene has also been observe to be involved in this macrophage‐mediated inflammation/infection responses. Álvaro de Mingo et al. reported that the CTSB gene regulates the inflammation responses by regulating the expression of SIRT1, which was identified to be associated with MRONJ. , ,

Although the GRM4 locus on chromosome 6 and VAV3 locus on chromosome 1 did not reach genomewide significance, we believe these two loci merit discussion. The GRM4 gene on chromosome 6 is expressed in bone cells, modulated osteoblast differentiation, and activity, and plays an important role in bone remodeling and homeostasis. The deficiency in GRM4 was shown to increase the risk of cancer metastasis to the bone and cancer bone metastasis was associated with higher risk of MRONJ. The VAV3 gene encodes the Vav protein, which belongs to the family of mammalian Rho Guanine nucleotide exchange factors. The Rho Guanine nucleotide exchange factors were associated with activation of the Rho GTPases. Rho family of GTPases plays an essential role in osteoclast differentiation. VAV3 is most abundant in osteoclasts during differentiation, while the mature osteoblasts have low VAV3 expression levels. VAV3 overexpression in osteoclasts stimulates osteoclasts differentiation and bone resorption, leading to bone loss.

Our pathway analysis result suggests that VAV3, BLK, FDFT1, and CTSB interact with each other through EGFR and SLC25A24. Karin Writzl et al. demonstrated that SLC25A24 regulated mitochondrial ATP synthesis and hyperpolarization. EGFR plays an essential role in osteoblast cell proliferation, differentiation, and survial. The osteoblast‐specific deletion of EGFR was associated with bone defects and low expression or lack of EGFR induced low BMD.

The PheWAS analysis for the target SNPs across all UK Biobank phenotypes released by Neale and colleagues indicated that the alternative allele T of the lead SNP rs2736308 was associated with an increased BMD. So, the minor allele C was associated with a decreased BMD of heel bone, which is positively associated with jawbone BMD (Figure  ). This is consistent with the TWAS result. Therefore, it is plausible that the minor allele C of the lead SNP rs2736308 on chromosome 8 is associated with lower BMD by regulating the expression of FDFT1, CTSB, and BLK, which may influence BMD through interaction with EGFR.

Our study has some limitations that need to be recognized. First, due to the imbalance of case and controls across different datasets, we combined all patients with cancer and randomly selected 70% for discovery and 30% for replication. Second, different datasets were genotyped using different GWAS platforms and only 126,146 variants were common across different platforms and were included in the GWAS genotype imputation. It is not clear how this might have affected the imputation quality. Last, there were not enough patients of other ancestries, we only focused on patients of European ancestry. Therefore, our findings might not be generalizable to individuals of other ancestries.

In summary, we identified a genomewide significant locus on chromosome 8 to be associated with higher odds of MRONJ in patients with cancer treated with i.v. BPs and this locus was validated in patients with osteoporosis treated with oral BPs. This was made possible by the unprecedented size of this international study. This locus is associated with differential mRNA expression of multiple genes important to the bone‐remodeling pathways. We also identified three loci at suggestive level of significance that may merit further investigation. Our study provides insights into potential mechanisms of MRONJ. If validated, our findings could provide basis for a Precision Medicine approach to antiresorptive therapy.

FUNDING

Y.G. is supported by University of Florida College or Pharmacy, University of Florida Health Cancer Center, University of Florida Informatics Institute, and University of Florida Clinical and Translational Science Institute, which is supported in part by the NIH National Center for Advancing Translational Sciences under award number UL1 TR001427. Some of the dataset(s) used for the analyses described were obtained from Vanderbilt University Medical Center’s BioVU which is supported by numerous sources: institutional funding, private agencies, and federal grants. These include the NIH‐funded Shared Instrumentation Grant S10RR025141; and CTSA grants UL1TR002243, UL1TR000445, and UL1RR024975. Genomic data are also supported by investigator‐led projects that include U01HG004798, R01NS032830, RC2GM092618, P50GM115305, U01HG006378, U19HL065962, R01HD074711; and additional funding sources listed at https://victr.vumc.org/biovu‐funding/.

CONFLICT OF INTEREST

The authors declared no competing interests for this work.

AUTHOR CONTRIBUTIONS

G.Y. wrote the manuscript. G.Y., S.S., J.L., C.M., I.H., L.S.H., D.W., T.L., J.K., P.L., B.B., J.K., S.V.D., D.P., W.F., G.P., J.M., and Y.G. performed the research. G.Y., C.M., and Y.G. analyzed the data. Y.G. and G.Y. designed the research.

DISCLAIMER

As an Associate Editor of Clinical Pharmacology & Therapeutics, Dr. Sara L. Van Driest was not involved in the review of decision process for this paper.

Supplementary Material

Click here for additional data file.