Exploring the Association Between Immune Cell Phenotypes and Osteoporosis Mediated by Inflammatory Cytokines: Insights from GWAS and Single-Cell Transcriptomics

Background

Osteoporosis is a common disease characterized by low bone mass and destruction of bone structure, leading to impaired bone strength and increased fracture risk.¹ The prevalence of osteoporosis in the world’s older adults was 21.7%; and was 24.3%, 16.7%, and 11.5% in Asia, Europe, and the United States, respectively, with Asia having the highest prevalence.² As the global population continues to age, the economic burden of osteoporosis is expected to increase, highlighting the need to identify risk factors for its development.

The immune system is a fascinating world, composed of interrelated organs, tissues, cells, and molecules. Immune cells are an essential component of the human immune system, playing a crucial role in maintaining health and defending against diseases. They also have associations with various systems within the body, significantly influencing a range of life processes.³ The actions of inflammatory mediators are highly complex, as they can function independently or as part of a highly organized, hierarchical cascade, extensively participating in various processes within the human body.⁴

Previous studies have indicated that there is a complex interplay between immune cells and bones, mediating a series of osteogenic and osteoclastic processes.⁵ This interaction plays a significant role in the formation and progression of osteoporosis.⁶ Additionally, inflammatory markers also play a significant role in osteoporosis, with a range of pro-inflammatory cytokines proven to stimulate osteoclastic bone resorption, including TNFα, IL-6, IL-11, IL-15, and IL-17A.^7–9 Many studies have also revealed some possible associations by which immune cells and inflammatory markers are involved in the development of osteoporosis. In observational studies, Tang et al found that increased systemic immune-inflammatory index (SII) is associated with an increased risk of low bone mineral density (BMD) and osteoporosis, and other inflammatory markers, particularly neutrophil-to-lymphocyte ratio (NLR) and the product of platelet count and neutrophil count (PPN), are inversely correlated with BMD and positively associated with osteoporosis risk.¹⁰ Di et al used Multivariable linear, Cox regression models, and mediation analysis to demonstrate that systemic inflammation and frailty phenotype are independently associated with an increased risk of osteoporosis and fractures.¹¹

However, current studies still have certain limitations. Firstly, traditional observational studies do not determine causation, and the association between immune, inflammation and osteoporosis reported in previous studies may still be interfered with by reverse causality and residual confounders.^12,13 Secondly, current evidence regarding the effects of immunity and inflammation on bone is limited and sometimes even controversial. Moreover, comprehensive and systematic studies evaluating the causal relationships among the three are lacking. Therefore, association between immunity, inflammation, and osteoporosis are not well understood, underscoring an urgent need for further clarification. Given the above limitations, leveraging the advantages of the Mendelian database, we conducted a study on immunity, inflammation, and bone to clarify the interactions among these factors and further elucidate the role and function of inflammation in bone regulation. To our knowledge, there have been no MR studies reporting on the effects of immune cell phenotypes and inflammatory markers on osteoporosis outcomes.

Mendelian randomization, as a widely used analytical method, aims to test the causal hypothesis between exposure factors and outcomes.¹⁴ It uses genetic variation, usually single nucleotide polymorphisms, as instrumental variables for hypothesizing risk factors.¹⁵ Because genetic variants are present at birth and remain stable throughout the lifespan, the conclusions of Mendelian randomization analyses are not disturbed by reverse causality and are less likely influenced by confounders.^16,17 Due to the high level of evidence for Mendelian randomization analysis, many scholars have used this method to explore the causal relationship between certain factors and diseases. Previous studies have shown that using genetic variants as instrumental variables to investigate the causal relationships between immunity, inflammation, and disease risk is feasible.^18–20 Therefore, we employ Mendelian randomization analysis to explore the causal relationship between immune cells and various inflammatory markers with the risk of osteoporosis.

Additionally, we analyzed single-cell RNA sequencing (scRNA-seq) data from an osteoporosis patient. This approach allows for a detailed characterization of individual cell types and their interactions, providing deeper insights into the mechanisms underlying osteoporosis.

Materials and Methods

Study Design

The application of Mendelian randomization is based on three key assumptions. The first is the relevance assumption, which posits that there is a robust and strong correlation between the genetic variants used as instrumental variables and the exposure of interest; the use of weak instrumental variables can introduce bias. The second is the independence assumption, which states that genetic variants are not associated with confounders that linked to both the exposure and the corresponding outcome. Lastly, the exclusivity assumption posits that there is no direct causal relationship between the genetic variants and the outcome, their connection is mediated solely through the exposure.^21,22 Based on the aforementioned principles, this study employs a two-sample Mendelian randomization approach to explore the potential causal relationship between immune cells and inflammatory factors with osteoporosis (Figure 1). Pleiotropy and heterogeneity analyses were conducted in the analysis to verify the robustness of the results. Subsequently, we analyzed single-cell sequencing data to investigate the biological functions of key cells. All studies involved have been approved by the relevant institutional review boards. We hope that our exploration will provide a relatively clear direction for future research and treatment.

