Study design

We used summary-level data from GWAS to conduct a two-sample MR analysis to investigate the causal relationship between circulating inflammatory proteins and the risk of OM. Figure 1 shows a schematic diagram of the study design. To accurately establish the causal effect, MR relies on three key assumptions: (1) the association assumption: instrumental variables (IVs) should strongly associate with the exposure; (2) independence assumption: the included IVs are unrelated to confounding factors associated with the exposure-outcome link; (3) exclusivity assumption: IVs only modify the outcome through their effect on the exposure. This study strictly adhered to the guidelines of Strengthening the Reporting of Observational Studies in Epidemiology using Mendelian Randomization (STROBE-MR) to report the study findings (Supplementary Table SI). All participants in the included studies were of European ancestry, avoiding population stratification-related statistical bias. The data were all obtained from publicly available databases, and each cohort received ethical approval during the study period; thus no additional ethical review was required.

Figure 1

Overall study design regarding the association between 91 circulating inflammatory proteins and osteomyelitis

Data sources

The latest data on circulating inflammatory proteins are derived from a study by Zhao et al. [18]. The researchers recruited 14,824 participants of European descent and performed whole-genome protein quantitative trait loci (pQTL) mapping on 91 plasma proteins measured from 11 cohorts. Inflammatory proteins were generated by measuring genome-wide genetic data and plasma proteomics data with the Olink Target Inflammation immunoassay panel. GWAS analysis within each cohort was performed by applying an additive genetic association model based on linear regression, and the impact of inflammatory protein was reported as a change in inverse-rank normalized protein level per dosage of the effect allele [19]. For further details regarding the research, please consult the original literature [18].

