Osteoporosis and osteoarthritis: a bi-directional Mendelian randomization study

Objective To investigate the causal relationship between low bone mineral density (BMD) and osteoarthritis (OA) using Mendelian randomization (MR) design. Methods Two-sample bi-directional MR analyses were performed using summary-level information on OA traits from UK Biobank and arcOGEN. Sensitivity analyses including MR-Egger, simple median, weighted median, MR pleiotropy residual sum, and outlier approaches were utilized in conjunction with inverse variance weighting (IVW). Gene ontology (GO) enrichment analyses and expression quantitative trait locus (eQTL) colocalization analyses were used to investigate the potential mechanism and shared genes between osteoporosis (OP) and OA. Results The IVW method revealed that genetically predicted low femoral neck BMD was significantly linked with hip (β = 0.105, 95% CI: 0.023–0.188) and knee OA (β = 0.117, 95% CI: 0.049–0.184), but not with other site-specific OA. Genetically predicted low lumber spine BMD was significantly associated with OA at any sites (β = 0.048, 95% CI: 0.011–0.085), knee OA (β = 0.101, 95% CI: 0.045–0.156), and hip OA (β = 0.150, 95% CI: 0.077–0.224). Only hip OA was significantly linked with genetically predicted reduced total bone BMD (β = 0.092, 95% CI: 0.010–0.174). In the reverse MR analyses, no evidence for a causal effect of OA on BMD was found. GO enrichment analysis and eQTL analysis illustrated that DDN and SMAD-3 were the most prominent co-located genes. Conclusions These findings suggested that OP may be causally linked to an increased risk of OA, indicating that measures to raise BMD may be effective in preventing OA. More research is required to determine the underlying processes via which OP causes OA. Supplementary Information The online version contains supplementary material available at 10.1186/s13075-023-03213-5.


Introduction
Osteoporosis (OP) is a common metabolic skeletal disease characterized by reduced bone mineral density (BMD), resulting in an increment of bone fragility among the elderly [1].It causes considerable emotional, physical, and financial stresses to patients, and often leads to disability and poor quality of life.Osteoarthritis (OA) is the most common joint disease in which degenerative changes in joint cartilage cause aseptic inflammation and involve the entire joint tissues [2].Its symptoms include joint pain, stiffness, and activity limitation.The etiology and pathological changes in OA remain largely unclear, resulting in the notable absence of curable therapies.
Does OA directly affect OP, or vice versa?Numerous observational studies have shown a strong link between OP and OA.Some have shown that patients with OA had a lower BMD, but others have reported that the progression of OA was accompanied by an increase in BMD [3].There is evidence of greater systemic BMD in OA patients at several joint locations, after adjustment for spinal osteophytes [4].A longitudinal study, however, has concluded that the development of pre-existing OA may be negatively correlated with BMD.Increased age-specific femoral neck BMD quartiles were linked to a lower risk of knee OA development throughout the course of the 8-year trial, according to research by Zhang et al. in the Framingham population [5].Additionally, Hart et al. noted a tendency for lower hip BMD in those with progressing knee OA compared to non-progressors in the Chingford Study [6].Although many studies have attempted to explore this controversy, previous observational data are limited for causal inference due to potential biases introduced by confounders.Mendelian randomization (MR), an emerging approach to evaluate causal links, has been popular in recent years.MR may reduce confounding effects and eliminate the bias of reverse causation since genotypes are independent of postnatal lifestyle and environmental variables and predate the illness process [7].
Pleiotropy may exist between OA-influencing variables such as body mass index (BMI) and genetic variables linked to BMD.Using a Mendelian randomization approach, April Hartley et al. found causal effects of hip and knee OA on BMD independent of BMI, while Liu Lin et al. found that OP reduced the incidence of OA (knee OA and hip OA).However, it is still unknown if there are site-specific causal links between low BMD and OA, and the potential signaling pathways involved are not clearly defined.Therefore, the aims of our study were to assess the causal relationship between OP and OA in the framework of a two-sample bi-directional MR study in which multiple OP measures (FN-BMD, LS-BMD, and TB-BMD) and OA measures (OA at any sites, knee OA and hip OA) will be used, and to further explore the potential signaling pathways and shared genes between OP and OA were identified.