Figure 1 Overview of the assumptions of the Mendelian randomization (MR) design and the study design.

Data Sources

We accessed Genome-Wide Association Study (GWAS) data within the scope of immunology to obtain information related to immune cells, with the summary statistics for each immune trait being publicly accessible through the GWAS Catalog, and the accession numbers ranging from GCST0001391 to GCST002121. We have included a total of 731 immune phenotypes, covering a variety of categories, providing ample data support for our research. All the data comes from 3757 Sardinians.²³ Data on inflammatory cytokines are sourced from two distinct studies. One study offers genome variant associations with 41 cytokines and growth factors in 8293 Finnish individuals.²⁴ The GWAS summary statistics for this study can be accessed and downloaded from the GWAS Catalog website (https://www.ebi.ac.uk/gwas/), with accession numbers ranging from GCST004420 to GCST004460. The other study has conducted a meta-analysis of 91 circulating inflammatory proteins from 11 cohorts, encompassing 14,824 participants of European descent.²⁵ The comprehensive GWAS summary statistics for each protein are also available for download from the GWAS Catalog website (https://www.ebi.ac.uk/gwas/), with accession numbers spanning from GCST90274758 to GCST90274848. The relevant data for outcome factors come from the FinnGen study. We used the R11 version of the FinnGen documentation and selected five groups of data, which are as follows: finngen_R11_DRUGADVERS_OSTEOPO (Drug-induced osteoporosis; Sample sizes: 453733, 346 cases and 453387 controls), finngen_R11_M13_OSTEOPOROSIS (Osteoporosis; Sample sizes: 438872, 9046 cases and 429826 controls), finngen_R11_OSTEOPOROSIS_FRACTURE_FG (Osteoporosis with pathological fracture; Sample sizes: 343460, 2085 cases and 341375 controls), finngen_R11_OSTPOPATFRACTURE (Drug-induced osteoporosis with pathological fracture; Sample sizes: 452082, 434 cases and 451648 controls) and finngen_R11_OSTPOPATFRCTURE_POSTEMENO (Postmenopausal osteoporosis with pathological fracture; Sample sizes: 254370, 1709 cases and 252661 controls). The data sources are: FinnGen. FinnGen Documentation of R11 release, 2024. Available at: https://finngen.gitbook.io/documentation/.

The diagnostic criteria for osteoporosis with pathological fracture, osteoporosis without pathological fracture, and secondary osteoporosis due to other classified diseases are categorized under M80, M81, and M82, respectively, in the International Classification of Diseases (ICD)-10. Dual-energy X-ray absorptiometry (DXA or DEXA) is the gold standard for diagnosing osteoporosis. A diagnosis of osteoporosis is established when bone mineral density (BMD) is at least 2.5 standard deviations below the young adult reference mean, expressed as a T-score. The classification based on T-scores is as follows: 1. Normal: T-score > −1.0; 2. Low bone mass: T-score between −1.0 and −2.5; 3. Osteoporosis:T-score < −2.5. If low bone mass leads to a fall or fracture, the condition may be classified as severe osteoporosis (or confirmed osteoporosis).

The single-cell data in our study were obtained from the Gene Expression Omnibus (GEO) database under accession GSE147287, which contains scRNA transcriptomic sequencing data from a 67-year-old postmenopausal osteoporosis patient’s bone marrow biopsy. CD271+ bone marrow-derived mononuclear cells (BM-MNCs) were first isolated from the tissue and then sequenced using the Illumina NovaSeq 6000 system.²⁶ The data are available at https://www.ncbi.nlm.nih.gov/geo/.

Selection of Genetic Instruments

In MR analysis, we used single nucleotide polymorphisms (SNPs) closely related to the exposure as instrumental variables (IVs). Based on previous research, we adopted a threshold of P < 1 x 10^-5 for SNP selection.^27–29 We then conducted further screening of these SNPs to eliminate those in linkage disequilibrium(LD), with the inclusion criteria set to an LD distance threshold of 10,000 kb and r² < 0.001, to ensure the independence among genetic variants.³⁰ The calculation of the F-statistic plays a vital role in Mendelian Randomization (MR) analysis, as it is used to evaluate the power of the MR analysis. F-statistic>10 suggests a strong association between the instrumental variables and the exposure.³¹ We calculate the F-statistics for each SNP individually and cumulatively, and eliminate those that do not meet the criteria to avoid the influence of weak instrument bias. Additionally, palindromic SNPs have been removed. And we also used proxy SNPs with high linkage disequilibrium, selecting those with R2 > 0.8.The Steiger test is used to verify the causal direction of each SNP, ensuring that the correct causal relationship is considered in the analysis and preventing interference from reverse causality.^32–34

MR Analysis

To test the causal effect of exposure on outcome, we conducted MR analysis using many methods, with the Inverse Variance Weighted (IVW) method as the primary one.³⁵ The IVW method is a valid analysis conducted under the fundamental premise that all genetic variations are effective instrumental variables. It has strong causal relationship detection capabilities and allows for heterogeneity among SNPs.^36–38

Additionally, on the one hand, the heterogeneity in this MR analysis is tested through Cochran’s Q statistic. If the P-value of Cochran’s Q test is less than 0.05, it indicates significant heterogeneity, suggesting that the fixed-effects model may not be appropriate.³⁷ In the presence of potential heterogeneity among SNPs, the random-effects IVW method is employed to avoid the bias caused by the fixed-effects IVW method.¹⁹ On the other hand, the MR Egger intercept test were used to evaluate the global horizontal pleiotropy of the instrumental variables.³⁵ P-values greater than 0.05 for both methods indicated no evidence of horizontal pleiotropy.^32,37,38

We also used other approaches, including the Bayesian Weighted Mendelian Randomization(BWMR), Debiased inverse−variance weighted method(DIVW), Maximum likelihood and Robust adjusted profile score (RAPS) etc. to reveal the robustness of our results.

Mediation Analysis

We performed a mediation MR analysis, with “immune cell phenotypes” as the exposure, “inflammatory cytokines” as the mediator, and “osteoporosis” as the outcome. We further conducted a two-step MR design for mediation analysis to explore whether the “Mediator” mediates the causal pathway from “Exposure” to “Outcome”. First, we used two-sample MR to assess the effect of “Exposure” on the “Mediator” (Figure 2a). Second, we used two-sample MR to evaluate the effect of mediating phenotypes, which were statistically significantly associated with “Exposure”, on the “Outcome” (Figure 2b). The total effect obtained from the previous MR analysis (Figure 2c) can be decomposed into an indirect effect (mediated by the mediator, a×b in Figure 2) and a direct effect (not mediated by the mediator, Figure 2c). By dividing the indirect effect by the total effect, we were able to determine the percentage of the total effect mediated by the mediator. Additionally, we computed the 95% confidence interval using the delta method.

Figure 2 Relationship between immune cell phenotypes and osteoporosis with inflammatory cytokines as mediators in the Mendelian randomization. Here, “c” is total effect; “ c” is direct effect of immune cell phenotypes on osteoporosis; “a” is the causal effect of immune cell phenotypes on inflammatory cytokines; “b” is the causal effect of inflammatory cytokines on osteoporosis.

All results are presented as effect sizes (ES) with corresponding 95% confidence intervals (CI). All analyses were conducted using the TwoSampleMR and MRPRESSO software packages in R version 4.3.3.

Furthermore, in this study, we conducted colocalization analyses of GWAS data corresponding to traits and diseases within exposures, mediators, and outcomes using the “coloc” R package. The aim was to assess whether they share the same genetic causal variants. This analysis is based on a Bayesian model and includes the Posterior Probability Hypothesis (PPH) of five assumptions. Conventionally, a PP4 value > 0.8 is considered the threshold for collocation significance. Given that “coloc” assumes the existence of a single causal variant, we also performed a “coloc.SuSiE” analysis, which accounts for the possibility of multiple causal variants.

Single-Cell RNA Sequencing

Processing osteoporosis scRNA-seq data using the R package Seurat.³⁹ First, low-quality cells and genes were filtered out by setting thresholds for the number of features and counts per cell (200 < nFeature < 10,000; 1000 < nCount < 100,000). Additionally, thresholds were applied to filter cells based on the proportions of mitochondrial, ribosomal, and hemoglobin genes (0 < pMT < 20, 0 < pRP < 55, 0 < pHB < 5) (Figure 3). And then, we performed Uniform Manifold Approximation and Projection (UMAP) for dimensionality reduction. To identify the optimal clustering results, we tested a range of resolution values, increasing from 0.1 to 2 in increments of 0.1. We then utilized the R package SingleR for automated cell type annotation. Subsequently, key genes were identified within specific cell types, followed by pathway enrichment analysis.

Figure 3 Single-Cell RNA Sequencing Analysis of Immune Cell Characteristics in a Postmenopausal Osteoporosis Patient. Low-quality cells and genes were filtered by setting thresholds for the number of features per cell (200 < nFeature < 10,000) and counts (1000 < nCount < 100,000). Additionally, mitochondrial gene ratio (0 < pMT < 20) and ribosomal gene ratio (0 < pRP < 55) thresholds were applied to further filter cell quality.

GSEA Analysis

This study utilized Gene Set Enrichment Analysis (GSEA) to annotate the functions of differentially expressed genes. Gene sets were obtained from the GSEA website (GSEA MSigDB) and included categories C2-C8: signaling pathway collections (eg, KEGG pathways and Reactome pathways, C2), regulatory factor target gene sets (eg, microRNA and transcription factor target genes, C3), oncogene collections (eg, oncogenic signatures and cancer-related genomic maps, C4), GO-based gene sets (including biological processes, molecular functions, and cellular components, C5), tumor-specific gene sets (C6), immune-related gene sets (eg, immune cell subtype-specific genes, C7), and cell type-specific gene sets (C8).

The GSEA analysis was conducted using the R package “clusterProfiler”. Predefined gene sets were used to assess the enrichment of ranked genes in functional pathways, with normalized enrichment scores (NES), p-values, and false discovery rates (FDR) calculated to evaluate statistical significance. The enrichment results included upregulated and downregulated pathways, representing activated or suppressed gene sets, respectively. The results were visualized through bar plots for better interpretation.

Animal Model

A cohort of 20 male C57BL/6J mice, aged 7 weeks, was procured from Jinan Pengyue Experimental Animal Breeding Co., Ltd. All animal experiments conducted in this study received approval from the Animal Ethics Committee of the Shandong Provincial Hospital Affiliated to Shandong First Medical University. A total of 20 male C57BL/6 mice (7 weeks old) were randomly assigned to 2 groups (n = 10 per group): control group, Glucocorticoid-induced osteoporosis (GIOP) group. After one week of acclimatization to the housing conditions (8 weeks old), the experiment was conducted. To induce GIOP in mice, dexamethasone (100mg/kg/day) was injected into both thigh muscles for 4 weeks. The control group was injected with normal saline in the same manner. The research team assessed health status by measuring body weight twice per week, recording food and water intake daily, and observing general indicators such as animal activity and fur condition. Mice with a body weight of 25 ± 5 g were selected for data analysis. Animals that died prematurely were excluded to prevent the collection of behavioral and histological data. After four weeks, the mice were sacrificed, and vertebral and femoral osteoporosis was evaluated using micro-CT and histomorphological analysis. The analysts were blinded to the group assignments.

Micro-CT Analysis

Bone scanning was performed using an in vivo micro-CT imaging system (Scanco, Switzerland). Briefly, micro-CT was used to scan the vertebrae and femur to reconstruct their three-dimensional structures, followed by quantitative analysis using the μCT evaluation program. Bone volume fraction (BV/TV, where BV is bone volume and TV is total volume), trabecular number (Tb.N), trabecular thickness (Tb.Th), and trabecular separation (Tb.Sp) were measured to assess bone quality.

IHC Staining

After fixation in 4% paraformaldehyde, tibial samples were embedded in paraffin and sectioned. The sections were then deparaffinized with xylene and gradually rehydrated using ethanol. To quench endogenous peroxidase activity, 3% hydrogen peroxide was applied. Subsequently, the tissue sections were blocked with 10% goat serum (Sigma-Aldrich, St. Louis, Missouri, USA) at 4°C for 1 hour before being incubated overnight with primary antibodies. The sections were then treated with biotin-labeled secondary antibodies (Proteintech, SA00004-6, 1:100), followed by staining with diaminobenzidine and counterstaining with hematoxylin. The primary antibody used was anti-IL-17C antibody (Cohesion, CQA2538, 1:200).

ELISA Analysis

Plasma from the two groups of mice was collected at 4 weeks after treatment. An ELISA kit (MEIKE, CHINA) was employed to detect the levels of IL-17C. In brief, the test plate containing 10 μL serum sample and 40 ul sample diluent was added to per well. Then, 50 μL HRP-conjugate reagent was added to each well and incubated at 37 °C for 60 min. Wash buffer was used to wash the wells five times (60 s, per wash). Then, 50 μL of substrate A solution and substrate B solution were mixed together at 37 °C for 15 min. Then, 50 μL of stop solution was added to each well. Finally, a spectrophotometer (Thermo Fisher, USA) was used to test the light absorbance.

Flow Cytometry

To identify and isolate IgD⁺CD24⁺ B cells, mouse blood was collected and anticoagulated with EDTA (Solarbio). Red blood cells were lysed using RBC lysis buffer (Solarbio), followed by centrifugation at 500g for 5 min. The supernatant was discarded, and the cell pellet was resuspended in 1 mL of PBS (Servicebio) supplemented with FBS (Gibco), then centrifuged again, and the supernatant was removed. The cell pellet was resuspended in 100 μL of PBS containing 1% BSA (Solarbio) and stained with CD24 (FITC Anti-Mouse CD24 Antibody, E-AB-F1179C, Elabscience) and IgD (PE Anti-Mouse IgD Antibody, E-AB-F1189D, Elabscience) at 4°C for 30 min. After staining, 1 mL of PBS containing 1% BSA was added, followed by centrifugation at 300g for 5 min. The cell pellet was then resuspended in 500 μL of PBS and analyzed using a flow cytometer (BD LSRFortessa).

Western Blotting Assay and Antibodies

Femur samples were soaked in liquid nitrogen for 10 minutes, then the bones were ground into powder using a mortar and pestle. The powder was collected and added to tissue lysate (RIPA Lysis buffer, Solarbio, Beijing, China) and Protease inhibitors (Phenylmethyl sulfonyl fluoride, Solarbio, Beijing, China) to extract protein samples. The protein concentration of the samples was determined using the BCA Protein Assay Kit (Solarbio, Beijing, China). The protein samples of different groups were electrophoresis on SDS-PAGE in equal amounts. Following electrophoresis, the PVDF membranes were transferred and subsequently blocked with 5% skim milk for 2 hours. The membranes were then incubated overnight at 4°C with primary antibodies. The PVDF membrane was placed in secondary antibody (Goat Anti-Rabbit IgG H&L, Abcam, ab6721, 1:5000) for 1h. TBST was washed again for three times, and Amersham Imager 680 ultra-sensitive multifunctional imager was used for visualization and detection of protein bands. Normalized with β-Actin. The following antibodies were used: anti-IL-17C (Cohesion, CQA2538, 1:1000). All experiments were conducted in triplicates.

Result

Exploration of the Causal Effect of Immune Cell Phenotypes on Osteoporosis

We employed the IVW method as the primary approach for MR analysis, conducting a two-sample MR to investigate the relationship between immune cell phenotypes and osteoporosis (Supplementary Table 1). We identified 32 immune cell phenotypes that show a causal relationship with osteoporosis, with a significance level of 0.05. Among these, 23 are risk factors for osteoporosis, while 9 serve as protective factors (Figure 4a). For instance, CD25++CD8+T cell%T cell serves as a protective factor against osteoporosis. These cells can be induced by osteoclasts, subsequently inhibiting bone resorption and forming a negative feedback loop.⁶ Results from other MR analysis methods for the immune cell phenotypes and osteoporosis are similar to those obtained by the IVW method. Cochran’s Q test and MR-Egger regression indicated no heterogeneity. MR Egger indicated no evidence of horizontal pleiotropy (all P-values > 0.05), supporting the robustness of the results. The results of the colocalization analysis indicate that the two do not share a common causal variant within the specified region (Figure 4b).

Figure 4 The causal effect of immune cell phenotypes on osteoporosis. (a)Forest plots show the causal effect of immune cell phenotypes on osteoporosis. (b) Colocalization analysis results for immune cell phenotypes on osteoporosis. Abbreviation: OR, odds ratio; CI, confidence interval.

Exploration of the Causal Effect of Inflammatory Cytokines on Osteoporosis

Subsequently, using inflammatory cytokines as the exposure and osteoporosis as the outcome, a two-sample MR analysis was conducted with IVW as the primary analytical method (Supplementary Table 2). We identified 38 inflammatory cytokines that show a causal relationship with osteoporosis, with a significance level of 0.05. Among these, 21 are risk factors for osteoporosis, while 17 serve as protective factors (Figure 5a). As a representative example, TNF-related apoptosis-inducing ligand (TRAIL) plays a crucial role in bone metabolism by promoting osteoclast formation, which subsequently disrupts bone homeostasis.⁴⁰ Results from other MR analysis methods for the inflammatory cytokines and osteoporosis are similar to those obtained by the IVW method. Cochran’s Q test and MR-Egger regression indicated no heterogeneity. MR Egger detected no horizontal pleiotropy (all P-values > 0.05), reinforcing the reliability of the findings. The results of the colocalization analysis indicate that the two do not share a common causal variant within the specified region (Figure 5b).

Figure 5 The causal effect of inflammatory cytokines on osteoporosis. (a) Forest plots show the causal effect of inflammatory cytokines on osteoporosis. (b) Colocalization analysis results for inflammatory cytokines on osteoporosis. Abbreviation: OR, odds ratio; CI, confidence interval.

Exploration of the Causal Effect of Immune Cell Phenotypes on Inflammatory Cytokines

Next, we employed the IVW method as the primary MR analysis to conduct a two-sample MR analysis examining the relationship between immune cell phenotypes and inflammatory cytokines (Supplementary Table 3). We identified 25 immune cell phenotypes that show a causal relationship with inflammatory cytokines, with a significance level of 0.05. Among these, 13 are positive factors, while 12 serve as negative factors (Figure 6a). Hodge reported that in NKT-like cells, which are predominantly CD4⁻CD8⁻, the loss of CD28 was associated with an increased population of NKT-like cells expressing granzyme B, IFN-γ, and TNF-α.⁴¹ This is consistent with our findings. Results from other MR analysis methods for the immune cell phenotypes and inflammatory cytokines are similar to those obtained by the IVW method. Cochran’s Q test and MR-Egger regression indicated no heterogeneity. MR Egger indicated an absence of horizontal pleiotropy, supporting the robustness of the findings. The results of the colocalization analysis indicate that the two do not share a common causal variant within the specified region (Figure 6b).

Figure 6 The causal effect of immune cell phenotypes on inflammatory cytokines. (a)Forest plots show the causal effect of immune cell phenotypes on inflammatory cytokines. (b) Colocalization analysis results for immune cell phenotypes on inflammatory cytokines. Abbreviation: OR, odds ratio; CI, confidence interval.

Inflammatory Cytokines Mediate the Association Between Immune Cell Phenotypes and Osteoporosis

We utilized the two-sample MR analysis method to conduct pairwise examinations of the datasets related to immunity, inflammation, and osteoporosis and conducted mediation MR analysis. Finally, we summarized the results of our mediation MR analysis and identified an mediating factors. We computed the total effect (TE), direct effect (DE), and indirect effect (IE), and ascertained the proportion of the indirect effect (IE).

IgD+ CD24+ B cell %lymphocyte demonstrated a risk role in relation to Interleukin-17C levels (β = 0.126, P =0.00015), with Interleukin-17C levels identified as a risk factor for osteoporosis (β=0.623, P =0.00118). Furthermore, IgD+ CD24+ B cell % lymphocyte emerged as a risk factor about osteoporosis. (β = 0.504, P = 0.00098). Therefore, IgD+ CD24+ B cell %lymphocyte contributes to the occurrence of osteoporosis by increasing Interleukin-17C levels. The mediated effect was 0.07837, accounting for 15.5% of the total effect, while the direct effect was 0.42610 (Figure 7).

Figure 7 Results of mediation analyses.

Single-Cell RNA-Sequencing

We analyzed the scRNA transcriptomic sequencing data from a bone marrow biopsy of a 67-year-old postmenopausal osteoporosis patient to investigate the immune characteristics of the disease, with immune cells clustered based on their unique molecular markers (Figure 8a). Subsequently, we extracted B cells and further annotated them, resulting in the identification of our target cell population: IgD+ CD24+ B cells (Figure 8b). We further screened for highly expressed genes in IgD+ CD24+ B cells and performed gene set enrichment analysis. The results indicated that the highly expressed genes in this cell population are enriched in pathways related to cancer, inflammatory immune responses, and musculoskeletal function. (Figure 8c).

Figure 8 Single-Cell RNA Sequencing Analysis of Immune Cell Characteristics in a Postmenopausal Osteoporosis Patient. (a) Immune cells were clustered based on unique molecular markers. (b) B cells were extracted and further annotated, identifying the target cell population: IgD+ CD24+ B cells. (c) Highly expressed genes in the IgD+ CD24+ B cell population were screened, followed by gene set enrichment analysis.

IgD+CD24+ B Cells and IL-17C Play Critical Roles in Osteoporosis

The classic glucocorticoid-induced osteoporosis model was established in mice.⁴² To assess the effectiveness of osteoporosis model, micro-CT was used to scan the vertebrae and femur. The results confirmed the successful establishment of the glucocorticoid-induced osteoporosis model, with a notable difference in bone mass between the osteoporosis group and the control group. (Figure 9a)

Figure 9 IgD+CD24+ B cells and IL-17C play critical roles in osteoporosis. (a)The micro CT results in osteoporosis group and control group. (b and c) The proportion of IgD+CD24+ B cells in lymphocytes of mice in osteoporosis group and control group. Horizontal axis: CD24, Vertical axis: IgD; (d) ELISA was used to measure IL-17C levels in the serum obtained from mouse eyeball blood; (e) Western blot was performed to detect IL-17C levels in mouse bone tissue. (f) Immunohistochemistry (IHC) was used to assess IL-17C levels in bone tissue.

Membrane IgD was employed to label a subset of B cells, and flow cytometry was used to measure the proportion of IgD+CD24+ B cells among lymphocytes in the blood of osteoporotic mice and control mice. The results demonstrated a significant increase in the proportion of IgD+CD24+ B cells in the blood of osteoporotic mice (Figure 9b and c). ELISA and Western blotting were used to measure IL-17C levels in blood and bone tissue, respectively. The results showed that IL-17C levels were elevated in both blood and bone tissue of osteoporotic mice (Figure 9d and e). Immunohistochemistry results showed increased levels of IL-17C in the bone tissue of osteoporotic mice (Figure 9f and g). These results indicate that IgD⁺CD24⁺ B cells and IL-17C play a crucial role in osteoporosis.

Discussion

The complicated connections and cross talk between immune cells, inflammatory factors, and bone metabolism are attracting more attention. In this study, we employed Mendelian randomization analysis to investigate the potential causal relationships between 731 immune cell phenotypes, 41 and 91 inflammatory markers, and the risk of osteoporosis. Our research utilizes multiple instrumental variables from three GWAS to enhance the statistical power of detecting causal relationships and conducts mediation analyses to provide a more accurate assessment of effect sizes. Our research results indicate that there are 23 immune cell phenotypes that serve as risk factors for osteoporosis, while 9 serve as protective factors. And 21 inflammatory cytokines have been identified as risk factors for the condition, while 17 serve as protective factors. Additionally, causal relationships were observed between 25 immune cell phenotypes and inflammatory cytokines.

Among these findings, IgD+ CD24+B cell %lymphocyte demonstrated significant risk associations in both Drug-induced osteoporosis with pathological fracture (OR=1.66, p=9.790 e-04) and Drug-induced osteoporosis (OR=1.45, p=0.046). Conversely, CD28+ CD4-CD8- T cel1%T cell exhibited protective effects against Drug-induced osteoporosis with pathological fracture (OR=0.76, p=0.008) and Drug-induced osteoporosis (OR=0.79, p=0.039). Moreover, different inflammatory factors play distinct roles in various types of osteoporosis (Figure 10). In particular, IgD+ CD24+ B cell %lymphocyte exacerbates the risk of developing osteoporosis by influencing the levels of interleukin-17C. In addition, single-cell RNA sequencing revealed various pathways involving IgD+ CD24+ B cells in osteoporosis, highlighting their roles in musculoskeletal-related pathways. These findings provide insights for a deeper understanding of this disease. Our experimental results also highlight the important role of IgD⁺CD24⁺ B cells and IL-17C in osteoporosis.

Figure 10 Different inflammatory factors play distinct roles in various types of osteoporosis.

Previous studies have laid the groundwork for us. Many studies have partially revealed the close relationship between the immune cells, inflammatory factors, and bone metabolism. For instance, the roles of neutrophils, mast cells, macrophages, TNF-α, IL-6, IFN-γ, and IL-27 have been partially elucidated. Chakravarti et al found that Neutrophils promote bone resorption by increasing the expression of mRANKL.⁴³ Buckley et al and Malone reported that mast cells trigger osteoclast generation by producing pro-inflammatory mediators such as histamine, TNF-α, and IL-6.^44,45 Zhang et al have shown that the primary role of M1 macrophages is to promote osteoclastogenesis through high levels of reactive oxygen species (ROS) and pro-osteoclastogenic cytokines such as TNF-α and IL-1β, while M2 macrophages have a bone-protective effect.⁴⁶ Yao et al said that the excessive release of TNF-α triggers a systemic hyper-inflammatory response, leading to the overactivation of osteoclasts and impaired bone formation.⁴⁷ These studies have elucidated the role of immune cells and inflammatory factors in bone degradation, consistent with our findings. In addition, Xu et al elaborated on the effects of cytokines on the remodeling of bone in osteoporosis, particularly concerning osteoblasts and osteoclasts. For instance, IL-6 can inhibit the differentiation of osteoblasts while directly and indirectly stimulating the formation of osteoclasts, and it also suppresses the differentiation of osteoclast precursors into osteoclasts. IFN-γ promotes the expression of genes related to osteoblast differentiation, thereby stimulating their development. And IL-27 can inhibit the apoptosis of osteoblasts and suppress the generation of osteoclasts.⁴⁸ These results support the protective role of cytokines in bone metabolism observed in our study, which may be attributed to the functional diversity of cytokines and the varying responses of different cells to specific cytokines.

The term IgD+ CD24+ %B cells refers to the proportion of B cell subgroups that simultaneously express IgD and CD24 within the total B cell population. IgD+ CD24+ B cells are a special subset of B lymphocytes, characterized by a unique immunophenotype. IgD is a membrane-associated immunoglobulin on B cells that interacts with antigens and plays a role in the B cell immune response.⁴⁹ CD24, a cell surface molecule, which plays a role in cell-to-cell interactions and the interaction between cells and their external microenvironment.⁵⁰ It is linked to the adhesion and migration capabilities of B cells.^51–53 Mensah et al found that, in myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) patients, the frequency and expression of CD24 on total B cells are increased compared to healthy controls (HC), and this is confined to IgD+ subsets.⁵⁴ This suggests that IgD+ CD24+ B cells also play an important role in other diseases. In addition, Fischer et al demonstrated that B lymphocytes are key regulators of osteoclast formation through the secretion of granulocyte colony-stimulating factor and the RANKL/osteoprotegerin system under estrogen-deficient conditions.⁵ Under inflammatory conditions, activated B cells secrete RANKL, which plays a crucial role in the activation of osteoclast formation.^50,55 This is consistent with our study. And CD24 is involved in the osteogenic differentiation process of bone marrow mesenchymal stem cells (BMSCs).^56,57 When BMSCs differentiate into osteogenic or adipogenic lineages, CD24 expression on the cell surface shows significant differences in expression levels.⁵⁸ Reduced CD24 expression is associated with impaired osteogenic differentiation,⁵⁸ it has been proposed as a selective biomarker for BMSC subpopulations with enhanced osteogenic potential.^50,51 This contrasts with our findings, highlighting the complexity of the biological functions of immune cells.

