Identification of Potentially Pathogenic Variants Associated with Recurrence in Medication-Related Osteonecrosis of the Jaw (MRONJ) Patients Using Whole-Exome Sequencing

Background: Bisphosphonates are antiresorptive and antiangiogenic drugs that prevent and treat bone loss and mineralization in women with postmenopausal osteoporosis and cancer patients. Medication-related osteonecrosis of the jaw (MRONJ) is commonly caused by tooth extraction and dental trauma. Although genetic and pathological studies about MRONJ have been conducted, the pathogenesis of MRONJ still remains unclear. Methods: We aimed to identify genetic variants associated with MRONJ, using whole-exome sequencing (WES). Ten MRONJ patients prescribed bisphosphonates were recruited for WES, and jawbone tissue and blood samples were collected from the patients. Results: The analysis of the WES data found a total of 1866 SNP and 40 InDel variants which are specific to MRONJ. The functional classification assay using Gene Ontology and pathway analysis discovered that genes bearing the MRONJ variants are significantly enriched for keratinization and calcium ion transport. Some of the variants are potential pathogenic variants (24 missense mutations and seven frameshift mutations) with MAF < 0.01. Conclusions: The variants are located in eight different genes (KRT18, MUC5AC, NBPF9, PABPC3, MST1L, ASPN, ATN1, and SLAIN1). Nine deleterious SNPs significantly associated with MRONJ were found in the KRT18 and PABPC3 genes. It suggests that KRT18 and PABPC3 could be MRONJ-related key genes.