Instrumental variables for BMD and OA
OP is defined clinically through the measurement of BMD (T score < − 2.5), which remains the single best predictor of fracture [12].In this research, BMD was used to characterize three phenotypes of OP measurements at various locations, including FN (n = 32,735), LS (n = 28,498), and systemic TB (n = 66,628).For the three BMD phenotypes, genetic differences in GWAS for poor BMD (− 2.5 T score as cutoff ) were employed as instrumental factors.Based on two recent GWAS, we selected instrumental variables for FN-BMD, LS-BMD, and TB-BMD in our main analysis as independent genetic variants at the study-specific genome-wide significance level (P = 5 × 10 −8 ) [8,9].The Genetic Factors for Osteoporosis (GEFOS) that published a meta-analysis of 32,735 FN-BMD and 28,498 LS-BMD in the European population in 2012 were included [8].The other was the largest meta-analysis to date on FN-and LS-BMD, including 17 GWAS and 32,961 individuals of European descent [13].Three BMD (per SD) phenotypes were adjusted for age, sex, weight, and height in the previous GWAS studies.Summary statistics for a GWAS meta-analysis research comprising 66,628 European participants were included in the GWAS dataset for TB-BMD [9].From the genomic-wide meta-analysis spanning the UK Biobank and arcOGEN with 384,838 Europeans, summary-level data for three OA phenotypes, including OA at any sites, knee OA, and hip OA, were collected [10].Sensitivity analysis was performed using the latest OA databases released in 2021, which conducted a GWAS meta-analysis containing 826,690 individuals (177,517 with OA).They discovered 100 risk mutations across 11 OA phenotypes that were independently linked with the illness, 52 of which were previously unrelated to the condition [11].
Genetic variations linked to low BMD were considered as instrumental SNPs in the forward stage (FN-BMD, LS-MD [8], and TB-BMD [9]).We used linkage disequilibrium (LD) [14] clustering based on r 2 > 0.001 and removed variations within 1 Mb distance from other SNPs with a greater connection to fulfill the independence across instrumental SNPs for each exposure.Additionally, we harmonized the impact of these instrumental SNPs where possible and eliminated those that were absent from the GWAS of the outcomes to make sure that all associated risk factors and result alleles were on the same strand.Then, SNPs associated with low BMD were selected (FN-BMD 47, LS-BMD 47, TB-BMD 43) (Table S1).In the reverse stage, genetic variants associated with OA were used as instrumental SNPs [10].Similarly, we used LD clumping based on r 2 > 0.001 and removed variations within 1 Mb distance from other SNPs with a strong correlation to guarantee independence across instrumental SNPs for each exposure.To ensure that all corresponding risk factors and outcome alleles were on the same strand, we harmonized the effect of these instrumental SNPs where possible, and those not present in the GWAS of the outcomes were removed.At the same time, SNPs that might have horizontal pleiotropy were removed and the frequency of efficient alleles was used to ensure that palindrome instruments were correctly aligned where possible.Finally, SNPs associated with OA were selected (OA at any sites 27, knee OA 10, hip OA 26) (Table S2).

Statistical analysis MR analysis
The causal relationship between each exposure and each outcome was assessed using the inverse variance weighted (IVW) approach with a fix effect model.We excluded instrumental variables (IVs) that were substantially linked with outcome.The IVW approach was often regarded as the most trustworthy indicator in MR analysis when there was no sign of directional pleiotropy (P for MR-Egger intercept > 0.05) [15].The IVW method used the log (OR)/β coefficient of the disease SNP divided by the log (OR)/β coefficient of the exposed SNP to obtain the Wald estimates for each SNP, and then combined these Wald estimates using a method similar to metaanalysis [16].When each genetic variation satisfied the IV hypothesis, the IVW method could provide a consistent estimate of the causal effect of exposure on the outcome.Cochran's Q statistics were used to evaluate the IV heterogeneity.In order to further confirm the reliability of MR estimations, we used the MR-Pleiotropy Residual Sum and Outlier techniques (MR-PRESSO), which can identify and eliminate probable pleiotropic IVs and offer the outlier-adjusted estimates, to the IVs for IVW analysis (P-value < 0.05).In addition, we used weighted median, simple median, and MR-Egger as complementary analysis methods to test the robustness of the IVW method using random-effect model estimation.Utilizing aggregate data, the weighted median estimate was produced to give protection against ineffective instruments and reliable estimations of causation if at least 50% of the weight originates from IVs [17].