Our results suggests that the percentage of IgD+ CD24+ B cells in lymphocytes increases the risk of osteoporosis by affecting the levels of interleukin-17C. Interleukin-17C is a member of the IL-17 cytokine family, which consists of multiple members that may exert pleiotropic effects in host defense, autoimmunity, and the pathology of inflammation.⁵⁹ IL-17C is associated with a variety of diseases. Garcia et al stated that IL-17C is a driver of harmful inflammation during Neisseria gonorrhoeae infection of the human Fallopian tube.⁶⁰ Ramirez-Carrozzi et al said IL-17C promotes inflammation in the imiquimod-induced skin inflammation model, while it exerts protective effects in dextran sodium sulfate-induced colitis.⁶¹ Butcher et al found that IL-17C levels are elevated in psoriatic lesions, significantly affecting the abundance of F4/80+ macrophages within psoriatic plaques.⁵⁹ Our research highlights the potential hazardous role of Interleukin-17C levels in IgD+ CD24+ B cell lymphocytes, and they all contribute to an increased risk of developing osteoporosis, which has substantial implications for guiding the diagnosis and treatment strategies for osteoporosis. In addition, Goswami et al demonstrated that IL-17 signaling can prevent estrogen deficiency-induced bone loss, exerting a protective effect on bone.⁶² This is in direct contrast to the results we obtained. The contradictory phenomenon may be related to the mouse strain used in the study and the type of IL-17 examined. However, Tyagi et al have shown that increased IL-17 production induces bone loss by elevating osteoclastogenic factors (including TNF-α, IL-6, and RANKL) in osteoblasts, while functional blockade of IL-17 can prevent bone loss.⁶³ And IL-17 plays a key causal role in ovariectomy (Ovx)-induced bone loss and may be considered a potential therapeutic target in the pathogenesis of postmenopausal osteoporosis.⁶³ This is consistent with our findings. Our results are also consistent with DeSelm et al, who found that antibodies targeting IL-17 cytokines can prevent Ovx-induced bone loss.⁶⁴ These consistent findings highlight the important role of IL-17 in osteoporosis, providing a more solid theoretical foundation for future research.

