1. Introduction
The 2022 American Association of Oral and Maxillofacial Surgeons (AAOMS) policy paper defines medication-related osteonecrosis of the jaw (MRONJ) based on the following criteria: no history of radiation therapy to the jaws or evident metastatic disease to the jaws, exposed bone or bone that can be probed through an intraoral or extraoral fistula in the maxillofacial region that has persisted for more than eight weeks, and current or past treatment with antiresorptive and/or antiangiogenic agents [ 1 ]. Antiresorptive agents (e.g., bisphosphonates and denosumab) are widely used to treat osteoporosis, to control hypercalcemia and skeletal pain from bone metastases, and as part of therapy for multiple myeloma, either as a standalone treatment or in combination with antiangiogenic or immunomodulatory agents [ 1 ]. At present, due to the influence of pharmaceutical agents targeting the osteoclastogenesis process, patients suffering from a variety of pathologies—including, but not limited to, giant cell tumors, rheumatoid arthritis, Paget’s disease of bone, fibrous dysplasia, and osteogenesis imperfecta—require these pharmaceutical agents to manage their respective conditions [ 2 ]. As a result, the number of individuals susceptible to the development of MRONJ has increased significantly. The impact of MRONJ on quality of life is substantial [ 3 ], particularly in advanced cases characterized by severe pain and impaired oral intake. In patients with oncologic diseases, the condition has the potential to be life-threatening. Consequently, elucidating the pathophysiology of MRONJ is imperative to enhance prevention and treatment strategies. Even though research on MRONJ has been ongoing since its initial report in 2003 [ 4 ], the precise etiology of the condition remains uncertain, and a universally accepted treatment plan is yet to be established. The pathogenesis of MRONJ is complex, involving an interacting array of genetic, metabolic, and immunological factors [ 5 ]. A multitude of hypotheses have been proposed to explain the development of MRONJ. However, the predominant hypothesis suggests that it is caused by excessive suppression of bone metabolism and remodeling triggered by anti-resorptive medications [ 6 ]. The anti-resorptive agent in question exerts its effect by inhibiting both the process of osteoclastic bone resorption and that of osteoblastic bone formation. This effect arises from an interruption in the coupling of osteoblasts and osteoclasts. As posited in the research by Kim, J.-M. et al. [ 7 ], anti-resorptive agents have been shown to reduce the release of osteogenic clastic factors by inhibiting osteoclast differentiation. In addition to this, the release of osteogenic growth factors embedded within the bone matrix is also reduced by anti-resorptive agents as a result of decreased bone resorption. Researchers also deliberate on hypotheses that encompass oral infection [ 8 ] and inflammation [ 9 ]. The antiangiogenic impact and soft tissue toxicity of bisphosphonates are also assessed [ 10 ]. The discussion concludes with a focus on the inhibition of lymphangiogenesis [ 11 ]. Osteonecrosis is a condition characterized by the death of local bone and marrow cells due to continuous stimulation by detrimental factors, resulting in aberrant osteoimmune function [ 12 ]. This phenomenon leads to the occurrence of a chronic inflammatory microenvironment that inhibits bone regeneration and repair. In a plethora of bone tissue diseases, innate immune cells and adaptive immune cells interact with bone cells. The effects of this interaction on bone metabolic homeostasis have attracted increasing attention, resulting in the establishment of a new field: osteoimmunology. A thorough examination has been carried out to identify the regulatory mechanisms that control the inflammatory response, revealing that immune cells are the principal components involved. As a result, osteoimmune dysfunction has attracted considerable interest from the scientific community as a potential causative factor of osteonecrosis [ 13 ]. The onset of MRONJ and the activity of particular genes may be related, according to earlier research. Scientists analyzed gene expression data from individuals with multiple myeloma derived from microarray experiments and identified seven hub genes ( ACTB , FN1 , GAPDH , JUN , PTPRC , STAT3 , and TNF ) linked to bisphosphonate-induced osteonecrosis of the jaw [ 14 ]. However, dysregulation of GAPDH expression in MRONJ in patients with multiple myeloma should be approached with caution, as GAPDH is implicated in multiple cancer-related biological processes, having recently emerged as a primary focus in cancer research [ 15 , 16 ]. The preceding analyses resulted in the identification of three salivary proteins ( AACT , HBD , and MMP9 ) [ 17 ]; 17 key genes; and hsa-mir-16-1, hsa-mir-21, hsa-mir-23a, hsa-mir-145, hsa-mir-186, hsa-mir-221, and hsa-mir-424 [ 18 ] associated with MRONJ.
Despite the noteworthy advancements in elucidating the mechanisms underlying MRONJ, challenges persist with regard to the early diagnosis and treatment of this condition due to its heterogeneous radiological manifestations and unpredictable course [ 19 ].
Traditional diagnostic methods often prove ineffective, limiting their ability to provide timely diagnosis, treatment, and prognosis for patients. Genetic testing can serve as an important tool in the diagnosis and treatment of MRONJ, increasing the accuracy and effectiveness of both diagnostic and therapeutic approaches [ 20 ]. The identification of genetic markers associated with MRONJ will enable earlier and more accurate diagnosis, facilitating timely treatment and improving patient care. Genetic screening tests have therefore become a key element of therapeutic practice, increasing the accuracy of diagnosis and facilitating personalized treatment strategies. Identifying disease-associated genes via genetic testing, coupled with the verification of relevant biomarkers, would enhance the speed and efficiency of MRONJ diagnosis and treatment, thereby improving both short-term clinical outcomes and long-term disease management.
This study will investigate genes, biomarkers, and molecular processes related to MRONJ through comprehensive bioinformatics analysis to improve MRONJ diagnosis.
2. Results
2.1. Screening Analysis of the DEGs
To investigate the genes associated with the onset of MRONJ, a retrospective analysis was conducted using PBMC expression profiles from the GSE7116 dataset. This dataset included samples from 11 MM_ONJ patients, 10 age-matched MM controls, and 5 Healthy controls.
Raw microarray data were normalized and preprocessed for downstream analysis ( Figure 1 A). Differential gene expression analysis identified 22,173 transcripts per comparison. Using a significance threshold of p < 0.05 and |log 2 fold change| > 1, a total of 124 DEGs were identified in the Healthy vs. MM_ONJ comparison, and 426 DEGs in the MM vs. MM_ONJ comparison ( Supplementary Materials: S1A, S1B ). A visual representation of these findings was achieved using a heat map, which highlights the top DEGs ( Figure 1 B). Furthermore, the implementation of volcano plots ( Figure 1 C,D) elucidated genes that exhibited significant upregulation, in contrast to those that demonstrated significant downregulation.
2.2. Gene Set Enrichment Analysis
GSEA was performed for both distinct comparisons, Healthy vs. MM_ONJ and MM vs. MM_ONJ, to gain more profound insight into the biological pathways associated with the expression of genes related to MRONJ. In the Healthy vs. MM_ONJ comparison, the DEGs were significantly enriched in those pathways that were related to protein synthesis and ribosomal function ( Figure 2 ). The top five upregulated pathways (based on positive Normalized Enrichment Score, NES, and p < 0.001) included Eukaryotic Translation Elongation, Eukaryotic Translation Initiation, Ribosome, Response of EIF2AK4 ( GCN2 ) to Amino Acid Deficiency, and rRNA Processing.
On the other hand, the top five downregulated pathways in the Healthy vs. MM_ONJ group ( Figure 2 ) were primarily associated with innate immune defense, in particular neutrophil and antimicrobial activities. These included Antimicrobial Peptides, Nicotinate and Nicotinamide Metabolism, Egr2 and SOX10 -Mediated Initiation of Schwann Cell Myelination, Neutrophil Degranulation, and Defensins. Our findings indicate that for multiple myeloma patients with MRONJ disease, a significant suppression of neutrophil activity and host defense mechanisms can be observed.
For the MM vs. MM_ONJ comparison ( Figure 3 ), both up- and downregulated pathway signatures were observed. The top five significantly upregulated pathways included T Cell Receptor (TCR) Signaling, FcεRI-Mediated NF-κB Activation, Leishmania Infection, Activation of mRNA upon Binding of the Cap-Binding Complex and eIFs, and Subsequent Binding to 43S, SRP-Dependent Cotranslational Protein Targeting to Membrane. At the same time, the five pathways with the most significant downregulation were: Keratinization, Voltage-Gated Potassium Channels, Olfactory Signaling Pathway, Receptor-Type Tyrosine–Protein Phosphatases, and Formation of the Cornified Envelope. All pathways mentioned above exhibited statistically significant negative NES and adj. p < 0.05, with many reaching thresholds of adj. p < 0.001, thereby reinforcing their potential biological relevance to MRONJ pathogenesis.
2.3. Weighted Gene Co-Expression Network Analysis
WGCNA is a systems biology method that operates under the assumption that genes exhibiting similar expression profiles may often be functionally linked, co-regulated, or involved in the same biological pathways. This approach facilitates the establishment of a gene similarity network, which can subsequently be utilized to identify clusters comprising genes that exhibit a high degree of correlation. Subsequently, these genes are organized into modules through a hierarchical clustering process. Module eigengenes (MEs), defined as the first principal component of each module, were computed and correlated with binary clinical traits (Healthy vs. MM_ONJ; MM vs. MM_ONJ). A sample dendrogram was generated for each comparison to visualize relationships among samples ( Figure 4 A and Figure 5 A). To determine the appropriate soft-thresholding power (β) for scale-free network topology, values of β = 15 (for Healthy vs. MM_ONJ) and β = 13 (for MM vs. MM_ONJ) were selected based on scale-free topology fit indices ( Figure 4 B and Figure 5 B). For the Healthy vs. MM_ONJ comparison, the soft-thresholding power β = 15 was selected because it provided a stable, high scale-free topology fit (R 2 plateau), reduced sensitivity to noise, and more biologically interpretable co-expression modules while simultaneously maintaining an acceptable mean connectivity. For the MM vs. MM_ONJ dataset, no soft-thresholding power reached the conventional R 2 ≈ 0.85, likely due to limited sample size and weaker global co-expression structure. In such situations, the power that best approaches scale-free topology while maintaining stable connectivity was chosen. Therefore, β = 13 was selected, which showed the highest attainable R 2 , the most stable connectivity pattern, and minimal noise—representing the optimal practical compromise when an ideal scale-free fit cannot be achieved. Network construction and module detection were carried out using the following parameters: minModuleSize = 20 and mergeCutHeight = 0.15. Modules with eigengene correlations greater than 0.75 were merged, resulting in 39 and 31 merged modules for the two comparisons, respectively ( Figure 4 C and Figure 5 C). Both raw and merged module dendrograms were visualized ( Figure 4 D and Figure 5 D), with modules color-coded for ease of distinction. According to ME-trait heatmap analysis, four modules showed statistically significant associations with the onset of MRONJ in the Healthy vs. MM_ONJ comparison: MEantiquewhite4, MEpalevioletred3, MEfirebrick4, and MElightcyan1 ( Figure 4 E), comprising a total of 1261 genes. In the MM vs. MM_ONJ comparison, five modules were significantly associated with MRONJ: MElightcyan, MEcyan, MEdarkmagenta, MElightsteelblue1, and MEskyblue2 ( Figure 5 E), comprising 2512 genes in total. Modules exhibiting the strongest correlations with disease status were selected for downstream analysis.
2.4. Integrated CytoHubba-Identified PPI Hubs and LASSO-Selected Predictive Genes
In the Healthy vs. MM_ONJ comparison, a total of 40 hub module DEGs were identified by intersecting the 1261 hub module genes from four WGCNA-derived modules with the 124 DEGs from the previous analysis (significance threshold: p < 0.05 and |log 2 fold change| > 1) ( Figure 6 A). In the MM vs. MM_ONJ comparison, 116 DEGs were identified by intersecting the 2512 hub genes from five WGCNA-derived modules with the 426 DEGs from the earlier screening ( Figure 6 B). The gene selection process was performed independently using two complementary approaches: LASSO regression and CytoHubba-based PPI network analysis, applied to both sets of hub genes under investigation.
In LASSO regression, λ min corresponds to the regularization parameter that minimizes CV error, while λ 1Se is the largest value within one standard error of this minimum, favoring simpler models with fewer predictors. For the Healthy vs. MM_ONJ comparison, LASSO regression with 10-fold CV was performed ( Figure 6 C,D). The optimal λ values were λ min ≈ e −4.1 and λ 1Se ≈ e −3.1 . At λ min , six non-zero coefficients were retained. When using λ 1Se, a more parsimonious model with five non-zero coefficients ( FCMR , ZAP70 , HBG2 , ADM , NAMPT ) was obtained. Subsequently, the parameter λ 1Se was selected for further analyses to optimize the balance between model simplicity and its predictive accuracy. For the MM vs. MM_ONJ comparison, LASSO regression with 10-fold CV was also performed ( Figure 6 E,F). The optimal λ values were λ min ≈ e −4.8 and λ 1Se ≈ e −3.1 . At the minimum value of the parameter λ, eight non-zero coefficients were retained. The parameter λ 1Se yielded a more parsimonious model with five non-zero coefficients ( IGKV4-1 , PPP6R2 , XPO7 , BTG2 , PDE4B ). As was the case in the preceding comparison, the value of the parameter λ 1Se was selected in order to prioritize the simplicity of the model while ensuring the maintenance of a high degree of predictive performance. The network topology of the 40 and 116 DEGs identified in both comparisons was analyzed using the CytoHubba plugin in Cytoscape. Highly connected hub genes were extracted from the main complex DEG network using CytoHubba local-based methods—Maximum Neighborhood Component (MNC), Maximal Clique Centrality (MCC), Density of Maximum Neighborhood Component (DMNC), and Degree (Deg)—as well as global-based methods, including Betweenness Centrality (BC), BottleNeck (BN), Closeness Centrality (Clo), EcCentricity (EC), Edge Percolated Component (EPC), Stress (Str), and Radiality (Rad). From this analysis, 16 genes were identified that appeared in the intersection of at least three ranking methods: BCL2L1 , BCL10 , CCL3 , CD44 , CREB1 , CXCL1 , CXCL3 , EIF4A2 , ETS1 , EZR , IFIH1 , IL1B , JAK1 , PPP1R15A , TNFAIP3 , and ZAP70 ( Figure 6 G).
2.5. Gene Ontology and Pathway Enrichment
GO enrichment analysis was performed on the integrated set of 25 DEGs obtained from the intersection of CytoHubba and LASSO results ( BCL2L1 , ADM , BCL10 , BTG2 , CCL3 , CD44 , CREB1 , CXCL1 , CXCL3 , EIF4A2 , ETS1 , EZR , FCMR , HBG2 , IFIH1 , IGKV4-1 , IL1B , JAK1 , NAMPT , PDE4B , PPP1R15A , PPP6R2 , TNFAIP3 , XPO7 , ZAP70 ) using the BiNGO plugin in Cytoscape. The Benjamini–Hochberg false discovery rate method, with a significance level of p < 0.05, was employed to correct for the multiple testing. The three distinct domains: biological process (GO_BP), molecular function (GO_MF), and cellular component (GO_CC) were implemented. The selection of these domains was made to provide illustrative examples of the functional roles and regulatory context of the identified gene sets ( Figure 7 A). The top three enriched GO_BP terms, ranked by adjusted p -value, were as follows: GO:0070486, leukocyte aggregation; GO:0022409, positive regulation of cell–cell adhesion; and GO:0030593, neutrophil chemotaxis. The following key terms were identified for GO_MF: GO:004237, chemokine receptor binding; GO:0005126, cytokine receptor binding; and GO:0008009, chemokine activity. In GO_CC, the top three terms were classified as such: GO:0001772, immunological synapse; GO:0031234, extrinsic component of cytoplasmic side of plasma membrane; and GO:0045121, membrane raft. Pathway enrichment analysis of the integrated set of DEGs, performed using both KEGG and Reactome databases, revealed the biological relevance of the hub genes within the following signaling and metabolic pathways: KEGG: hsa04064 NF-kappa B signaling pathway, hsa04621 NOD-like receptor signaling pathway, hsa04668 TNF signaling pathway; Reactome: R-HSA-6783783 Interleukin-10 signaling, R-HSA-449147 Signaling by interleukins, R-HSA-380108 Chemokine receptors bind chemokines. All pathways met the significance threshold of p < 0.05 ( Figure 7 B,C). In parallel, a combined KEGG and Reactome pathway enrichment analysis was performed using the ClueGO plugin in Cytoscape to further elucidate the biological relevance of the hub genes. The Clue-GO network, including the hub genes, is shown in Figure 7 D. The most significant group of pathways was that with Interleukin-10 pathway hub.
2.6. Violin Plot and ROC Curve Analysis
RMA-normalized gene expression data were aggregated by gene symbol, and both the expression matrix and metadata were reshaped into long format for visualization. The expression patterns of 25 genes of interest across three experimental groups Healthy, MM, and MM_ONJ ( BCL2L1 , ADM , BCL10 , BTG2 , CCL3 , CD44 , CREB1 , CXCL1 , CXCL3 , EIF4A2 , ETS1 , EZR , FCMR , HBG2 , IFIH1 , IGKV4-1 , IL1B , JAK1 , NAMPT , PDE4B , PPP1R15A , PPP6R2 , TNFAIP3 , XPO7 , and ZAP70 ) were then visualized ( Figure 8 A). The assessment of differences between groups was conducted with the aid of the Wilcoxon rank-sum test, followed by false discovery rate (FDR) correction to account for multiple comparisons. This analysis revealed group-specific expression trends, with several genes showing statistically significant differences in expression between comparison groups Healthy vs. MM_ONJ and MM vs. MM_ONJ ( Table 1 ). In addition, to assess the diagnostic potential of individual genes, ROC curve analysis was performed. Genes with an AUC > 0.9 and an adj. p < 0.003 ( EIF4A2 , ETS1 , FCMR , IGKV4-1 , JAK1 , PDE4B , XPO7 ) were considered highly discriminative and selected for further investigation ( Figure 8 B).
2.7. Transcriptional Signature and Predictive Modeling of MRONJ-Associated Gene Expression
The generation of heatmaps enabled further evaluation of the expression patterns of candidate genes associated with osteonecrosis of the jaw. The comparisons used were Healthy vs. MM_ONJ and MM vs. MM_ONJ, with the use of z-score normalized expression values. In both comparisons, seven selected genes ( PDE4B , ETS1 , FCMR , JAK1 , EIF4A2 , IGKV4-1 , and XPO7 ) demonstrated consistent upregulation in the MM_ONJ group relative to controls ( Figure 9 ). In the Healthy vs. MM_ONJ heatmap, clear separation between groups was observed. The higher expression of all seven genes was demonstrated in MM_ONJ samples ( Figure 9 A). This trend also remained evident in the MM vs. MM_ONJ comparison, although with slightly reduced contrast ( Figure 9 B). The hierarchical grouping of samples demonstrated expression profiles that were specific to individual groups, thereby highlighting the discriminatory power of the genes that were selected. The most significant increase in expression was observed in the PDE4B , FCMR , XPO7 , and EIF4A2 genes in the MRONJ samples. These findings support the hypothesis that these genes may serve as potential biomarkers for MRONJ. Furthermore, they imply the presence of a distinct transcriptional signature associated with MRONJ development, separate from the transcriptional changes generally associated with MM.
The predictive accuracy of gene expression signatures was evaluated by minimizing bias resulting from small sample sizes and imbalanced groups. For this purpose, Firth’s logistic regression with reduced bias within nested CV was employed. The approach entailed the implementation of an outer loop, the function of which was to evaluate generalization performance. An inner loop was also incorporated, the purpose of which was to select features based on one-dimensional t-tests. In each iteration of the outer fold, the five genes exhibiting the lowest p -values were selected for model training.
Two binary classification tasks were evaluated: MM vs. MM_ONJ and Healthy vs. MM_ONJ ( Figure 10 ). In a comparative analysis involving MM and MM_ONJ groups, an AUC of 0.982 was achieved ( Figure 10 A). This result represents an almost flawless distinction between these two different groups. A comparable result was observed in the comparison of the Healthy and MM_ONJ groups, where an AUC value of 0.945 was achieved ( Figure 10 B), which also indicates exceptional classification performance. These high AUC values suggest that MRONJ cases are molecularly distinct from both MM patients without MRONJ and from healthy controls. These findings collectively validate the strong discriminatory ability of the chosen gene signatures in detecting MRONJ.
We also evaluated the discriminatory ability of fixed LASSO-selected gene sets for Healthy vs. MM_ONJ and MM vs. MM_ONJ comparisons using PCA followed by ridge-penalized logistic regression within a safe stratified 0.632+ bootstrap framework. For each comparison, the dataset was restricted to the predefined four-gene panel, transformed to a samples × genes matrix, and analyzed across 200 bootstrap iterations with 200 matched stratified permutation tests to generate null AUC distributions. Model safety checks ensured adequate class representation in all resampling stages and folds. Healthy vs. MM_ONJ: Using the LASSO-selected panel ( FCMR , IGKV4-1 , JAK1 , EIF4A2 ), the bootstrap-derived OOB mean AUC was 0.990, with a median of 1.000 and a 95% CI of [1.000, 1.000]. The 0.632+ mean AUC was 1.000. The empirical p -value was 0.0405, and the effect size (Cohen’s d) was 3.99, indicating an enormous separation between bootstrap and permutation distributions. The permutation AUC distribution was centered near the chance level (~0.5) with no overlap with the bootstrap distribution ( Figure 11 A,B). MM vs. MM_ONJ: For the LASSO-selected panel ( PDE4B , EIF4A2 , JAK1 , XPO7 ), the bootstrap-derived OOB mean AUC was 0.979, with a median of 1.000 and a 95% CI of [0.875, 1.000]. The 0.632+ mean AUC was 0.991. The empirical p -value was 0.0152, and Cohen’s d was 4.78, again indicating an extremely large separation between distributions. While discrimination remained very strong, the bootstrap distribution showed slightly greater spread compared to the Healthy vs. MM_ONJ analysis, suggesting higher biological heterogeneity in this comparison ( Figure 11 C,D). Both gene panels achieved near-perfect classification performance under rigorous bootstrap evaluation. The complete or near-complete separation between bootstrap and permutation distributions, together with large effect sizes, confirms that these results are highly unlikely under the null hypothesis and robust to resampling variability. Across both comparisons, the complete separation between bootstrap and permutation distributions confirms that the models capture biologically meaningful differences in gene expression rather than artifacts of sampling or overfitting.
3. Discussion
In the present study, an integrative transcriptomic and machine-learning pipeline was implemented for the objective of identifying and validating systemic biomarkers of MRONJ in patients diagnosed with multiple myeloma. Through the combination of differential gene expression analysis, WGCNA, LASSO regression, bootstrap and permutation validation, and CytoHubba network ranking, we delineated molecular signatures implicating immune dysregulation, inflammatory signaling, and, as a hypothesis extension, maladaptive stress responses that may exacerbate bone necrosis.
3.1. Immune and Inflammatory Pathways in MRONJ
Our DEG and enrichment analyses revealed that the Healthy vs. MM_ONJ comparison exhibited consistent upregulation of translational and ribosomal machinery pathways (e.g., eukaryotic translation elongation, initiation, and rRNA processing) alongside suppression of neutrophil degranulation, antimicrobial peptide signaling, and defensin pathways. This pattern suggests an enhanced translational drive accompanied by weakened innate immune and antimicrobial responses, potentially compromising mucosal defense and jaw tissue resilience [ 21 , 22 , 23 , 24 ]. In contrast, analyses of MM vs. MM_ONJ comparison revealed an augmented state of adaptive immune activation, as indicated by increased regulation of T-cell signaling (TCR), FcεRI-mediated NF-κB activation, and mRNA cap-binding/eIF pathways, along with diminished regulation of keratinization, potassium channels, and epithelial barrier processes [ 25 , 26 ]. It is hypothesized that the aforementioned complementary characteristics of MRONJ imply a disturbance in immune-bone homeostasis, wherein uncontrolled immune stimulation is inadequately counterbalanced by protective innate and barrier mechanisms. Intersecting DEGs with WGCNA modules refined these findings to 40 hub module DEGs (Healthy vs. MM_ONJ) and 116 hub module DEGs (MM vs. MM_ONJ). Focusing on genes embedded within coherent co-expression networks rather than isolated DEGs improves biological interpretability and highlights modules likely reflecting systemic immunological perturbations associated with MRONJ susceptibility [ 27 , 28 , 29 ].
3.2. LASSO Selection and Model Validation
LASSO regression with ten-fold cross-validation distilled compact, interpretable gene panels. In Healthy vs. MM_ONJ, the λ 1 ₛₑ model retained five predictors— FCMR , ZAP70 , ADM , NAMPT , and HBG2 —whereas in MM vs. MM_ONJ, the selected set comprised IGKV4-1 , PPP6R2 , XPO7 , BTG2 , and PDE4B . As HBG2 is a globin gene that is likely indicative of pre-analytical contamination, it was excluded from further analyses. Notably, the classifier exhibited consistent performance in the absence of HBG2 , thereby confirming the biological validity of the results, which extend beyond the confines of technical artifacts. Complementary validation strategies reinforced these results. Nested cross-validation using Firth logistic regression achieved AUC = 0.945 for Healthy vs. MM_ONJ and AUC = 0.982 for MM vs. MM_ONJ. PCA + ridge logistic regression, embedded in a safe stratified 0.632+ bootstrap, yielded median bootstrap AUCs of 1.00 in both comparisons (Cohen’s d ≈ 3.99–4.78). Permutation null distributions were completely separated from bootstrap AUC distributions, indicating that discrimination was not driven by random variation. These findings collectively affirm the reproducibility and generalizability of the gene panels. While the model achieved AUC > 0.9 in internal validation, external cohort validation and tissue-level assays will be required to confirm the diagnostic utility of the proposed biomarkers. However, it is crucial to acknowledge that the bootstrap and permutation methods employed to assess the accuracy of the model’s performance are not inherently immune to bias. An AUC of 1.00 may reflect limited cohort diversity and should be interpreted cautiously.
3.3. Hub Network and Biological Interpretation
CytoHubba network analysis identified sixteen highly connected hub genes in both comparisons, including JAK1 , IL1B , TNFAIP3 , ETS1 , CXCL1 , and CXCL3 . Many of these genes have well-described roles in cytokine signaling, NF-κB regulation, and immune-cell activation, which are central pathways in osteoimmunology. JAK1 is a key component of several cytokine receptor systems (for example, IL-6 and γc cytokines) [ 30 , 31 ] and can influence osteoclast activity by modulating RANKL expression [ 32 ]. TNFAIP3 ( A20 ) acts as a negative regulator of NF-κB and helps maintain limits on inflammatory signaling and osteoclast differentiation [ 33 ]. IL1B , CXCL1 , and CXCL3 promote neutrophil recruitment and inflammatory cell movement [ 34 , 35 ]; these processes support mucosal healing but when prolonged may contribute to tissue injury. Persistent chemokine signaling in the oral mucosa is consistent with the non-healing ulcers seen in MRONJ. ETS1 is a transcription factor involved in immune-cell differentiation and inflammatory gene expression. It can also affect bone remodeling by increasing RANKL and MMP transcription in fibroblasts [ 36 ], linking inflammation with tissue breakdown and impaired repair [ 37 ]. Although direct evidence in MRONJ is limited, these mechanisms support their relevance to immune-bone interactions.
Overall, the network findings point to sustained cytokine and chemokine dysregulation in MRONJ, where unresolved inflammation may drive local tissue damage and impaired bone healing.
3.4. Hypothesis: Stress Responses as Co-Modifiers
Although immune dysregulation in MRONJ is well recognized [ 38 , 39 ], our enrichment analysis also pointed to changes in translation, ribosomal biogenesis, and eIF signaling. These findings are consistent with activation of the integrated stress response (ISR). We suggest that impaired stress-handling pathways—combining altered translational control, oxidative stress, and the unfolded protein response—may worsen immune-mediated tissue injury. There is mechanistic support for this idea. ISR signaling ( PERK → eIF2α → ATF4 ) helps cells adapt to metabolic and oxidative stress. In osteoblasts, ATF4 is required for normal matrix production and stress tolerance, and its disruption can limit osteogenesis and hinder bone repair [ 40 , 41 ]. Oxidative stress also promotes osteoclast differentiation and bone resorption [ 42 ] and increased reactive oxygen species have been shown to intensify bisphosphonate-associated necrosis in animal studies [ 43 ].
Considering these findings, we hypothesize that in patients with MM who are treated with antiresorptive agents, persistent immune activation and the reduced capacity to manage cellular stress, which are concomitant with the treatment, increase the susceptibility to jaw osteonecrosis. To validate this model, future research should assess ISR-related markers (e.g., phosphorylated eIF2α and ATF4 -dependent genes), oxidative stress markers, and immune cell phenotypes in patient PBMCs and affected bone tissue.
3.5. Limitations and Future Directions
The present study has limitations. The modest sample size and limited clinical diversity of the dataset analyzed constrain the generalizability of the results. Although resampling and bias-reduction techniques (nested CV, Firth logistic regression, and 0.632+ bootstrap) were implemented to mitigate overfitting, external validation in larger, prospectively collected, and demographically heterogeneous cohorts remains essential to confirm reproducibility. PBMC transcriptomes serve as systemic surrogates for the jaw microenvironment. Thus, tissue-level validation of the identified biomarkers is essential to establish causality. Furthermore, the proposed stress-response hypothesis remains inferential and requires mechanistic experimentation.
The findings should therefore be viewed as exploratory. Future studies should include larger, prospectively collected cohorts and incorporate single-cell or spatial transcriptomics of jaw lesions to confirm cell-type-specific expression of the hub genes. We plan to validate our seven-gene MRONJ signature in a prospective multi-center study using independent patient cohorts recruited from collaborating institutions. Standardized protocols for sample collection and processing will be implemented to minimize inter-site variability. In parallel, we will perform targeted transcript quantification using the NanoString nCounter ® platform (either a Panel Plus–enhanced immune panel or a custom CodeSet) to confirm differential expression with a highly reproducible, amplification-free technology. These combined efforts will provide robust external validation and assess the clinical utility of the proposed biomarker panel in broader patient populations.
Preclinical experimentation involving the modulation of JAK—STAT, NF-κB, or ISR signaling in MRONJ animal models may facilitate the elucidation of causal pathways. Also, the integration of proteomic and phospho-signaling datasets has the potential to offer insights into post-transcriptional regulation and signaling activity.
Although PBMC transcriptomes reflect systemic immune and stress-related dysregulation, they may not fully capture the local molecular landscape of necrotic bone tissue. However, several hub genes identified here ( JAK1 , TNFAIP3 , ETS1 , and CXCL family members) have been independently reported as upregulated in MRONJ jawbone lesions or gingival tissues (e.g., [ 17 , 21 ]), suggesting partial correspondence between systemic and local expression profiles. In addition to tissue-level studies, salivary extracellular vesicles (EVs) represent an emerging, non-invasive source of molecular information, as their proteomic and transcriptomic cargo can mirror local inflammatory and osteolytic processes occurring in the oral cavity. Integrating salivary EV profiling with PBMC-based signatures may therefore help bridge systemic and site-specific pathology. Future work employing single-cell or spatial transcriptomics of jaw lesions, combined with EV-based molecular profiling, will be essential to establish cell-type-specific expression of these hub genes and clarify systemic-local signaling relationships in MRONJ.
3.6. Clinical Implications
The seven-gene signature ( PDE4B , ETS1 , FCMR , JAK1 , EIF4A2 , IGKV4-1 , XPO7 ) offers a practical foundation for developing a multiplex PCR or digital droplet PCR (ddPCR) assay. Future translational steps will involve optimizing assay sensitivity and specificity, defining diagnostic thresholds based on ROC-derived cutoffs, and evaluating cost efficiency relative to current imaging-based screening. Integration with clinical predictors (e.g., duration of bisphosphonate therapy, comorbidities) may enable composite risk scores suitable for clinical decision support. Ultimately, panel refinement and validation could precede prospective evaluation in a biomarker-guided clinical trial.
The hub genes and molecular patterns identified here may have clinical relevance. A peripheral blood transcriptomic panel that includes JAK1 , TNFAIP3 , CXCL1/3 , ETS1 , and FCMR could potentially serve as a minimally invasive tool for assessing MRONJ risk before starting antiresorptive or antiangiogenic treatment. The implicated pathways—such as NF-κB, JAK–STAT, and chemokine signaling—also suggest opportunities for therapeutic modulation, for example, through cytokine inhibitors or JAK inhibitors, to limit inflammatory bone loss.
If the stress-response hypothesis is confirmed, treatments aimed at improving proteostasis or reducing oxidative stress (such as antioxidants or ISR-modulating agents) may complement current preventive approaches. Ultimately, biomarker-guided patient selection combined with pathway-targeted intervention could lower MRONJ risk and improve treatment safety in oncology and osteoporosis care.
4. Materials and Methods
4.1. Data Resource Collection
Gene expression microarray data from patients with osteonecrosis of the jaw [ 44 ] were obtained from the Gene Expression Omnibus (GEO) database (accession GSE7116 ) (National Institutes of Health, Bethesda, MD, USA; http://www.ncbi.nlm.nih.gov/geo , accessed on 10 June 2025) using the R programming language. The dataset was obtained utilizing the keywords ‘osteonecrosis of the jaw’, ‘ Homo sapiens ’, and ‘ GPL570 platform’. It comprises transcriptional profiles of peripheral blood mononuclear cells (PBMCs) from 11 multiple myeloma patients with MRONJ (MM_ONJ: GSM170982 .CEL, GSM170984 .CEL, GSM170985 .CEL, GSM170987 .CEL, GSM170988 .CEL, GSM170989 .CEL, GSM170990 .CEL, GSM170991 .CEL, GSM170993 .CEL, GSM170994 .CEL, GSM170995 .CEL), compared with 10 age-matched multiple myeloma controls (MM: GSM170996 .CEL, GSM170997 .CEL, GSM170998 .CEL, GSM170999 .CEL, GSM171000 .CEL, GSM171001 .CEL, GSM171002 .CEL, GSM171003 .CEL, GSM171004 .CEL, GSM171005 .CEL) and five healthy volunteers (Healthy: GSM171006 .CEL, GSM171007 .CEL, GSM171008 .CEL, GSM171009 .CEL, GSM171010 .CEL). Patients (aged 57–81 years) had received pamidronate (n = 3), zoledronate (n = 4), or both (n = 4) for an average of 38.7 months. The GPL570 platform (Affymetrix GeneChip Human Genome U133 Plus 2.0 Array, Affy-metrix, Santa Clara, CA, USA) was used.
4.2. Data Normalization and Differential Expression Analysis
The raw microarray data were processed with the Affy package (v 1.68.0; Bio-conductor) in R (v 4.4.2; CRAN). The dataset was obtained from the GEO database through the GEOquery package (v 2.74.0). To ensure the acquisition of high-quality, standardized expression data, a series of procedures were implemented. Background correction and normalization were first performed using the Robust Multi-array Average (RMA) algorithm, followed by quantile normalization. The analysis of differential expression was performed utilizing the limma package (v 3.62.2). The lmFit and eBayes functions were utilized to model expression data and calculate moderated statistics for differential gene expression across groups. A comparative analysis was performed involving patients diagnosed with multiple myeloma and MRONJ (MM_ONJ) alongside two control groups: patients with multiple myeloma without MRONJ (MM) and healthy individuals (Healthy). A cutoff criterion of absolute log2-fold change (|log 2 FC|) > 1 and p < 0.05 was utilized to identify significantly differentially expressed genes (DEGs). The results were visualized using a heatmap that displayed the top 30 DEGs according to their (by |log 2 FC|) values, alongside volcano plots illustrating both upregulated and downregulated genes. All visualizations were generated using the ggplot2 (v 4.0.0) and pheatmap (v 1.0.13) packages in the R programming environment.
4.3. Gene Set Enrichment Analysis
A Gene Set Enrichment Analysis (GSEA) was performed on pre-ranked gene lists using the clusterProfiler package (v 4.14.6) in R or the GSEA (v 4.3.3) desktop application (Broad Institute) using the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome pathway curated databases.
To evaluate the statistical significance of the observed results, enrichment scores were computed utilizing a Kolmogorov–Smirnov-like statistic with 100,000 permutations. This study identified pathways with an adjusted p -value ( adj. p < 0.05) as significantly enriched. The ten pathways identified as significantly overexpressed, and the ten pathways identified as significantly underexpressed were subsequently ranked based on their Normalized Enrichment Score (NES). To enhance understanding, these rankings were later illustrated through bar plots, created with the R statistical software (v 4.4.2).
4.4. Weighted Gene Co-Expression Network Analysis
A Weighted Gene Co-expression Network Analysis (WGCNA) [ 45 ] was carried out to identify modules of co-expressed genes that are associated with clinical traits (WGCNA, v1.73). The raw microarray data underwent a series of standard procedures, including normalization using the Robust Multi-array Average (RMA) algorithm and log 2 -transformation of expression values. The low-variance genes were filtered out by keeping the top 75% of the most variable probes. The quality of the samples was then assessed using hierarchical clustering, a method that detects and excludes outliers.
The selection of a soft-thresholding power (β) was made using the scale-free topology criterion (R 2 ≥ 0.85), with the fit indices being computed over a range of powers. An unsigned network was constructed using the blockwiseModules function. To assess connectivity among genes, the adjacency and topological overlap matrices (TOM) were calculated. Initially, modules were defined by dynamic tree cutting. Subsequently, they were merged if their eigengene dissimilarity was less than 0.25 (i.e., eigengene correlation > 0.75). The mergeCloseModules function was used to carry out these merges. Module eigengenes (MEs), which correspond to the initial principal component of each module, were correlated with binary clinical traits (Healthy vs. MM_ONJ, MM vs. MM_ONJ) employing Pearson’s correlation. Subsequently, statistically significant module-trait associations ( p < 0.05) were represented visually via labeled heatmaps. Probe identifiers were mapped to gene symbols using the hgu133plus2.db annotation package.
4.5. Integrated Hub Gene Identification: LASSO and CytoHubba
To determine dominant hub genes, two distinct approaches were used: Least Absolute Shrinkage and Selection Operator (LASSO) regression and CytoHubba-based protein–protein interaction (PPI) network analysis. These analyses were conducted separately for the Healthy vs. MM_ONJ and MM vs. MM_ONJ comparisons. LASSO was utilized for the selection of predictive genes from expression data, whereas Cyto-Hubba identified topological hub genes based on network connectivity.
4.5.1. LASSO Regression
The gene expression data were divided into subsets and normalized prior to modeling. For each comparison, a penalized logistic regression model was fitted using the glmnet package (v 4.1.10) in R, and cross-validation (CV) via cv.glmnet was used to determine the optimal penalty parameter (lambda.min). Genes with non-zero coefficients at lambda.min were identified as LASSO-derived hub genes. When only a single predictor remained, standard logistic regression was applied instead.
4.5.2. CytoHubba
Protein–protein interaction networks were constructed and analyzed using Cytoscape (v 3.10.2) [ 46 ]. The STRING database was used as the source for interaction data, and the StringApp (v 2.1.1) [ 47 ] extension facilitated data import and visualization. The CytoHubba (v 0.1) plugin [ 48 ] was employed to identify hub genes using all available scoring algorithms (e.g., MCC, Degree, EPC). To this aim, an UpSet plot was generated to identify genes that appeared in at least three of the ranking methods, thus highlighting consistently top-ranked hub candidates.
4.6. Gene Ontology and Pathway Enrichment
Gene Ontology (GO) enrichment analysis was performed using the BiNGO (v 3.0.5) extension [ 49 ] in Cytoscape to evaluate the overrepresentation of hub genes within hierarchical GO categories. The investigation was executed across three domains: biological process (BP), molecular function (MF), and cellular component (CC), with the objective of elucidating the functional roles and regulatory context of the identified gene sets. Multiple testing correction was executed by implementing the Benjamini–Hochberg false discovery rate methodology, employing a significant threshold of p < 0.05. Parallel to this, pathway enrichment analysis was conducted using the ClueGO/CluePedia (v 2.5.10) [ 50 ] plugin in Cytoscape. This approach comprised the integrated interrogation of both KEGG and Reactome pathway databases, with the intention of further elucidating the biological relevance of the hub genes within the framework of recognized signaling and metabolic pathways.
4.7. Violin Plot and ROC Curve Analysis
RMA-normalized gene expression data were aggregated by gene symbol, and both expression values and corresponding metadata were reshaped into a long format for the purpose of further analysis and visualization. The violin plots were constructed with the aid of the ggplot2 (v 4.0.0) and ggpubr (v 0.6.1) packages in R. Group differences were assessed using the Wilcoxon rank-sum test, followed by false discovery rate (FDR) correction to control for multiple comparisons. The discriminatory power of individual genes was evaluated using receiver operating characteristic (ROC) analysis implemented using the pROC package. Genes that demonstrated an area under the curve (AUC) > 0.9 and an adjusted p -value < 0.003 were considered to possess significant discriminatory power and were therefore selected for further investigation.
4.8. Nested Cross-Validation and Firth Logistic Regression
A nested CV framework was applied for evaluation of the predictive power of gene expression profiles to distinguish between MM_ONJ and control groups (Healthy or MM). The outer loop was used exclusively for model evaluation, while the inner loop performed feature selection. In each outer fold, statistical filtering was applied to the training data to identify significantly differentially expressed genes. For this, a two-sample t-test followed by false discovery rate (FDR) adjustment was used. Then, as predictors, the top five genes with the lowest adjusted p -values ( adj. p < 0.05) were selected. The logistf (v 1.26.1) package enabled us to employ these features to train Firth-penalized logistic regression models to address the challenges posed by limited data samples and mitigate potential biases in parameter estimation. Subsequently, an evaluation of the model’s performance was conducted by employing out-of-fold predictions from the outer CV loop. The ROC curves were generated from these predictions using the pROC (v 1.19.0.1) package. The AUC was calculated for each comparison. To facilitate this comparison, ROC curves were generated for the Healthy vs. MM_ONJ and MM vs. MM_ONJ groups. The AUC values were reported for each curve.
4.9. PCA + Ridge Regression with Safe Stratified 0.632+ Bootstrap
A safe stratified 0.632+ bootstrap framework was implemented to reduce optimism bias in AUC estimation and to assess statistical significance against the null distribution. The pipeline combined Principal Component Analysis (PCA) with ridge-penalized logistic regression. For each comparison (Healthy vs. MM_ONJ and MM vs. MM_ONJ), the dataset was restricted to the fixed four-gene LASSO-selected panel specific to that comparison. Genes were intersected with available expression features, and the resulting expression matrix was transposed to samples × genes format. A total of 200 bootstrap iterations were performed. In each iteration, stratified sampling with replacement selected 63.2% of the samples for the training set, leaving the remainder as out-of-bag (OOB) test samples. Safety constraints included the following: 1. Both outcome classes present in both training and test sets, 2. ≥2 samples per class in the training set before inner cross-validation, 3. Inner 3-fold stratified cross-validation with ≥2 samples per class in each fold. PCA was fitted only on the training set, retaining the minimum number of components required to explain ≥95% of variance, capped at three components. These components were used to train a ridge-penalized logistic regression model (cv.glmnet, α = 0). A safe CV routine was used, defaulting to “deviance” as the type.measure if “auc” was unstable. The following AUC metrics were recorded: Apparent AUC (training set), OOB AUC (test set), and 0.632+-corrected AUC, calculated with error = 1 − AUC and no-information error = 0.5. To assess significance, 200 stratified permutation iterations were performed. In each iteration, outcome labels were randomly shuffled while preserving the feature matrix, and the same PCA + ridge pipeline was applied to generate a null AUC distribution. The empirical p -value was computed as the proportion of permutation AUCs greater than or equal to the median 0.632+ AUC from the bootstrap. The 95% confidence interval (CI) for OOB AUCs was calculated from empirical quantiles. Cohen’s d effect size was computed between bootstrap and permutation AUC distributions.
5. Conclusions
In summary, the integrative pipeline that was utilized in this study yielded a convergent molecular portrait of MRONJ. This condition is driven by immune-inflammatory dysregulation, with stress responses hypothesized as critical modifiers. These findings provide both biomarker hypotheses and mechanistic insights that align with the intention to link molecular discovery and translational medicine to facilitate biomarker-guided risk prediction in MRONJ.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Supplementary Materials
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ijms262411788/s1 .
Author Contributions
Conceptualization, G.L. and I.T.; Data curation G.L.; Writing—Original Draft Preparation, G.L.; Writing—Review & Editing, G.L. and I.T.; Visualization, G.L.; Supervision G.L., Funding Acquisition, I.T.; Project Administration, I.T. and J.S. All authors have read and agreed to the published version of the manuscript.
Institutional Review Board Statement
Not applicable.
Informed Consent Statement
Not applicable.
Data Availability Statement
The data presented in this study are available in NCBI GEO database at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE7116 , accessed on 10 June 2025, reference number GSE7116 . These data were derived from the following resources available in the public domain: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE7116 .
Conflicts of Interest
The authors declare no conflicts of interest.
References
1. Ruggiero S.L. Dodson T.B. Aghaloo T. Carlson E.R. Ward B.B. Kademani D. American Association of Oral and Maxillofacial Surgeons’ Position Paper on Medication-Related Osteonecrosis of the Jaws—2022 Update J. Oral. Maxillofac. Surg. 2022 80 920 943 10.1016/j.joms.2022.02.008 35300956
2. Migliorati C.A. Oral Complications in Cancer Patients—Medication-Related Osteonecrosis of the Jaw (MRONJ) Front. Oral Heal. 2022 3 866871 10.3389/froh.2022.866871 PMC9086704 35558547
3. Kim C.M. Lee S. Hwang W. Son E. Kim T.W. Kim K. Kim Y.H. Obesity and Periodontitis: A Systematic Review and Updated Meta-Analysis Front. Endocrinol. 2022 13 999455 10.3389/fendo.2022.999455 36353241 PMC9637837
4. Marx R.E. Pamidronate (Aredia) and Zoledronate (Zometa) Induced Avascular Necrosis of the Jaws: A Growing Epidemic J. Oral Maxillofaci. Surg. 2003 61 1115 1117 10.1016/S0278-2391(03)00720-1 12966493
5. Bassan Marinho Maciel G. Marinho Maciel R. Linhares Ferrazzo K. Cademartori Danesi C. Etiopathogenesis of Medication-Related Osteonecrosis of the Jaws: A Review J. Mol. Med. 2024 102 353 364 10.1007/s00109-024-02425-9 38302741
6. Allen M.R. Burr D.B. The Pathogenesis of Bisphosphonate-Related Osteonecrosis of the Jaw: So Many Hypotheses, So Few Data J. Oral Maxillofaci. Surg. 2009 67 61 70 10.1016/j.joms.2009.01.007 19371816
7. Kim J.-M. Lin C. Stavre Z. Greenblatt M.B. Shim J.-H. Osteoblast-Osteoclast Communication and Bone Homeostasis Cells 2020 9 2073 10.3390/cells9092073 32927921 PMC7564526
8. Kusumoto J. Muraki Y. Iwata E. Matsumura M. Furudoi S. Akashi M. Involvement of Resistant Bacteria in the Severity of Refractory Osteonecrosis of the Jaw Clin. Oral Investig. 2025 29 454 10.1007/s00784-025-06547-3 PMC12441049 40956356
9. George E.L. Truesdell S.L. Magyar A.L. Saunders M.M. The Effects of Mechanically Loaded Osteocytes and Inflammation on Bone Remodeling in a Bisphosphonate-Induced Environment Bone 2019 127 460 473 10.1016/j.bone.2019.07.008 31301402
10. Bullock G. Miller C.A. McKechnie A. Hearnden V. A Review into the Effects of Pamidronic Acid and Zoledronic Acid on the Oral Mucosa in Medication-Related Osteonecrosis of the Jaw Front. Oral Heal. 2022 2 822411 10.3389/froh.2021.822411 PMC8865370 35224540
11. Kuroshima S. Yamashita J. Chemotherapeutic and Antiresorptive Combination Therapy Suppressed Lymphangiogenesis and Induced Osteonecrosis of the Jaw-like Lesions in Mice Bone 2013 56 101 109 10.1016/j.bone.2013.05.013 23727433
12. Konarski W. Poboży T. Konarska K. Śliwczyński A. Kotela I. Krakowiak J. Exploring the Impact of Novel Anti-Cancer Therapies on Jaw Osteonecrosis and Other Bones: A Comprehensive Review J. Clin. Med. 2024 13 1889 10.3390/jcm13071889 38610654 PMC11012550
13. Zhou F. Zhang G. Wu Y. Xiong Y. Inflammasome Complexes: Crucial Mediators in Osteoimmunology and Bone Diseases Int. Immunopharmacol. 2022 110 109072 10.1016/j.intimp.2022.109072 35978515
14. Balachandran K. Ramli R. Karsani S.A. Abdul Rahman M. Identification of Potential Biomarkers and Small Molecule Drugs for Bisphosphonate-Related Osteonecrosis of the Jaw (BRONJ): An Integrated Bioinformatics Study Using Big Data Int. J. Mol. Sci. 2023 24 8635 10.3390/ijms24108635 37239981 PMC10218302
15. Chen C. Hu J.-F. Liu B.-Y. Sun Y.-M. Xiao B.-L. Transcription Factor ZNF207 Drives Aerobic Glycolysis and Facilitates Malignancy Progression in Hepatocellular Carcinoma Cell. Signal. 2025 135 112009 10.1016/j.cellsig.2025.112009 40684964
16. Jóźwik M. Sidorkiewicz I. Wojtkiewicz J. Sulkowski S. Semczuk A. Jóźwik M. Selecting Optimal Housekeeping Genes for RT-qPCR in Endometrial Cancer Studies: A Narrative Review Int. J. Mol. Sci. 2025 26 8610 10.3390/ijms26178610 40943534 PMC12429465
17. Lorenzo-Pouso A.I. Bravo S.B. Carballo J. Chantada-Vázquez M.D.P. Bagán J. Bagán L. Chamorro-Petronacci C.M. Conde-Amboage M. López-López R. García-García A. Quantitative Proteomics in Medication-related Osteonecrosis of the Jaw: A Proof-of-concept Study Oral Dis. 2023 29 2117 2129 10.1111/odi.14201 35377498
18. Laputková G. Talian I. Schwartzová V. Medication-Related Osteonecrosis of the Jaw: A Systematic Review and a Bioinformatic Analysis Int. J. Mol. Sci. 2023 24 16745 10.3390/ijms242316745 38069068 PMC10706386
19. Walton K. Grogan T.R. Eshaghzadeh E. Hadaya D. Elashoff D.A. Aghaloo T.L. Tetradis S. Medication Related Osteonecrosis of the Jaw in Osteoporotic vs Oncologic Patients—Quantifying Radiographic Appearance and Relationship to Clinical Findings Dentomaxillofacial Radiol. 2019 48 20180128 10.1259/dmfr.20180128 29952657 PMC6398905
20. Shi Y. Li M. Yu Y. Zhou Y. Zhang W. Hua H. Wang S. Stress Response in Periodontal Ligament Stem Cells May Contribute to Bisphosphonate—Associated Osteonecrosis of the Jaw: A Gene Expression Array Analysis Mol. Med. Rep. 2020 22 2043 2051 10.3892/mmr.2020.11276 32705175 PMC7411417
21. Thiel Y. Ghayor C. Lindhorst D. Essig H. Weber F. Rücker M. Schumann P. Antimicrobial Peptide Gene Expression in Medication-Related Osteonecrosis of the Jaw Pathol.-Res. Pract. 2020 216 153245 10.1016/j.prp.2020.153245 33065485
22. Kalyan S. Chandrasekaran V. Quabius E.S. Lindhorst T.K. Kabelitz D. Neutrophil Uptake of Nitrogen-Bisphosphonates Leads to the Suppression of Human Peripheral Blood Γδ T Cells Cell. Mol. Life Sci. 2014 71 2335 2346 10.1007/s00018-013-1495-x 24162933 PMC11114071
23. Zheng J. Yao Z. Xue L. Wang D. Tan Z. The Role of Immune Cells in Modulating Chronic Inflammation and Osteonecrosis Front. Immunol. 2022 13 1064245 10.3389/fimmu.2022.1064245 36582244 PMC9792770
24. Stockmann P. Wehrhan F. Schwarz-Furlan S. Stelzle F. Trabert S. Neukam F.W. Nkenke E. Increased Human Defensine Levels Hint at an Inflammatory Etiology of Bisphosphonate-Associated Osteonecrosis of the Jaw: An Immunohistological Study J. Transl. Med. 2011 9 135 10.1186/1479-5876-9-135 21843332 PMC3163206
25. Zheng Y. Wang X. He Y. Chen S. He L. Zhang Y. Exosomes from Adipose-Derived Mesenchymal Stromal Cells Prevent Medication-Related Osteonecrosis of the Jaw by Inhibiting Macrophage M1 Polarization and Pyroptosis Int. J. Nanomed. 2024 19 12675 12693 10.2147/IJN.S482849 PMC11608005 39619056
26. De Oliveira C. Brizeno L. De Sousa F. Mota M. Alves A. Osteonecrosis of the Jaw Induced by Receptor Activator of Nuclear Factor-Kappa B Ligand (Denosumab)—Review Med. Oral Patol. Oral Cir. Bucal 2016 21 e431 e439 10.4317/medoral.21044 26827069 PMC4920456
27. Giryes S. Wong C. Bridgewood C. Harland M. Altaie A. McDermott N. Abacar K. Rao A. Khan A. McMillan T. Paradoxical Activation of Entheseal Myeloid Cells by JAK1 and Tyk2 Inhibitors via Interleukin—10 Antagonism Arthritis Rheumatol. 2025 77 1373 1386 10.1002/art.43210 40320905 PMC12479184
28. Koizumi K. Saitoh Y. Minami T. Takeno N. Tsuneyama K. Miyahara T. Nakayama T. Sakurai H. Takano Y. Nishimura M. Role of CX3CL1/Fractalkine in Osteoclast Differentiation and Bone Resorption J. Immunol. 2009 183 7825 7831 10.4049/jimmunol.0803627 19923448
29. Zupan J. Jeras M. Marc J. Osteoimmunology and the Influence of Pro-Inflammatory Cytokines on Osteoclasts Biochem. Med. 2013 23 43 63 10.11613/BM.2013.007 23457765 PMC3900089
30. Spinelli F.R. Colbert R.A. Gadina M. JAK1: Number One in the Family; Number One in Inflammation? Rheumatology 2021 60 ii3 ii10 10.1093/rheumatology/keab024 33950229 PMC8599761
31. Damerau A. Gaber T. Ohrndorf S. Hoff P. JAK/STAT Activation: A General Mechanism for Bone Development, Homeostasis, and Regeneration Int. J. Mol. Sci. 2020 21 9004 10.3390/ijms21239004 33256266 PMC7729940
32. Murakami K. Kobayashi Y. Uehara S. Suzuki T. Koide M. Yamashita T. Nakamura M. Takahashi N. Kato H. Udagawa N. A Jak1/2 Inhibitor, Baricitinib, Inhibits Osteoclastogenesis by Suppressing RANKL Expression in Osteoblasts in Vitro PLoS ONE 2017 12 e0181126 10.1371/journal.pone.0181126 28708884 PMC5510865
33. Martens A. Hertens P. Priem D. Rinotas V. Meletakos T. Gennadi M. Van Hove L. Louagie E. Coudenys J. De Muynck A. A20 Controls RANK—Dependent Osteoclast Formation and Bone Physiology Embo Rep. 2022 23 e55233 10.15252/embr.202255233 36194667 PMC9724664
34. Bao Y. Tong C. Xiong X. CXCL3: A Key Player in Tumor Microenvironment and Inflammatory Diseases Life Sci. 2024 348 122691 10.1016/j.lfs.2024.122691 38714265
35. Sawant K.V. Sepuru K.M. Lowry E. Penaranda B. Frevert C.W. Garofalo R.P. Rajarathnam K. Neutrophil Recruitment by Chemokines Cxcl1/KC and Cxcl2/MIP2: Role of Cxcr2 Activation and Glycosaminoglycan Interactions J. Leukoc. Biol. 2021 109 777 791 10.1002/JLB.3A0820-207R 32881070 PMC8296306
36. Yan M. Komatsu N. Muro R. Huynh N.C.-N. Tomofuji Y. Okada Y. Suzuki H.I. Takaba H. Kitazawa R. Kitazawa S. ETS1 Governs Pathological Tissue-Remodeling Programs in Disease-Associated Fibroblasts Nat. Immunol. 2022 23 1330 1341 10.1038/s41590-022-01285-0 35999392
37. Yang Y. Han X. Sun L. Shao F. Yin Y. Zhang W. ETS Transcription Factors in Immune Cells and Immune-Related Diseases Int. J. Mol. Sci. 2024 25 10004 10.3390/ijms251810004 39337492 PMC11432452
38. Roato I. Mauceri R. Notaro V. Genova T. Fusco V. Mussano F. Immune Dysfunction in Medication-Related Osteonecrosis of the Jaw Int. J. Mol. Sci. 2023 24 7948 10.3390/ijms24097948 37175652 PMC10177780
39. Zhang W. Gao L. Ren W. Li S. Zheng J. Li S. Jiang C. Yang S. Zhi K. The Role of the Immune Response in the Development of Medication-Related Osteonecrosis of the Jaw Front. Immunol. 2021 12 606043 10.3389/fimmu.2021.606043 33717086 PMC7947359
40. Tetradis S. Allen M.R. Ruggiero S.L. Pathophysiology of Medication—Related Osteonecrosis of the Jaw—A Minireview JBMR Plus 2023 7 e10785 10.1002/jbm4.10785 37614299 PMC10443081
41. Xiao Y. Xie X. Chen Z. Yin G. Kong W. Zhou J. Advances in the Roles of ATF4 in Osteoporosis Biomed. Pharmacother. 2023 169 115864 10.1016/j.biopha.2023.115864 37948991
42. Muehlebach M.E. Holstein S.A. The Role of the Unfolded Protein Response Pathway in Bone Homeostasis and Potential Therapeutic Target in Cancer-Associated Bone Disease Bone Res. 2025 13 76 10.1038/s41413-025-00457-6 40877255 PMC12394429
43. Tamaoka J. Takaoka K. Hattori H. Ueta M. Maeda H. Yamamura M. Yamanegi K. Noguchi K. Kishimoto H. Osteonecrosis of the Jaws Caused by Bisphosphonate Treatment and Oxidative Stress in Mice Exp. Ther. Med. 2019 17 1440 1448 10.3892/etm.2018.7076 30680026 PMC6327453
44. Raje N. Woo S.-B. Hande K. Yap J.T. Richardson P.G. Vallet S. Treister N. Hideshima T. Sheehy N. Chhetri S. Clinical, Radiographic, and Biochemical Characterization of Multiple Myeloma Patients with Osteonecrosis of the Jaw Clin. Cancer Res. 2008 14 2387 2395 10.1158/1078-0432.CCR-07-1430 18413829
45. Langfelder P. Horvath S. WGCNA: An R Package for Weighted Correlation Network Analysis BMC Bioinform. 2008 9 559 10.1186/1471-2105-9-559 19114008 PMC2631488
46. Shannon P. Markiel A. Ozier O. Baliga N.S. Wang J.T. Ramage D. Amin N. Schwikowski B. Ideker T. Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks Genome Res. 2003 13 2498 2504 10.1101/gr.1239303 14597658 PMC403769
47. Doncheva N.T. Morris J.H. Gorodkin J. Jensen L.J. Cytoscape StringApp: Network Analysis and Visualization of Proteomics Data J. Proteome Res. 2019 18 623 632 10.1021/acs.jproteome.8b00702 30450911 PMC6800166
48. Chin C.-H. Chen S.-H. Wu H.-H. Ho C.-W. Ko M.-T. Lin C.-Y. cytoHubba: Identifying Hub Objects and Sub-Networks from Complex Interactome BMC Syst. Biol. 2014 8 S11 10.1186/1752-0509-8-S4-S11 25521941 PMC4290687
49. Maere S. Heymans K. Kuiper M. BiNGO: A Cytoscape Plugin to Assess Overrepresentation of Gene Ontology Categories in Biological Networks Bioinformatics 2005 21 3448 3449 10.1093/bioinformatics/bti551 15972284
50. Bindea G. Mlecnik B. Hackl H. Charoentong P. Tosolini M. Kirilovsky A. Fridman W.-H. Pagès F. Trajanoski Z. Galon J. ClueGO: A Cytoscape Plug-in to Decipher Functionally Grouped Gene Ontology and Pathway Annotation Networks Bioinformatics 2009 25 1091 1093 10.1093/bioinformatics/btp101 19237447 PMC2666812