Subgroup analysis of age
Considering IVs of TB-BMD are from GWAS analyses of all age groups which may lead to a decrease in the accuracy of MR analyses [9], we did a detailed analysis by age groups.We used the IVW method of the random-effects model to evaluate the causal effects of each exposure and each outcome.TB-BMD data on the basis of age was divided into four stages, which were 15 to 30 (n = 4,180), 30 to 45 (n = 10,062), 45 to 60 years old (n = 18,805), and over 60 years old (n = 22,504), respectively [9].MR analysis method was used to estimate the correlation between OA and the corresponding phase TB-BMD.

Gene ontology enrichment analysis
The gene ontology (GO) [18] which involves biological process, molecular function, and cellular component were analyzed to further validate whether the potential targets were indeed related to OA.We connected the causative BMD lead SNPs identified in several OA symptoms to the adjacent genes.After using the Benjamini-Hochberg method for adjustment [19], P < 0.05 was used as the significance cutoff in our study.

Expression quantitative trait locus analysis
We tested the relationships between OA and OP features using a gene-based strategy and extrapolated GTEx V746 expression levels for the musculoskeletal tissues.The elastic net model is used by MetaXcan [20] to impute the cis-genetic component of expression into a much bigger set for a collection of reference participants for whom both gene expression and genetic variation have been assessed.It then did a transcriptome-wide association study to identify meaningful expression-trait connections and corroborated the imputed gene expression to the trait of interest.With 20,000 genes spread over 48 tissues, we employed a conservative Bonferroni correction to account for the gene-tissue pairings, which resulted in a significance threshold of 5.20 × 10 −8 .We calculated the likelihood of each GWAS and expression quantitative trait locus (eQTL) signal colocalizing in each significant MetaXcan results using Coloc47 to lessen the impact of LD confounding on the MetaXcan results when various causal SNPs were affecting expression levels and phenotypic traits in a GWAS [21].R 3.5.3and STATA 14.0 were used to conduct all statistical analyses.Unless otherwise stated, statistical significance was defined as P < 0.05.Because the objective of our study was the causal relationship between different types of BMD and OA, we did not apply a multiple testing correction to the reported P. GO enrichment analysis was performed by websites tool "FUMA" at http:// fuma.ctglab.nl.

Stage 1: forward MR analysis of the effects of BMD on OA
By applying a two-sample MR approach, we first investigated the causal associations of low BMD on OA.We identified 47 (FN-BMD), 47 (LS-BMD), and 43 (TB-BMD) IVs from GWAS that met genome-wide significance level (P < 5 × 10 −8 ) based on the elimination of certain missing data (r 2 < 0.001).The heterogeneity test showed no significant heterogeneity among selected IVs (Q _ P > 0.05, Table 1), indicating our MR results were not biased by heterogeneity or horizontal pleiotropy.Table S1 shows the effect sizes of a few IVs.

Stage 2: reverse MR analysis of the effects of OA on BMD
Using OA-associated SNPs as IVs, we further investigated the causative relationships between OA and low BMD in the reversal direction.We obtained 26 (OA at any sites), 10 (knee OA), and 26 (hip OA) IVs without effects of LDindependent (r 2 < 0.001), reaching P < 1 × 10 −5 from GWAS.When examining the causal relationship of OA at any sites on FN-BMD and TB-BMD, we deleted one of the SNPs (rs528981060) because of its significant heterogeneity on outcomes.The heterogeneity test showed no significant heterogeneity (Q _ P > 0.05, Table 2) in selected IVs of OA on low BMD, demonstrating that neither horizontal pleiotropy nor heterogeneity affected our MR findings.The negative control analysis showed that OA was not associated with selected IVs, suggesting that the selected IVs for exposures in this study were appropriate (Table S3).The IVW method showed that genetically predicted three OA exposures were not associated with three BMD outcomes (β ranged from 0.018 to 0.061).The estimates from MR-Egger, simple median, and weighted median were largely similar (β ranged from 0.025 to 0.181) (Table 2).