The summary-level GWAS for OM was obtained from the study by Hamilton et al. and the FinnGen R11 database [20, 21]. The FinnGen R11 database is a large biological resource project aimed at collecting and analyzing genetic data of the Finnish population to study the influence of genetic and environmental factors on health and disease [20]. R11 is the updated public data set organized by the project for June 2024. In the discovery phase, the FinnGen R11 database included 431,951 cases of European ancestry (2125 cases of OM and 429,826 controls), and other information can be accessed through the FinnGen website (https://www.finngen.fi/en). In the replication phase, we used the pooled summary statistics from the IEU OpenGWAS project of the UK Biobank (https://gwas.mrcieu.ac.uk/). The UK Biobank is a large biomedical cohort study database initiated and established by the UK government. The database recruited 500,000 participants from across the UK and collected genetic, lifestyle, and health data from a large sample of individuals through various forms, including biological samples and questionnaires [21]. The study included a total of 4836 cases of OM and 481,648 controls. The sources of exposure and outcome data and the characteristics of the data are presented in Table I.

Selection of genetic instrumental variables

To meet the MR assumptions, SNPs (R² = 0.001 and genetic distance = 10,000 kb) that were strongly and independently associated with exposure at the genome-wide level (p < 5 × 10^–8) were included in the discovery phase, a strict threshold designed to prioritize the reliability of instrumental variables and minimize the risk of weak instrument bias [22]. In the replication phase, we used a less stringent exclusion criterion for linkage disequilibrium (R² = 0.01 and genetic distance = 10,000 kb) to enhance our method and statistical power [23]. This adjustment strategy has been validated in other studies to balance the rigor of the discovery phase and the stability of the replication phase, and eventually include more relevant SNPs to enhance the reproducibility of the results. To prevent allelic confounding and bidirectional causality between circulating inflammatory proteins and bone marrow inflammation risk, we excluded inverted and incompatible SNPs. Furthermore, we calculated the F statistic to use the formula to assess the strength of IVs (R² = 2 × EAF × (1 – EAF) × beta²; F = R² × (N – 2)/(1 – R²)), and excluded SNPs with an F statistic < 10 (to exclude weakly associated SNPs) [24].

Statistical analysis

MR analysis was conducted using the TwoSampleMR package in the R 4.4.1 version software. The MR-Egger, weighted median, inverse variance weighted (IVW), simple mode, and weighted mode methods were used as the default analysis methods. Among these, IVW is widely accepted as the primary method due to its calculation of a weighted average of all IV effect sizes, providing relatively reliable results [25]. Other methods are only used as supplementary results to confirm the general direction due to poorer statistical efficiency. The causal relationship between exposure and outcome is evaluated based on p-values (p < 0.05 indicates statistical significance). The IVW method and Egger regression are used to detect heterogeneity [26]. Cochran ‘s Q test is used to assess heterogeneity between IVs. The intercept in the Egger regression is an effective indicator of whether the MR analysis results are affected by horizontal pleiotropy (p < 0.05 indicates significant heterogeneity). The “leave-one-out” method was employed to conduct sensitivity analysis, whereby each SNP was sequentially eliminated and the effect size of the remaining SNPs was calculated. If the exclusion of any single SNP significantly impacts the results, it indicates the presence of sensitivity [27]. Meta-analysis was conducted using Stata (17.0 version), and the OR and its related 95% confidence interval were combined.

Phenome-wide association analysis

To further assess whether inflammatory protein factors have beneficial or harmful effects on other traits, as well as to detect pleiotropy not identified by the MR-Egger test, we conducted a PheWAS using the AstraZeneca PheWAS Portal database (https://azphewas.com/) [28]. The research team studied the relationships between rare protein-coding variants and 17,361 binary and 1,419 quantitative phenotypes using exome sequencing data from 269,171 European-ancestry UK Biobank participants. PheWAS results can be interpreted as genetic predictions of protein expression associated with specific diseases or traits. The significant p-value threshold (1e–8) was established according to the study of Wang et al. [28] and corresponds to a false positive rate of 0.1%. The suggestive threshold (1e–6) is adjusted by a single phenotypic collapsing model to preserve conservative control for p < (0.05/18500 genes) [19].

Results of instrumental variable screening

In the MR analysis, the effect of SNPs on exposure and the effect of SNPs on outcomes must correspond to the same allele. In the discovery phase, after reliable and independent (R² = 0.001 and genetic distance = 10,000 kb) adjustment of the coordinated exposure and outcome data, duplicate and negatively correlated SNPs were removed, resulting in the identification of 263 usable SNPs (Supplementary Table SII). In the replication phase, 412 usable SNPs were identified with an R² = 0.01 and genetic distance = 10,000 kb criterion (Supplementary Table SV).

MR discovery analysis of 91 inflammatory proteins on OM

In the discovery phase, a total of 4 circulating inflammatory proteins were identified as potentially causally associated with OM using IVW-based MR analysis (Figure 2). The SNP information and the MR analysis process for these 4 circulating inflammatory proteins are presented in Supplementary Tables SII–SIV. Osteoprotegerin (OPG) (OR = 1.484, 95% CI = 1.139–1.933, p = 0.003), monocyte chemoattractant protein-4 (MCP-4) (OR = 1.223, 95% CI = 1.029–1.453, p = 0.022), adenosine deaminase (ADA) (OR = 1.209, 95% CI = 1.056–1.384, p = 0.006), and CCL4 (OR = 1.121, 95% CI = 0.997–1.261, p = 0.05) all showed statistically significant results (p < 0.05) in the MR analysis. Consistent with the results of the IVW method, the weighted median, simple mode, and weighted mode produced directionally consistent effect estimates, confirming the stability of our discovery phase research results (Figure 3). The sensitivity analysis results showed no significant heterogeneity or horizontal pleiotropy (Supplementary Tables SVIII and SIX). In the visualization results, we plotted scatter plots to demonstrate the impact of individual SNPs on the causal estimate. (Supplementary Figure S1). Furthermore, the leave-one-out plot analysis for the 4 circulating inflammatory proteins and OM also clearly showed that the MR analysis results were stable and not affected by any single SNP (Supplementary Figure S2).

Figure 2

Circular heat map of 91 circulating inflammatory proteins on osteomyelitis MR Results in the discovery cohort

Figure 3

Causal relationship between 4 inflammatory proteins and osteomyelitis in the discovery cohort. The effect of the four inflammatory proteins on osteomyelitis is represented by the OR value

MR replication analysis of 91 inflammatory proteins on OM

In the replication phase, the results obtained by IVW as the main MR method show that 6 inflammatory proteins have a potential causal association with OM (Figure 4). The SNP information and the process of MR analysis for the six circulating inflammatory proteins are presented in Supplementary Tables SV–SVII. Monocyte chemoattractant protein-3 (MCP-3) (OR = 0.822, 95% CI = 0.692–0.976, p = 0.025), interleukin-15 receptor subunit alpha (IL15RA) (OR = 1.088, 95% CI = 1.001–1.183, p = 0.047), C-X-C motif chemokine 5 (CXCL5) (OR = 0.865, 95% CI = 0.759–0.986, p = 0.03), CCL4 (OR = 1.108, 95% CI = 1.023–1.2, p = 0.012), C-C motif chemokine 19 (CCL19) (OR = 0.843, 95% CI = 0.723–0.984, p = 0.03), and artemin (ARTN) had a potential causal association with osteomyelitis (OR = 1.471, 95% CI = 1.031–2.098, p = 0.034). In the other four default methods, the same trend was observed despite the lack of statistical significance (Figure 5). Cochran’s Q test showed that the heterogeneity of all exposure-outcome analyses was not significant, indicating that the effect size of SNPs was consistent, and there was no need to further exclude specific SNPS. The MR-Egger intercept test again did not reveal significant horizontal pleiotropy (Supplementary Tables SVIII and SIX) In addition, the overall effect remained stable after excluding individual SNP results in the “leave-one-out” analysis.

Figure 4

Circular heat map of 91 circulating inflammatory proteins on osteomyelitis MR Results in the replication cohort

Figure 5

Causal relationship between inflammatory proteins and osteomyelitis in the discovery cohort. The effect of the six inflammatory proteins on osteomyelitis is represented by the OR value

Meta-analysis

Using the aforementioned data, we found that the statistical significance of CCL4 being associated with OM remained significant after validation cohort analysis. Sensitivity analysis indicated no data heterogeneity or horizontal pleiotropy. We conducted a meta-analysis using a random effects model for the two cohorts, with the combined result showing an OR_IVW = 1.11 (95% CI = 1.04–1.19) (Figure 6). The genetic prediction of a one standard deviation increase in CCL4 was associated with a 1.11-fold increased risk of OM, indicating a robust causal association.

Figure 6

Meta-analysis of the causal relationship between CCL4 and osteomyelitis. Random effects model pooled effect sizes for IVW outcomes in both cohorts are shown

Phenotype-wide association analysis

As shown in Supplementary Figures S3–S12, apart from IL15RA, which provides suggestive support for continuous traits in inflammatory diseases, four other drug targets were significantly associated with any phenotype. This suggests that the potential side effects of drugs targeting these targets and the pleiotropy present in these genes may be minimal, further supporting the reliability of the study findings.

Association with previous studies

To our knowledge, this is the first MR study to investigate the causal relationship between 91 circulating inflammatory proteins and the risk of OM. In fact, multiple countries and disease organizations classify OM as an imminent risk for post-traumatic fracture or injury prognosis and consistently state that elucidating its inflammatory development mechanisms and epidemiological characteristics is crucial for reducing its threat to human health [29]. A recent study confirmed that CCL4 is closely related to inflammatory responses and blood glucose level regulation, making it a risk factor or pathogenic factor for cardiovascular diseases [30]. Gan et al. also reported that an increase in blood CCL4 was associated with an increased risk of hospitalization for COVID-19, which may be related to its role in the overexpression of inflammatory factors and lung inflammation injury. These findings support the importance of this inflammatory protein in cardiovascular and respiratory inflammatory diseases and indirectly support the results of this study. Specifically, CCL4 mediates inflammatory responses in osteomyelitis by binding to the co-receptor CCR5, as well as by activating the mTORC1 pathway. After binding to the receptor CCR5, CCL4 recruits macrophages and CD8+ T cells to the site of infection, promotes the release of pro-inflammatory factors (such as IL-6 and TNF-α), and intensifies local inflammation [31]. In chronic osteomyelitis, CCL4 inhibits autophagy by activating the mTORC1 pathway, leading to the accumulation of inflammatory factors and enhanced bone resorption [32].

Possible explanations for the causality

Chronic OM is typically caused by bacterial infections, with SA being the most common pathogen. These infectious microorganisms disrupt the microecological balance of bone formation and resorption, primarily by being ingested and internalized by osteoblasts, persisting within them, recruiting immune cells to kill osteoblasts, and increasing the number of osteoclasts.

Osteoblasts ingest and internalize SA, which can secrete various inflammatory cytokines and chemokines [33]. While there have been extensive discussions on inflammatory cytokines such as tumor necrosis factor-α (TNF-α), interleukin (IL)-6, and IL-18 in previous studies, the results have been mixed [33–35]. On the one hand, only a few studies have reported the presence of these cytokines or their coding genes with increased expression in osteoblasts infected with Staphylococcus aureus. On the other hand, only increased expression of inflammatory cytokines was detected at the mRNA level in different studies, without any changes detected at the protein level [36]. Furthermore, due to the specificity of in vitro studies, it is difficult to simulate the actual situation of osteoblasts in vivo. These limitations have led to slow progress in this field, and the role of inflammatory cytokines secreted by osteoblasts in the inflammatory process of OM is still unclear.

It is worth noting that osteoblasts infected with Staphylococcus aureus can also secrete various chemokines to recruit and activate innate and adaptive immune cells. The chemokines belong to two families (CXC and CC), which are key signals for circulating cell migration and localization to various tissues and play an important role in bone metabolism [37]. Osteoblast-secreted chemokines, such as IL-6, CCL2, CCL3, and CCL5, participate in the recruitment and activation of macrophages [11, 38]. CCL4, CCL19, and CXCL5 may also participate in this process, leading to cell proliferation and differentiation, cell polarization, cell apoptosis and degranulation, nitric oxide induction, and reactive oxygen species production, as well as actin reorganization, through the Jak-STAT signaling pathway, MAPK signaling pathway, PLC/PKC signaling pathway, etc.

Osteoblasts produce NF-κB activator ligand (RANK-L) and OPG. RANK-L can bind to the ligand RANK on monocytes/macrophages, promoting their differentiation into osteoblasts. OPG is a decoy receptor that can compete with RANK for binding to RANK-L, thereby inhibiting the differentiation of osteoclasts. When Staphylococcus aureus infects osteoblasts, RANK-L synthesis increases and OPG synthesis decreases, reducing the inhibition of osteoclast differentiation. This change also indirectly leads to an increase in the number of osteoclasts. In addition to RANK-L, cellular factors such as IL-15 and IL-17 have also been shown to be involved in the activation and renewal of osteoclasts. The absence of an IL-15 signal impairs the activity of osteoclasts, and IL15RA is an important part of the IL-15 pro-inflammatory signal. The combination of IL-15 and IL-15RA activates the JAK/STAT signaling pathway, and promotes the survival and proliferation of NK cells and CD8+ T cells. These cells play a key role in the resistance to bacterial infection, and their insufficient activity may lead to the persistence of S. aureus in bone tissue and aggravate the progression of OM [39, 40]. In the pathological process of OM, MCP-3 mainly binds to its receptor CCR1/CCR3, which induces monocytes and macrophages to migrate to the infection site, releases inflammatory factors such as IL-6 and TNF-α, and aggravates the local inflammatory response and bone destruction. Furthermore, in vitro studies have shown that MCP-3 can exacerbate bone resorption by stimulating RANKL production [12]. However, in the study by Votta et al., no chemokine, including MCP-3 and MCP-4, showed chemotactic activity toward primary osteoblasts or osteoblasts derived from human bone grafts [41]. In summary, the mechanisms by which chemokines and inflammatory factors participate in bone metabolism are complex and diverse, and current research is still inadequate. Future studies should further explore mechanisms underlying the progression of OM from the perspective of bone immunology.

It is easy to understand that ADA can convert adenosine to inosine, which participates in purine metabolism. When the count of CD4 T cells is high, ADA activity increases. In addition, inflammatory bone pain associated with bone pathology or disease is an important clinical problem in which ARTN seems to play an important role. By binding to the receptor GFRα3, ARTN activates TRPV1 channels in the dorsal root ganglion (DRG) and enhances nociceptive signaling, which may lead to mechanical and thermal hyperalgesia in osteomyelitis patients [42]. In addition, ARTN may prolong the inflammatory cycle by inhibiting autophagy or promoting the sustained release of inflammatory factors such as IL6, resulting in difficult infection control [43]. However, current studies have not clarified the direct correlation between ARTN and TRPV1 channels and osteoclast activation. In the future, animal models are needed to verify the changes in the expression of ARTN in osteomyelitis, which may provide benefits for the treatment of osteomyelitis.

Limitations of the study

The study has the following limitations. First, the completeness of GWAS datasets in both cohorts, such as the diagnosis of OM and the completeness of medical records, directly affects the accuracy of the study. Second, the study is limited to individuals of European ancestry, and caution is needed when interpreting the results for other populations. Third, although we conducted a series of sensitivity analyses and PheWAS, it should be noted that heterogeneity cannot be completely eliminated due to inherent methodological biases of MR and differences in genetic data. Finally, further research is needed to investigate the effects of gender, lifestyle, or nutrients on the associations, which should be validated in larger GWAS samples or clinical settings.

In conclusion, our study findings support a potential causal association between 9 circulating inflammatory proteins and OM, with a particular emphasis on the role of CCL4 in the immune progression of osteoblasts. However, the current results can only be regarded as supporting evidence for a potential causal relationship, which needs to be further verified by combining experimental and clinical studies and in multi-ethnic populations.