In recent years, accumulating evidence has underscored the critical roles of various immune cell types in bone metabolism and demonstrated the bidirectional regulatory interactions between immune cells and osteoclasts. By identifying ligand-receptor interactions, Xu et al discovered that CD8-TEM cells and osteoclasts may crosstalk through the CD160-TNFRSF14 ligand-receptor interaction.⁶⁵ And Feng et al provided compelling evidence that various immune cells and osteoclasts participate in molecular crosstalk within the bone microenvironment through key ligand-receptor pairs.⁶⁶ Moreover, Fischer stated that during inflammatory conditions, osteoclasts have been shown to influence CD4+ T lymphocytes. Notably, the so-called “inflammatory osteoclasts”, which originate from dendritic cells rather than monocytic cells, regulate CD4+ T lymphocytes in an antigen-dependent manner and modulate their TNF-α production.⁵ We also found that the genes of IgD+ CD24+ B cells in osteoporosis were enriched in musculoskeletal-related pathways. These studies have revealed the important roles of different types of immune cells in bone metabolism and bone homeostasis and broadened the perspective for research in osteoimmunology. However, the impact of IgD+ CD24+ B cells on the bone microenvironment and the precise mechanisms underlying this effect remain unclear. It is possible that they also influence the bone microenvironment through key ligand-receptor interactions involved in molecular crosstalk or indirectly regulate bone metabolism through cytokine signaling, which requires further investigation.