Subgroup analysis by age
Subgroup analyses by age were showed in Table S4.In those with 30-45 years of age, TB-BMD only had a strong causal effect on OA at any sites (P = 0.017) and hip OA (P = 0.021).TB-BMD had a strong causal effect on all the OA phenotypes in participants between 45 and 60 years of age (OA at any sites, P = 0.001; knee OA, P = 0.001; hip OA, P = 0.009).In those above 60 years of age, TB-BMD also had a strong causal effect on OA at any sites, knee OA and hip OA.However, we did not find valid data for MR estimates (data not shown) in participants below 30 years of age.Due to the absence of FN-BMD and LS-BMD subgroup data in different age groups, we were unable to complete the relevant MR analysis.

GO enrichment analysis
The findings of the GO enrichment study are shown in Fig. 2. BMD on different OA sites revealed high enrichment in a number of important regulatory pathways according to GO enrichment analyses.For knee OA, 10 GO biological processes were observed to be involved, such as ossification, cell signaling by WNT, and bone mineralization (Fig. 2A).For hip OA, 68 GO biological processes were observed to be involved and the most significant ones were presented in Fig. 2, such as cell signaling by WNT, cell surface receptor signaling pathway, ossification, skeletal system development and neurogenesis (Fig. 2B).For OA at any sites, 15 GO biological processes were observed to be involved, such as regulation of gene expression, RNA biosynthetic process, tissue morphogenesis and osteoblast (Fig. 2C).

eQTL co-location analysis
The eQTL colocation analysis showed that four BMD susceptibility genes (ANAPC1, GALNT3, LIN7C, and DDN) had variants co-located on musculoskeletal tissues (Table 3).We found that gene GALNT3 was widely presented in all three types of BMD, whose variants co-located on LS-BMD were the most prominent (rs1346004, rs11680288).Among all susceptibility genes with variants co-located discovered, the effect of gene DDN was the largest (P < 0.001).We also found four hip OA susceptibility genes (LTBP3, MLXIP, SMAD3, and MAPT) that had variants co-located on musculoskeletal tissues.For hip OA co-located variants, the effect of gene SMAD3 was the largest (P < 0.001), followed by MAPT, LTBP3, and MLXIP (all P < 0.001).No genes had variants co-located on knee OA or OA at any sites.