Introduction
Medication-related osteonecrosis of the jaw (MRONJ) refers to an uncommon medical condition that can occur in patients receiving antiresorptive or antiangiogenic therapy to treat osteoporosis and diverse malignant conditions such as multiple myeloma, breast cancer, and prostate cancer [1]. When the first case of MRONJ was documented in 2003, only bisphosphonates were associated with the pathophysiology, and thus the condition was widely known as bisphosphonate-related osteonecrosis of the jaw (BRONJ) [2]. However, in light of the numerous recent reports of osteonecrosis cases associated with the use of various other antiresorptive and antiangiogenic agents such as denosumab, the current universal nomenclature for the disease has been modified to MRONJ. The antiangiogenic agents or targets of the VEGF pathway, such as the anti-VEGF monoclonal antibodies (Bevacizumab) and tyrosine kinase inhibitors (Sunitinib, Sorafenib), that are prescribed for the treatment of cancer-related conditions or organ rejection in renal transplant patients (Sirolimus) have also been indicated in the pathophysiology of the disease [3]. The key hypotheses of rejection in renal transplant patients (Sirolimus) have also been indicated in the pathophysiology of the disease [3]. The key hypotheses of disease specificity for the jaws include bone remodeling inhibition, angiogenesis inhibition, inflammation, innate or acquired immune dysfunction, soft tissue toxicity, and genetic predisposition. Several hypotheses can explain the overall pathophysiology of MRONJ [4,5]. Antiresorptive and antiangiogenic reagents, including bisphosphonates and denosumab, inhibit bone remodeling by affecting osteoclast formation and differentiation, or they directly inhibit angiogenesis by reducing blood flow to the bone [5,6]. Genetic factors such as single nucleotide polymorphisms (SNPs) and mutations are significantly associated with MRONJ development by regulating bone remodeling, angiogenesis, and immune responses [7,8]. Prior attempts to determine the efficacy of a drug holiday in preventing MRONJ have yielded inconclusive results, and thus the evidence to support or refute the use of such a holiday remains limited [7]. The hallmark clinical feature of MRONJ is the presence of exposed necrotic bone in the oral cavity that persists for more than eight weeks without any signs of healing. This is frequently accompanied by concurrent infection signs such as pain, swelling, fistula formation, and in its severest form-pathological fractures ( Figure 1). In some MRONJ patients, a specific symptom, known as the 'numb chin syndrome' may also manifest [8]. According to the current AAOMS guidelines, MRONJ can be classified into four distinct stages (stages 0 to III), for which the suggested methods of treatment differ. For early stages (stage I) of MRONJ patients, conservative treatment, including antibiotic and septic administration, and an antibacterial oral rinse, results in pain control and prevention of necrosis. Surgery may not be needed in the absence of disease progression, and if the patient has an adequate quality of life [4]. For patients suffering from advanced stage (stage II-III) MRONJ, invasive surgery, including jaw resection, can yield optimal mucosal healing; however, up to 25% of patients can suffer from severe postoperative morbidity following this type of extensive surgery. Therefore, to minimize morbidity by establishing novel prevention or early detection methods, recent studies have focused on determining the pathophysiology of MRONJ.  The incidence of MRONJ is relatively lower in the oral bisphosphonate group for osteoporosis than in the intravenous bisphosphonate group for cancer; hence, to date, many hypotheses concerning the pathogenesis of MRONJ have been raised by many authors. Previous studies reported that genetic polymorphisms in several candidate genes (TGFb1, MMP2, CYP2C8, VEGF, COL1A1, RANK, OPG, OPN, and PPARG) involved in the activation and differentiation of osteoclast and bone remodeling are significantly associated with MRONJ [1,[9][10][11]. In addition, many researchers recently found that MRONJ is caused by multiple biological pathways involving immunoglobulin, cell adhesion, and cytoskeleton [12]. However, despite the numerous genetic studies, such as the genome-wide association studies (GWAS), whole-exome studies (WES), and the candidate gene studies (CGS), that have been performed to ascertain the condition's exact pathogenic mechanism, the pathogenesis remains yet to be determined. Furthermore, although potential associations between the development of MRONJ and several genes have been documented in the previous literature, the GWAS, WES, and CGS have all failed to elucidate a single gene as the sole risk factor for MRONJ development.
Our study aimed to identify genetic variants associated with MRONJ by analyzing somatic and germline mutations through WES. By conducting an integrated analysis with the preceding research results, eight genes were determined to be closely related to MRONJ, and mutations in the genes were investigated in the present study. The result found that twelve and ten genetic mutations of KRT18 and PABPC3 genes, respectively, are predicted to modify the protein structure of the genes in MRONJ patients.

Ethics Statement
Experiments and analyses involving human subjects or human data in this study were approved by the Institutional Review Board of Dankook University Hospital (IRB No. 2018-01-009). All clinical investigations were performed in accordance with the Declaration of Helsinki. Written informed consent was obtained from all patients included in this study.

Patient Selections
Ten MRONJ patients (stage II-IIII) from the oral and maxillofacial surgery department of Dankook University Hospital participated in this study. The patients visited the Dankook University Hospital for surgical intervention from March 2018 to December 2018 and received bisphosphonates (BPs) therapy, either risedronate (9 cases) or alendronate (1 case). Patients were considered to have MRONJ if all of the following three characteristics were present according to the position paper issued by the American Association of Oral and Maxillofacial Surgeons (AAOMS) in 2014; (i) current, previous treatment with an antiresorptive or antiangiogenic therapy, (ii) non-healing, exposed bone in the maxillofacial region that had persisted for more than eight weeks, and (iii) no history of radiation therapy to the jaw. We also analyzed the characteristics of the clinical findings (medications, smoking, sex, stage) in a panoramic radiographic study ( Figure 1).

DNA Preparation
Peripheral venous blood (2 mL) and routine jawbone specimens were obtained from the ten MRONJ patients in the Department of Oral and Maxillofacial Surgery of Dankook University Hospital. Genomic DNA was extracted from the blood and the lesion jawbone tissues immediately after surgical treatment of ten MRONJ patients using a DNeasy Blood & Tissue Kit (Qiagen, Germany) in accordance with the manufacturer's protocols, and stored at −80 • C until use.

Whole-Exome Sequencing (WES) and Individual Variant Calling
For NGS library construction, each genomic DNA (1 µg) was sheared to 180-280 bp fragments using a DNA shearing system (Covaris, Woburn, MA, USA). Sequencing libraries were generated using a SureSelect Human All Exon Kit (Agilent Technologies, Santa Clara, CA, USA) following the manufacturer's instructions. DNA that was ligated with adapters were enriched during PCR reaction. The WES libraries were captured and purified using streptavidin-coated magnetic beads and the AMPure XP system (Beckman Coulter, Beverly, MA, USA). WES was performed on the Illumina NovaSeq 6000 with a 2 × 150 bp paired end (PE) read.
Sequencing reads were mapped to the human reference genome (GRCh37/hg19) using the Burrows-Wheeler Alignment tool (BWA 0.7.8) with default parameters [13]. PCR duplicates were removed using Picard (http://picard.sourceforge.net/index.shtml, accessed on 8 November 2019). Variant calling was performed using the Genome Analysis Tool Kit (GATK v3.8) [14]. The functional annotation of variants was conducted using ANNOVAR (2015Dec14) [15]. To obtain genetic variants highly related to MRONJ, we compared two datasets with the variant candidates (Korean variant archive; KOVA and GSK project). KOVA was used as the control dataset to exclude the common Korean variants. The GSK project investigating the MRONJ-related variants found in the saliva of MRONJ patients was used to collect variants shared between MRONJ diagnosed patients.

Functional Gene Classification and Pathway Analysis
To investigate the biological relevance of the candidate genes, we performed gene ontology (GO) enrichment analysis to classify them into three classes, including biological process (BP), cellular component (CC), and molecular function (MF), by using Metascape software (https://metascape.org/gp/index.html, accessed on 19 February 2020) [16]. Then, the enrichment of candidate genes was conducted to identify pathways using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/ pathway.html, accessed on 19 February 2020). The significant function of genes was filtered by statistical cut-off using the adjusted p-value (FDR; False Discovery Rate).

Clinical Findings
The clinical characteristics of the ten MRONJ patients are shown in Table 1. One man and nine women were enrolled in this study, and none of them were smokers. The mean age was 77.8 ± 5.64 years. Of the ten MRONJ patients, the mandible bone was involved in seven (70%), and the rest were maxillary bones. One of the ten MRONJ patients was spontaneous, eight had teeth extractions, and the other had root canal therapy prior to the onset of MRONJ. Only one patient was taking alendronate, and the others were taking risedronate (Table 1). Ten patients developed the disease after long-term administration for at least four years.

Whole-Exome Sequencing and Quality Controls
WES data were generated with an average 11.4 Gb of sequences per sample on the Illumina NovaSeq6000 platform. Sequencing quality for Q20 and Q30 was an average of 96.75% and 91.99%, respectively (Supplementary Table S1). Of the WES reads, 97.34% were mapped to the human reference genome. The blood and lesion samples (i.e., the lesion jawbone tissues) were sequenced to an average of 104.16× depth on target, and an average of 99.57% of exonic regions were covered with at least 10× depth on target (Supplementary Table S1).

Identification of Single Nucleotide Polymorphism (SNPs) Related to MRONJ
From WES data, an average of 433,602 germline SNP variants and 367,806 somatic SNP variants were identified in ten patients (Supplementary Table S2). To identify variants associated with MRONJ, we used a strict filtering method that collects common variants from all ten patients. We identified a total of 39,720 germline SNP variants and 32,225 somatic SNP variants shared by all ten patients. Then, to rule out the common variants in Korean, we compared the total variants with the KOVA dataset containing 293,049 variants in 1055 healthy Korean individuals [19]. The 219,772 MRONJ-associated variants discovered from the GSK project with 142 saliva samples of MRONJ patients [20] were used for variant recalibration and separation of somatic mutations from germline mutations ( Figure 2). Through the analysis, we identified a total of 4171 SNPs in 2321 genes associated with MRONJ, including 3876 SNPs detected in all three datasets (saliva, blood, and lesion), 216 SNPs detected in both blood and saliva datasets, and 79 SNPs detected in both saliva and lesion datasets. Additionally, we collected 1151 lesion-specific SNPs to identify somatic variants (Table 2, Figure 2). We found 1020 nonsynonymous SNP variants, including the missense (n = 946), nonsense (n = 8), and unknown mutations (n = 66) in the exonic and splicing region (Table 3, Supplementary Table S3).

Identification of Insertion/Deletions (InDels) Variants Related to MRONJ
As for InDel variants, we identified an average of 72,754 germline InDels and 72,715 somatic InDels from WES data of ten patients (Supplementary Table S2). A total of 4901 germline InDels and 3892 somatic InDels are common in all ten MRONJ patients ( Figure 2). Finally, a total of 1641 InDels in 1142 genes were only discovered in MRONJ patients, except for common Korean variants as we described above, which include 265 InDels detected in all three datasets (saliva, blood, and lesion), 241 InDels detected in both saliva and lesion dataset, and 46 variants detected in both saliva and blood dataset. Additionally, we included 1089 lesion-specific InDels, which were somatic variants. Of these, we collected 132 exonic InDel variants to identify variants that affect MRONJ, including common InDels in the saliva and blood (n = 4), common InDels in the saliva and lesion (n = 42), common InDels in the saliva, blood, and lesion (n = 40), and lesion-specific InDels (n = 46). They contain 32 frameshift InDels, 53 non-frameshift InDels, 46 unknown mutations, and one nonsense mutation (Table 3, Supplementary Table S4).

Enrichment and Pathway Analysis
To understand the functional relevance of the isolated genetic variants focusing on MRONJ, we performed Gene Ontology (GO) and Kegg pathway analyses on 665 genes containing 1020 nonsynonymous SNP variants. As shown in Figure 3, the most significant gene set enrichment in the biological processes (BP) was 'detection of stimulus', 'keratinization', and 'calcium ion transport'. Interestingly, the 'protein kinase C signaling pathway' was also elected in the biological processes ( Figure 3; Supplementary Table S5). As referred to in the previous study, the bisphosphonates inhibit protein kinase C signaling, resulting in the accumulation of calcium and degradation of osteoclast [33]. The most significant enriched function was 'ECM structural constituent' for the molecular function (MF), followed by 'olfactory receptor activity', and 'calcium ion binding' (Figure 3; Supplementary Table S5). For the cellular component (CC), the function of genes associated with 'keratin filament' was most significantly enriched ( Figure 3; Supplementary Table S5). For the Kegg pathway, the top-ranked significant pathways were 'olfactory transduction', 'autoimmune thyroid disease', and 'ECM-receptor interaction' (Supplementary Table S5). One of the primary mechanisms, the extracellular matrix (ECM), comprises over 90% of type I collagen, the most abundant collagen, in human bone. It is a quite meaningful result because Type I collagen is synthesized by osteoblasts, which drives bone tissue remodeling during wound healing.
Furthermore, we examined the interaction between the classified functions of MRONJspecific variants by function network analysis using ClueGO (p-value < 0.05; Benjamini-Hochberg). The integrative function network analysis showed that the eight parent functions, including 'integral component of luminal side of endoplasmic reticulum membrane', 'O-glycan processing', 'sensory perception of the chemical stimulus', 'keratin filament', 'calcium ion transport', 'protein kinase A signaling', 'phenylpropanoid metabolic process', and 'post-synaptic specialization organization' were mainly clustered with 665 genes (Supplementary Figure S2).

Pathogenicity Analysis
For the private or rare variants (MAF < 0.01), we also utilized three databases, including the 1000 genomes project database, NHLBI exome sequencing project (ESP), and ExAC database. We filtered in the candidate variants with allele frequencies lower than 0.01 or no allele frequency in all three databases. A total of 24 candidate SNPs and seven candidate InDels were selected in four genes (KRT18, MUC5AC, NBPF9, and PABPC3) and five genes (MST1L, ASPN, ATN1, PABPC3, and SLAIN1), respectively (Table 4). associated with 'keratin filament' was most significantly enriched (Figure 3; Supplementary Table S5). For the Kegg pathway, the top-ranked significant pathways were 'olfactory transduction', 'autoimmune thyroid disease', and 'ECM-receptor interaction' (Supplementary Table S5). One of the primary mechanisms, the extracellular matrix (ECM), comprises over 90% of type I collagen, the most abundant collagen, in human bone. It is a quite meaningful result because Type I collagen is synthesized by osteoblasts, which drives bone tissue remodeling during wound healing. Figure 3. GO functional analysis of candidate genes. GO enrichment analysis of WES genes was retrieved using Metascape software. Significantly (p < 0.05) enriched GO terms involved in biological process (pink), cellular component (blue), and molecular function (yellow) and branches are presented. The adjusted statistically significant values were negative 10-base log-transformed. The X-axis represents the enrichment scores of the terms −log10(p-value), and the Y-axis represents the enriched GO terms in the pathway. P-value was adjusted by Benjamini-Hochberg FDR. GO, gene ontology. The adjusted statistically significant values were negative 10-base log-transformed. The X-axis represents the enrichment scores of the terms −log10(p-value), and the Y-axis represents the enriched GO terms in the pathway. P-value was adjusted by Benjamini-Hochberg FDR. GO, gene ontology. In addition, we analyzed 24 rare SNP variants for the annotation of potentially deleterious effects using SIFT and polyphen2. Of these, three SNP variants (p.G38C, p.G43R in KRT18, and p.K231E in PABPC3) were predicted to be 'deleterious' in the SIFT, and 'possibly damaging or probably damaging' in Polyphen-2 (Table 4). Six SNP variants (p.R27W, p.P28Q, p.A32S, p.S34R in KRT18, and p.A181T, p.L207F in PABPC3) were predicted to be either deleterious or benign in the SIFT and Polyphen-2 (Table 4). KRT18 encoding the type 1 intermediate filament chain regulates protein synthesis, organelle positioning, and protection from apoptosis and necrosis. PABPC3 is a member of RNA binding proteins involved in the cytoplasmic regulatory process of mRNA metabolisms, such as RNA stability and translation initiation. RBMS3 (the c-myc gene single-strand binding protein), a paralog for PABPC3, regulates myc/ras cooperative transforming activity, apoptosis, DNA replication, and cytoplasmic RNA metabolism. GWAS studies showed that two SNPs (rs17024608 and rs10510628) in the RBMS3 were associated with a bisphosphonates-induced reduction in collagen formation, osteonecrosis disruption, and bone turnover [34,35].
To predict protein structural alteration by amino acid (AA) changes, protein conformational changes in WT and mutants of KRT18 and PABPC3 were predicted using the AlphaFold protein structure database, phyre2 web tool. Nine rare SNP variants (p.R27W, p.P28Q, p.A32S, p.S34R, p.G38C, and p.G43R in KRT18, and p.A181T, p.L207F, and p.K231E in PABPC3) were structurally changed ( Figure 4A). In addition, we calculated the distance between the atoms of two residues to predict the effect of the mutations on the distances between each pair of residues using the UCSF chimera tool. The distances between wild-type and mutations were changed in KRT18 and PABPC3 ( Figure 4B). Notably, we identified pairs of charged residues that were significantly skewed from 8.68 Å to 4.67 Å by the R27W mutation ( Figure 4B).
Furthermore, we used DynaMut2 to analyze the entropy energy change between wild-type and mutant structures (in kcal/mol/K). All six substitutions (p.R27W, p.P28Q, p.A32S, p.S34R, p.G38C, and p.G43R), located in the head domain, which is a primary phosphorylation site of KRT18, have weak effects (−0.5 < ∆∆G < 0.5 kcal/mol) on KRT18 stability ( Figure 5). The phosphorylation of KRT18 at serine 33 and 34 regulate filament stability, keratin reorganization, binding to 14-3-3 protein during mitosis, and keratin protein turnover during apoptosis [36][37][38]. The glycosylation at serine 30 and 31 protects epithelial tissue from injury [37]. The loss and mutations of KRT18 contribute to the phenotypes of several human diseases, including skin disorders, chronic renal disease, and liver cirrhosis [38].
In addition, three substitutions (A181T, L207F, and K231E) are located in the RNA recognition motifs (RRMs) of PABPC3 that bind the poly(A)tail of single-stranded RNAs ( Figure 4A). Unfortunately, the entropy change of PABPC3 could not be analyzed because the intact protein structure of PABPC3 could not be obtained from the protein data bank (PDB). However, a p.K231E substitution located in the RNP1 motif ((K/R)-G-(F/Y)-(G/A)-F-V-X-(F/Y)), a highly conserved region, may interact with the translation initiation factor and regulatory proteins, affecting translation initiation and nonsense-mediated decay. The adjacent amino acid is noted by black words, and the predicted distance between two residues in the protein structure is noted by orange words. Furthermore, we used DynaMut2 to analyze the entropy energy change between wild-type and mutant structures (in kcal/mol/K). All six substitutions (p.R27W, p.P28Q, p.A32S, p.S34R, p.G38C, and p.G43R), located in the head domain, which is a primary phosphorylation site of KRT18, have weak effects (−0.5 < ΔΔG < 0.5 kcal/mol) on KRT18 stability ( Figure 5). The phosphorylation of KRT18 at serine 33 and 34 regulate filament stability, keratin reorganization, binding to 14-3-3 protein during mitosis, and keratin protein turnover during apoptosis [36][37][38]. The glycosylation at serine 30 and 31 protects epithelial tissue from injury [37]. The loss and mutations of KRT18 contribute to the The adjacent amino acid is noted by black words, and the predicted distance between two residues in the protein structure is noted by orange words.
In addition, three substitutions (A181T, L207F, and K231E) are located in the RNA recognition motifs (RRMs) of PABPC3 that bind the poly(A)tail of single-stranded RNAs ( Figure 4A). Unfortunately, the entropy change of PABPC3 could not be analyzed because the intact protein structure of PABPC3 could not be obtained from the protein data bank (PDB). However, a p.K231E substitution located in the RNP1 motif ((K/R)-G-(F/Y)-(G/A)-F-V-X-(F/Y)), a highly conserved region, may interact with the translation initiation factor and regulatory proteins, affecting translation initiation and nonsense-mediated decay.

Discussion
Bisphosphonates, known as a primary drug for the onset of MRONJ, is the antiresorptive agent that prevents loss of bone density. They are used to treat osteoporosis, osteogenesis imperfecta, primary hyperparathyroidism, and cancer. The bisphosphonates bind preferentially to calcium ions in the bones and accumulate in high concentrations. However, the bisphosphonates cause mucosal ulceration and exposure to the necrotic bone, leading to osteonecrosis of the jaw. Previous studies reported that the bisphosphonates induce senescence and apoptosis of the cells by blocking the cholesterol biosynthetic pathways (mevalonate pathway) that indirectly inhibit the blood supply of the bone [20,39,40]. Bisphosphonates, which are closely related to anti-angiogenesis, affect the structure of soft tissues in the jaw, causing soft tissue cell damage. In this study, we identified the candidate genes and sorted out the variants closely associated with MRONJ in lesion and blood samples through WES. Although previous studies have constantly attempted to screen the candidate genes and variants that show a significant correlation with MRONJ using GWAS, WES, and CGS, the single genes and variants associated with MRONJ development have not been identified; therefore, the need for in-depth research on the onset, development, diagnosis, and treatment of MRONJ is drawing attention by securing various genetic data according to various races, drugs, living environments, and clinical aspects.
In function enrichment analysis for MRONJ-specific variants, the protein kinase C (PKC) signaling pathway, that can either read out lipid signals alone or combine the

Discussion
Bisphosphonates, known as a primary drug for the onset of MRONJ, is the antiresorptive agent that prevents loss of bone density. They are used to treat osteoporosis, osteogenesis imperfecta, primary hyperparathyroidism, and cancer. The bisphosphonates bind preferentially to calcium ions in the bones and accumulate in high concentrations. However, the bisphosphonates cause mucosal ulceration and exposure to the necrotic bone, leading to osteonecrosis of the jaw. Previous studies reported that the bisphosphonates induce senescence and apoptosis of the cells by blocking the cholesterol biosynthetic pathways (mevalonate pathway) that indirectly inhibit the blood supply of the bone [20,39,40]. Bisphosphonates, which are closely related to anti-angiogenesis, affect the structure of soft tissues in the jaw, causing soft tissue cell damage. In this study, we identified the candidate genes and sorted out the variants closely associated with MRONJ in lesion and blood samples through WES. Although previous studies have constantly attempted to screen the candidate genes and variants that show a significant correlation with MRONJ using GWAS, WES, and CGS, the single genes and variants associated with MRONJ development have not been identified; therefore, the need for in-depth research on the onset, development, diagnosis, and treatment of MRONJ is drawing attention by securing various genetic data according to various races, drugs, living environments, and clinical aspects.
In function enrichment analysis for MRONJ-specific variants, the protein kinase C (PKC) signaling pathway, that can either read out lipid signals alone or combine the ability to read out simultaneous calcium ions, is highly enriched in the biological processes. The calcium ion is a second messenger that regulates differentiation and activation of osteoclasts according to calcium concentration [41]. The calcium ion signaling regulates the activation of PKC and NFATc1 in osteoclasts, which is involved in cell proliferation and differentiation [33]. Previous studies reported that treatments of bisphosphonates, such as alendronate and risedronate, were accumulated in the bone matrix. They trigger osteoclast apoptosis, resulting in cytoskeletal regulation [42,43]. Human parathyroid hormone (PTH) is highly associated with the regulation of calcium metabolism; thus, when calcium ions increase in the serum, the osteoclast activates, which degrades bone. These results suggest that the calcium ion is closely related to MRONJ development, regulating cell proliferation and differentiation. Through our study, we were able to observe variations in the following regulatory factors: FLT4, MC1R, ULK4, ADGRV1, TTN, MYOM1, and AIP, closely related to PKC, which respond immediately to calcium signals in osteoclast activity.
Through the systematic pathogenicity analysis, we identified 24 candidate SNPs and seven candidate InDels in four genes (KRT18, MUC5AC, NBPF9, and PABPC3) and five genes (MSTIL, ASPN, ATN1, PABPC3, and SLAIN1), respectively. KRTs are intermediate filament proteins found in skin and epithelial tissues that play an important role in constructing the cytoskeleton of the cells. KRTs maintain the structure of skin cells and other epithelial tissues by regulating electrolyte transport, post-translational modification, and prevention of degradations; however, disruption of KRTs leads to skin and mucosal diseases. As shown in the results of previous studies, the genetic effect of KRT4 and KRT13 mutations, a member of the keratin family, is mainly implicated in the pathogenesis of white sponge nevus (WSN) disorder, a rare autosomal dominant disorder in oral mucosa [44,45].
KRT18 is also involved in the cytoskeletal signaling pathway, estrogen signaling pathway, and staphylococcus aureus infection pathway [46][47][48][49]. Estrogen generated from the mevalonate pathway contributes to lipid metabolism, producing cholesterol, fatty acids, and steroid hormones, inducing osteoblast proliferation, osteoclast apoptosis, and stimulating matrix synthesis associated with MRONJ [50,51]. Bisphosphonates could influence the production of estrogen, a type of steroid hormone, by inhibiting the mevalonate pathway, and reduced estrogen exacerbates age-related bone loss in women [52]. In addition, a previous study reported that the microbial films of predominant bacteria, including Fusobacterium, Actinomyces, Staphylococcus, Streptococcus, Selenomonas, and Trepenemes, were present in osteonecrosis of MRONJ patients [53][54][55][56]. Disruption of KRT18 could be made to fail to interact with CDH1 (calcium-dependent cell adhesion proteins) and EPCAM (epithelial cell adhesion molecule), which are involved in the proliferation of epithelial cells and defense against mucosal infection, respectively, which affects filament reorganization; therefore, disruption of KRT18 may influence the cytoskeleton of the skin and mucosal, leading to oral mucosa diseases, including MRONJ.
Here, we accomplished the first rare-variant association study for MRONJ in lesions isolated from human jawbone tissue and blood samples. Our genetic expectation demonstrated that the KRT18 and PABPC3 could be potential candidate genes related to MRONJ and developing the MRONJ phenotype. This information could help clinicians to exclude patients who may develop MRONJ during dental procedures and take special preventive measures. Furthermore, it would be an innovative action to prevent drug side effects caused by genotypes through drug sensitivity diagnosis in patients, before the prescription of osteoporosis treatment drugs.
As a limitation of the current research, first, the number of samples in this study was not large enough to detect statistically significant genetic variants. If more robust data is accumulated from patient cohorts with similar clinical environments and conditions, it is believed that it would provide a strong basis from which to derive genetic implications. Second, examining the functional association with germline and somatic variants focusing on MRONJ might provide a crucial message to identify alterations that can be issued in the transcription or translation stage.

Conclusions
We performed whole-exome sequencing to identify MRONJ-associated variants. We identified 954 nonsynonymous SNP variants and 40 InDels accompanied with 691 genes, shared by MRONJ patients in this study. Of these, we found significant genetic variations between MRONJ patients and control sets, and we examined their function through analysis of GO classification and KEGG pathways. The genes bearing the variants were highly related to biological functions of the cytoskeleton, calcium-ion binding, and keratinization. Notably, nine variants located in KRT18 and PABPC3 genes were considered pathogenic variations, influencing immoral protein changes. Our finding suggests that KRT18 and PABPC3 may be implicated in MRONJ pathophysiology, and that they may be used as a therapeutic target in the treatment of osteonecrosis. It is clear that further studies using more cases and controls are needed to confirm the molecular biological association between KRT18 and PABPC3 in the development and occurrence of MRONJ, but we have demonstrated new proven genetic changes to understand the pathogenesis of MRONJ.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/jcm11082145/s1, Figure S1: Outline of experimental design and analysis workflow for this study; Figure S2: Molecular enrichment networks analysis of genes associated with MRONJ; Table S1: Summary of sequencing statistics; Table S2: Number of variants identified from ten MRONJ patients; Table S3: List of SNP candidates; Table S4: List of InDel candidates; Table S5: Enrichment of GO terms and KEGG pathways of MRONJ-associated genes. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Written informed consent has been obtained from the patient(s) to publish this paper.

Data Availability Statement:
The WES data generated in this study are available from the corresponding author upon request.