Our study has several strengths. First, we used Mendelian randomization to evaluate the causal relationships among immune cell phenotypes, inflammatory factors, and osteoporosis, and through mediation analysis, we identified factors mediating osteoporosis risk, providing valuable clinical insights. Second, we employed multiple MR analysis methods, enhancing the reliability of our conclusions. Third, we incorporated single-cell sequencing analysis to further investigate the characteristics of immune cell phenotypes associated with osteoporosis. However, our study also has some limitations. Firstly, While this study utilized genetic data from publicly available GWAS databases, these data may not encompass all races and populations (primarily focusing on European populations). Therefore, whether our findings can be generalized to all populations remains to be determined. Secondly, while our analysis considers 91 and 41 inflammatory factors (which may have some overlap) as mediators, the complexity of the cytokine network suggests that there may be other mediating factors that have not been taken into account. Similarly, although we included the 731 immune cell phenotypes, the immune cell phenotypes is highly intricate, and many additional types may not have been considered. Thirdly, we did not employ bidirectional Mendelian analysis to assess potential causal relationships between the factors of interest, meaning we did not evaluate the effects of osteoporosis on immune cell phenotypes and inflammatory factors.

In summary, the relationship between the immune system and osteoporosis is complex, involving interactions among various immune cells and inflammatory factors. This MR study provides new evidence supporting the causal relationship between immunity, inflammation, and osteoporosis. The study found that IgD+ CD24+ B cell %lymphocyte may increase the risk of osteoporosis by modulating the levels of interleukin-17C. In addition, in osteoporosis, the gene expression of IgD+ CD24+ B cells is enriched in immune inflammation and musculoskeletal-related pathways. Experimental results also highlight the important role of IgD⁺CD24⁺ B cells and IL-17C in osteoporosis. Overall, these results emphasize the importance of immune cells and inflammatory factors in the risk of osteoporosis, providing new perspectives for further understanding the role of immunity in bone metabolism. Furthermore, the immune cells and inflammatory factors identified in our study provide a potential immunological basis for the diagnosis of osteoporosis and offer directions and possibilities for targeted therapy or immunotherapy.

Conclusions

Our results indicate that inflammatory cytokines play a critical role in the pathogenesis of immunity-related osteoporosis. In particular, an increase in the percentage of IgD+ CD24+ B cells in lymphocytes may enhances the risk of osteoporosis by modulating interleukin-17C levels. These results provide evidence supporting the link between immunity and osteoporosis, suggesting that targeting inflammatory pathways to intervene in immune-mediated osteoporosis may offer novel therapeutic strategies for the treatment of osteoporosis. Further clinical and experimental research is needed in future to confirm these findings.