Discussion
Using the bi-directional two-sample MR method, we found evidence for a causal association of low BMD on OA, and the association was more prominent in people over 45 years old.No clear evidence of a causal relationship from OA to BMD was found.The eQTL co-location analysis revealed that four BMD susceptibility genes (ANAPC1, GALNT3, LIN7C, and DDN) and four hip OA susceptibility genes (LTBP3, MLXIP, SMAD3, and MAPT) had variants co-located on musculoskeletal tissues, suggesting these genes may play roles in regulating the development of OP and OA.
Previous studies reported that OP may be one of the etiologies to promote OA development.A meta-analysis showed that the worldwide incidence of OP increased with aging, and could promote OA through a number of mechanisms [22].For example, a prospective study showed older participants with radiographic hip and knee OA had higher total hip bone loss over 2.6 years [23].Developing OP has been shown to aggravate cartilage lesions in an experimental model of OA in rabbits [24].Similar findings were made by Bellido M. et al., who discovered that the fragility and poor quality of subchondral bone may also increase cartilage injury in addition to promoting bone remodeling [25].An increased rate of bone turnover has been shown in animal models of early OA [26].In cadaver specimens of people with early OA, the elastic modulus of trabecular subchondral bone from proximal tibiae decreased despite an increase in bone volume [27].This local bone softening may accompany a drop in BMD, which may be the result of inadequate mineralization brought on by increased bone remodeling [27].In summary, our findings are generally consistent with the above observational studies that low BMD may be associated with an increased risk of OA.In contrast, there were some observational studies showing that higher femoral neck and total body BMD were associated with an increased risk of incident OA [28].A cross-sectional and longitudinal study reported a positive correlation between subchondral BMD and OA severity in patients [29].In general, these studies indicated the detrimental effect of higher BMD on OA development.
Regarding the possible causal relationship from OA to OP, results varied among studies.The MR Study by Hartley et al. suggested that a reduction in BMD would increase the incidence of OA, while Lin et al. suggested that a reduction in bone mineral density had a protective effect on OA [30,31].The possible explanations for the discrepancies include (1) different sources of exposure and differences in the final selection of SNPs, (2) although we tried to minimize phenotypic heterogeneity, it could not disappear, (3) different methods for assessing OA and/or OP [27], and (4) lack of adjustment for essential confounders such as age, gender, weight, and/or BMI [32].The protective effect of high BMD on OA may be through its effect on reducing the risk of joint space loss [5], and the mechanism of OP causing OA could be the bone loss in the subchondral bone resulting in the collapse of the articular surface and leading to uneven stress on the articular cartilage, which leads to secondary osteophyte proliferation and cartilage damage [33].It is worth noting that Lin et al. [31] only analyzed the causal relationship from OP to OA, and the summary data of exposure and outcome in their study were both from UK   [30], indicating that improving BMD might help to prevent OA from a genetic standpoint.Kindly noted the previously study assessed BMD by heel ultrasound as the only exposure, while our study used dual-energy X-ray assessed FN-BMD, LS-BMD, and TB-BMD as the exposures.
In our study, we found that the causal effect of FN-BMD on knee OA was stronger than that of hip OA, while LS-BMD had a strong causal effect on OA at any sites, knee OA and hip OA.Hackinger et al. [34] discovered a genetic link between LS-BMD (but not FN) and OA, which is consistent with common biological pathways causing both BMD and OA.On the other hand, TB-BMD only had a weak causal effect on hip OA, possibly because of the special trabecular structure of the hip [35].Moreover, although we used two databases for verification, we cannot rule out the possibility of false positives.In the reverse MR analyses, no evidence for a causal effect of OA on BMD was found.This could be due to the influence of the genetic variation on the result is not entirely mediated by the binary exposure (OA definition), power estimations are probably conservative in this case.However, estimating approaches for a binary exposure call for strong assumptions that are unlikely to be physiologically tenable in typical MR situations, making it difficult to perform tests for causal effects without employing exposure information [36].
GO enrichment analysis found that a large number of GO biologic processes play key roles in the potential causal relationship from low BMD to OA, such as ossification [37], cell surface receptor signaling pathway [38], and cell signaling by WNT [39].Using OP mice mode, Wu et al. discovered that nine cellular elements, including the cytosol, nucleus, cytoplasm, neuronal cell body, protein complex, caveola, and endoplasmic reticulum, were involved in the anti-OP effects on OA [37].In addition, our GO analysis revealed that the primary ways in which OA affected OP were via inflammatory response, aging, responses to estradiol, glucocorticoids, and hydrogen peroxide [37].The study by Guan et al. suggested that XianLing GuBao Capsule (XLGB) might treat OP through the PI3K/AKT and MAPK pathways [38].Similarly, Xiao et al. found that PI3K/AKT/NF-κB and MAPK pathways would be the potential mechanism for XLGB in the treatment of OA [39].Additionally, aberrant bone remodeling and reduced mineralization have been linked to abnormal OPG and RANKL expression in OA osteoblasts [40].WNT protein is considered as a risk factor for the onset and development of OA.There are two primary types of WNT signaling pathways: canonical and non-canonical.Of these, the canonical WNT signaling pathway promotes osteogenesis.Osteocytes generate the inhibitor of this process, sclerostin, which prevents osteogenesis [41].Since we substituted BMD SNPs for OP, it is reasonable to assume that OP also has these biological processes and WNT signaling could be potential pathways to be targeted for treating both diseases.
Through eQTL colocation analysis [21], 3 SNPs (rs1346004, rs11680288, and rs7586085) were co-located at gene GALNT3.An observational study found that the GALNT3 gene significantly correlated with OP and the low expression of the GALNT3 gene can promote the occurrence and deterioration of OP [42].Our data showed that three newly discovered genes had eQTL effects with OP (ANAPC1, LIN7C, and DDN).In addition, previous reports deemed that the LIN7C gene significantly correlates with knee OA [43].Our study also found four genes with eQTL effect on hip OA (LTBP3, MLXIP, SMAD3, and MAPT).The LTBP3 gene may be involved in the occurrence and development of OA by regulating the TGF-β signaling pathway [44].SMAD3 gene is the target gene of circ-RNA (Circ0083429), which regulates the mRNA level of SMAD3 through the sponging of miRNA-346, thus affecting OA.Studies have found that the SMAD3 gene plays an important role in cartilage bone reconstruction and maintenance, and have confirmed the differential expression of SMAD3 in intact and degraded knee and hip cartilage [34].Top-level differentially expressed genes in OA bone, such as those in the WNT signaling pathway (TWIST1, IBSP, S100A4, MMP25, RUNX2, and CD14) and TGF-β/SMAD3 signaling pathway, are known to function in osteoblasts, osteocytes, and osteoclasts (ADAMTS4, ADM, MEPE, GADD45B, COL4A1, and FST) [45].MLXIP and MAPT have not been reported to be related to the development of hip OA.Therefore, these newly discovered genes need further confirmation.
The strengths of the current study lie on: First, to minimize possible false-positive results, summary data from two different organizations were used for outcomes, and the latest summary data of OA were also used for verification [10,11], proving the credibility of our conclusions.Second, four methods of MR analysis on 18 groups of BMD GWAS and OA GWAS summary data were performed, and subgroup analyses by age were performed for TB-BMD.Third, we used GO Enrichment analysis and eQTL to search for potential enrichment genes and to determine whether there are possible co-acting pathways.
Limitations of our study include: First, the majority of the GWAS summary data stem from people of European origin, therefore, we should be cautious when applying our findings to people with other racial and ethnic backgrounds.Second, we evaluated a linear correlation between OP and OA in our MR analysis, but we did not account for the potential of different shapes of association [15].Third, because summary statistics rather than raw data were used in the analysis, it was not possible to perform subgroup analyses, such as stratified by sex or ethnicity.Last, the results of GO and eQTL were based only on bioinformatic analyses, in vitro and in vivo experiments are needed in the future to validate our findings.

Conclusions
These findings suggested that OP may be causally linked to an increased risk of OA, indicating that measures to raise BMD may be effective in preventing OA.More research is required to determine the underlying processes via which OP causes OA.

Fig. 1
Fig. 1 Workflow of bi-directional MR analysis.A, the fundamental idea of MR analysis: If we cannot randomize the exposure, we can find a randomized instrumental variable to disentangle.① Instrumental variables were highly correlated with exposure factors.② Instrumental variables were independent of the outcome.③ Instrumental variables were not correlated with confounding factors.B, Workflow of our bi-directional MR analysis.MR: Mendelian randomization; BMD: bone mineral density; OA: osteoarthritis; FN: femoral neck; LS: lumbar spine; TB: total body

Fig. 2
Fig. 2 Gene Ontology enrichment analysis of the casual effect of Bone Mineral Density on Osteoarthritis.A, Gene Ontology enrichment analysis of the casual effect of Bone Mineral Density on knee Osteoarthritis.B, Gene Ontology enrichment analysis of the casual effect of Bone Mineral Density on Hip Osteoarthritis.C, Gene Ontology enrichment analysis of the casual effect of Bone Mineral Density on Osteoarthritis at any sites

Table 1
Mendelian randomization estimates for bone mineral density on osteoarthritisBMD bone mineral density, OA osteoarthritis, FN femoral neck, LS lumbar spine, TB total body, IVs instrumental variables, IVW inverse variance weighted

Table 2
Mendelian randomization estimates for osteoarthritis on bone mineral density BMD bone mineral density, OA osteoarthritis, FN femoral neck, LS lumbar spine, TB total body, IVs instrumental variables, IVW inverse variance weighted a IVs are 26, SNP (rs528981060) overlaps the outcome

Table 3
Significant eQTL signals for osteoporosis and osteoarthritis from musculoskeletal tissues based on the GTEx data eQTL expression quantitative trait locus, BMD bone mineral density, OA osteoarthritis, CHR chromosome, POS position, P-value: the statistical analysis of enrichment score effect size, which is used to characterize the credibility of enrichment results Biobank, which may cause bias in the results.Based on two separate OA summary datasets, our findings largely agreed with Hartley et al. 's