Rheumatoid arthritis and osteoporosis: a bi-directional Mendelian randomization study
Abstract
Many observation studies have demonstrated a close relationship between rheumatoid arthritis (RA) and osteoporosis (OP). However, the causal genetic correlation between RA and OP remains unclear. In this study, we performed bi-directional Mendelian randomization (MR) analyses to explore causal inference between these two traits. The instrumental variables for RA were selected from a large-scale genome-wide association study (GWAS) (1,523 cases and 461,487 controls). Bone mineral density (BMD) at five different sites (heel (n=265,627), forearm (FA) (n=8,143), femoral neck (FN) (n=32,735), lumbar spine (LS) (n=28,498), and total body (n=28,498)) were used as phenotypes for OP. The inverse variance weighted (IVW) method did not detect any causal effect of BMDs on RA except heel BMD (beta = -7.57 × 10-4, p = 0.02). However, other methods (MR-Egger, weighted median, weighted mode, MR-PRESSO, and MR-RAPS) showed no causal association between heel BMD and RA. Likewise, we did not find a causal effect of RA on BMD at any sites. In conclusion, we found no evidence that RA is causally associated with OP/BMD, or vice versa. We suggested that the associations found in previous observational studies between RA and OP/BMD are possibly related to secondary effects such as antirheumatic treatment and reduced physical activity.
Introduction
Rheumatoid arthritis (RA) is a kind of common autoimmune disease with hyperplasia of joint tissue and synovial inflammation that can finally lead to some serious systematic disorders, such as cardiovascular, pulmonary, skeletal disorders, and psychological [1]. It affects around 1.3 million people in the USA [2], and 0.32%-0.36% population in China [3]. One of the most severe comorbidities of RA is osteoporosis (OP), which is a chronic metabolic skeletal disease leading to an increased risk of low trauma fracture [4]. Osteoporosis can be characterized by microarchitectural deterioration of bone tissue and low bone mass. Epidemiology studies indicate that about 60-80% of RA patients have a comorbidity of OP [5]. The most commonly used measurement for OP is bone mineral density (BMD) [6].
Both RA and OP have a strong genetic component. Previous studies have suggested that the heritability of RA is approximately 60% [7], while the heritability of OP is up to 50-85% [8]. To date, genome-wide association studies (GWASs) have successfully identified more than 200 single-nucleotide polymorphisms (SNPs) for OP and explained approximately 5% of the genetic heritability [9]. GWASs for RA have identified more than 100 susceptibility loci and explained approximately 12% of the genetic heritability [10]. These two kinds of complex diseases may share some common genetic mechanisms and biological processes. For example, proinflammatory cytokines including TNF-α, IL-17, IL-6, and IL-1 have been reported to be closely associated with OP [11], and they also play important roles in the development of RA [12].
In addition, many observation studies demonstrate a close relationship between OP and RA with strong evidence. Focal or generalized bone involvement occurs frequently in RA patients. Güler-Yüksel M et al. reported that osteoporosis was found in the spine and hip in 11% and 25% of 381 recently diagnosed active RA patients [13]. Synovial membrane inflammation will lead to periarticular cortical bone loss and marginal bone erosion, while generalized osteoporosis involving the appendicular and axial skeletons probably occurs before the onset of articular disease [14]. Follow-up studies conducted by the same group demonstrated that joint damage and joint damage progression were independently associated with high bone mineral density (BMD) loss both in hip/spine and hands after 1 year of treatment [15]. Longer duration and severity of RA were also indicated as independent risk factors for vertebral fractures [16].
Does RA have a direct effect on OP, or vice versa? Due to the potential bias introduced by confounding factors, these prior observation data were limited for causal inference. The gold standard method for identifying causality is randomized controlled trial (RCT) design, but it consumes considerable time and money. In recent years, Mendelian randomization (MR) has been widely performed as an alternative method to assess causal relationships in observational data [17]. Genetic variants are used as instrumental variables (IVs) in the MR method to leverage the random assortment of genetic variants during gamete formation. If we assume that there are no genetic mating restrictions upon the population (panmixia), then the genotype distribution of this population should be unrelated to the confounding factors that typically affect observational epidemiology studies. In this regard, MR can be thought of as a “natural” RCT. Moreover, the assignment of genotypes is not affected by age, sex, lifestyle, or other environmental factors. Hence, the greatest advantage of MR is avoiding the effect of potential confounders or reverse causality compared with the general RCT design [18]. A genetic variant can be considered as an instrumental variable for a given exposure if it satisfies the 3 assumptions: 1) They are strongly associated with exposure. 2) They are independent of any known confounders. 3) They are conditionally independent of given exposure, outcome, and potential confounders, meaning that it does not affect the outcome except via the exposure, and it is not associated with the outcome due to confounding (Figure 1A).
Figure 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 (B) Workflow of our bi-directional MR analysis. MR: Mendelian randomization; BMD: Bone mineral density; RA: Rheumatoid arthritis.
Two-sample MR uses GWAS summary statistical data of both exposure and outcome traits to infer the causal association between exposure and outcome. Hence, it is not necessary to obtain the effect of the instrumental variable-exposure/-outcome association from the same sample of participants. In other words, two-sample MR allows us to perform MR between two traits using only independent GWAS summary statistics. In addition, there are some advantages to obtain summary statistical data from two different groups of participants. For example, the “winners’ curse” is unlikely to happen in two-sample MR, while it can underestimate true causal effects in one-sample MR [19]. Likewise, the weak instrument bias that biases effects towards the confounded multivariable regression result has a great impact in one-sample MR, but it is towards null in two-sample MR. The main advantage of using summary data from large GWASs in two-sample MR is the increasing of statistical power, particularly in testing effects on binary disease outcomes. Regardless of many successes of MR in investigating the potential causality of environmental factors-to-diseases or diseases-to-diseases, whether there is a causal effect between OP and RA is still unclear, though extensive evidence from observational studies showed a strong correlation between these two diseases. Therefore, we performed bi-directional MR analyses for causal inference between RA and OP in the present study. Since BMD is the most often used measurement for diagnosing OP [20, 21], summarized GWAS data for five different BMD measurements were included: the 1st~3rd BMD measurements were DXA measuring-based BMD at forearm (FA), low spine (LS), and femoral neck (FN), which were most commonly used for diagnosing of OP and OP-related fracture [22]; the 4th BMD measurement is DXA measuring-based BMD of total body, which has been reported with a strong correlation with BMD at low spine (LS) and femoral neck (FN), and with the largest collection of DXA measuring-based BMD [23]; the 5th BMD measurement is heel BMD measured by quantitative heel ultrasounds from UK biobank [24]. Although the heel BMD estimated from quantitative heel ultrasounds is not as standard as DXA-based measurements, it is also strongly associated with DXA-based BMD. In addition, benefit from its convenience and low cost, the study of eBMD in UK biobank is much larger than any DXA-based study.
The bi-directional MR study design is shown in Figure 1B. In the first stage, we examined whether BMD measurements have causal effects on RA. In the second stage, we detected whether RA is causally associated with BMD measurements.
Results
Stage1: Influence of BMD traits on RA
Respectively, we obtained 342, 3, 14, 14 and 50 IVs without effects of linkage disequilibrium (LD, r2 < 0.001) that reached genome-wide significance level (p < 5 × 10-8) from GWASs for heel, FA-, FN-, LS- and TB-BMD (Supplementary Table 1). The heterogeneity test showed no significant heterogeneity among selected IVs (Q_p value > 0.05, Table 1), except IVs of FN-BMD (Q_IVW = 24.73, Q_pval_IVW=0.02). All tests for MR Egger regression and leave-one-out analysis were negative (p for MR-Egger intercept > 0.05) (Table 1 and Supplementary Figures 1–5), indicating that our MR results were not biased by heterogeneity or horizontal pleiotropy. To demonstrate the power of selected IVs, we presented F statistics. The values of F statistics for selected IVs and the variance explained by them for heel BMD, FA-BMD, FN-BMD, LS-BMD, and TB-BMD were 97.05, 58.20, 52.33, 58.39, and 66.32, respectively. All of the F statistics were larger than 10, demonstrating that the selected IVs were strong enough to decrease any potential bias of the causal analyses. As expected, the negative control analyses presented that heel BMD, FA-BMD, FN-BMD, LS-BMD, and TB-BMD were not associated with myopia, suggesting that IVs of exposures we selected for this study were appropriate (Supplementary Tables 2, 3).
Table 1. Mendelian randomization estimates for BMD on RA.
| Exposure | Outcome | No. of IVs | Heterogeneity tests | Directional horizontal pleiotropy test | MR results |
| Methods | Cochran's Q (p) | MR-Egger intercept (p) | Method | Beta | P |
| Heel BMD | RA | 342 | MR Egger | 370.51 (0.12) | -1.87e-05 (0.20) | Inverse variance weighted | -7.57E-04 | 0.02 |
| | | Inverse variance weighted | 372.28 (0.12) | | MR Egger | -1.48E-05 | 0.98 |
| | | | | | Weighted median | -2.48E-04 | 0.62 |
| | | | | | Weighted mode | 6.49E-05 | 0.93 |
| | | | | | MR-Presso | -5.44E-04 | 0.08 |
| | | | | | MR-Raps | -5.49E-04 | 0.07 |
| FA-BMD | RA | 3 | MR Egger | 2.94 (0.09) | -2.36e-04 (0.65) | Inverse variance weighted | -5.64E-04 | 0.49 |
| | | Inverse variance weighted | 4.08 (0.13) | | MR Egger | 1.16E-03 | 0.76 |
| | | | | | Weighted median | -1.87E-04 | 0.78 |
| | | | | | Weighted mode | -7.27E-05 | 0.93 |
| FN-BMD | RA | 14 | MR Egger | 18.26 (0.11) | -5.13e-04 (0.06) | Inverse variance weighted | -3.03E-04 | 0.70 |
| | | Inverse variance weighted | 24.74 (0.02) | | MR Egger | -9.08E-03 | 0.06 |
| | | | | | Weighted median | -3.40E-04 | 0.67 |
| | | | | | Weighted mode | 1.17E-04 | 0.93 |
| LS-BMD | RA | 14 | MR Egger | 19.92 (0.07) | -2.36e-04 (0.41) | Inverse variance weighted | -1.16E-04 | 0.86 |
| | | Inverse variance weighted | 21.15 (0.07) | | MR Egger | -3.66E-03 | 0.40 |
| | | | | | Weighted median | -3.24E-04 | 0.65 |
| | | | | | Weighted mode | -5.09E-04 | 0.65 |
| TB-BMD | RA | 50 | MR Egger | 50.46 (0.38) | 9.29e-05 (0.09) | Inverse variance weighted | -3.49E-04 | 0.35 |
| | | Inverse variance weighted | 53.64 (0.30) | | MR Egger | -2.08E-03 | 0.06 |
| | | | | | Weighted median | -6.23E-04 | 0.25 |
| | | | | | Weighted mode | -6.54E-04 | 0.38 |
| MR-PRESSO and MR-RAPS were performed for pairs of exposure-outcomes with more than 100 IVs. BMD: bone mineral density; RA: rheumatoid arthritis; FA: forearm; FN: femoral neck; LS: lumbar spine; TB: total body; IVs: instrumental variables; MR-PRESSO: Mendelian Randomization Pleiotropy RESidual Sum and Outlier; MR-RAPS: Mendelian Randomization Robust Adjusted Profile Score. |
The inverse variance weighted (IVW) method supported a causative association between heel BMD and RA (beta = -7.57 × 10-4, p = 0.02), but MR Pleiotropy RESidual Sum and Outlier (MR-PRESSO) did not detect any potential pleiotropic IVs for BMD, and the corrected MR causal association between BMD and RA was negative (beta = -5.44 × 10-4, p = 0.08) (Table 1). As there were 342 IVs for heel BMD, we performed MR Robust Adjusted Profile Score (MR-RAPS) to test whether BMD affects RA through many weak instruments (beta = -5.49 × 10-4, p = 0.07) (Table 1). MR-RAPS did not show that heel BMD has a causal effect on RA. Furthermore, IVW analysis showed that FA-, FN-, LS- and TB-BMD were all negatively associated with RA (beta range from -1.16 × 10-4 to -5.64 × 10-4). In addition, MR-Egger, weighted median, and weighted mode methods did not identify any causal effect of heel, FA-, FN-, LS- and TB-BMD on RA (Table 1). Combined with results from different MR methods, we concluded that BMD has no causal effect on RA.
Stage2: Influence of RA on BMD traits
In total, we obtained 6 LD-independent (r2 < 0.001) instrumental variables (IVs) with p < 1 × 10-5 from GWAS for RA (Supplementary Table 4). The heterogeneity test showed significant heterogeneity (Q_p value < 0.05, Table 2) in selected IVs of RA on heel BMD and FN-BMD. All tests for MR Egger regression and leave-one-out analysis were negative (p for MR-Egger intercept > 0.05) (Table 2 and Supplementary Figures 6–10), indicating that our MR results were not biased by heterogeneity or horizontal pleiotropy. The value of F statistics for selected IVs is 71.71 (larger than 10), demonstrating that the IVs we selected in this study were powerful enough. As expected, the negative control analyses presented that RA was not associated with myopia, suggesting that the IVs of exposures we selected in this study were appropriate (Supplementary Tables 2, 3).
Table 2. Mendelian randomization estimates for RA on BMD.
| Exposure | Outcome | No. of IVs | Heterogeneity tests | Directional horizontal pleiotropy test | MR results |
| Methods | Cochran's Q (p) | MR-Egger intercept (p) | Method | Beta | P |
| RA | Heel BMD | 4 | MR Egger | 14.87 (0.001) | 3.00e-03 (0.63) | Inverse variance weighted | -2.52 | 0.20 |
| | | | | MR Egger | -4.33 | 0.38 |
| | Inverse variance weighted | 17.23 (0.001) | | Weighted median | -2.96 | 0.00 |
| | | | | Weighted mode | -3.26 | 0.04 |
| | MR Egger | 5.13 (0.16) | 4.00e-02 (0.61) | Inverse variance weighted | 3.93 | 0.80 |
| RA | FA-BMD | 5 | | | | MR Egger | -62.00 | 0.63 |
| | | Inverse variance weighted | 5.67 (0.22) | | Weighted median | 17.68 | 0.31 |
| | | | | | Weighted mode | 23.93 | 0.39 |
| | | MR Egger | 11.82 (0.02) | 8.00e-03 (0.58) | Inverse variance weighted | -4.66 | 0.50 |
| RA | FN-BMD | 6 | | | | MR Egger | -11.87 | 0.45 |
| | | Inverse variance weighted | 12.90 (0.02) | | Weighted median | -7.91 | 0.13 |
| | | | | | Weighted mode | -7.78 | 0.20 |
| | | MR Egger | 3.53 (0.47) | -6.00e-03 (0.53) | Inverse variance weighted | -4.96 | 0.32 |
| RA | LS-BMD | 6 | | | | MR Egger | 0.49 | 0.96 |
| | | Inverse variance weighted | 4.00 (0.55) | | Weighted median | -3.20 | 0.58 |
| | | | | | Weighted mode | -3.19 | 0.66 |
| | | MR Egger | 4.75 (0.31) | 6.34e-05 (0.99) | Inverse variance weighted | -1.74 | 0.56 |
| RA | TB-BMD | 6 | | | | MR Egger | -1.79 | 0.77 |
| | | Inverse variance weighted | 4.75 (0.45) | | Weighted median | -2.08 | 0.53 |
| | | | | | Weighted mode | -2.08 | 0.58 |
| BMD: bone mineral density; RA: rheumatoid arthritis; FA: forearm; FN: femoral neck; LS: lumbar spine; TB: total body; IVs: instrumental variables. |
Through IVW analyses, we did not detect any evidence for a causal effect of RA on BMD at any site (beta range from -4.96 to 3.93) (Table 2). The estimates from MR-Egger were consistent with these results (beta range from -62.00 to -2.18). The weighted median (WME) and weighted mode (WMO) analysis detected a significant causal effect of RA on heel BMD (p_WME = 9.93 × 10-4, p_WMO = 0.04). When there is absent evidence of directional pleiotropy (p for MR-Egger intercept > 0.05), the IVW method is considered the most reliable indicator in MR analyses [25, 26]. Therefore, we concluded that the causal association between RA and BMD is negative. For the other BMD traits, the results from weighted median and weighted mode analysis were consistent with IVW (Table 2). In summary, combined with the results from different MR methods, we concluded that there is no causal effect of RA on BMD.
Discussion
In the present study, we used bi-direction Mendelian randomization to figure out whether RA is causally associated with OP, or the other way around. Despite using the largest available public GWAS meta-analyses data, we were unable to demonstrate an association between genetic instruments for RA and BMD, or vice versa. Thus, there was no evidence for a causal relationship between RA and OP according to our MR analyses.
Previous observation studies have shown powerful evidence of an association between active RA and low BMD [27–30]. According to a population-based study, the prevalence of OP had a twofold increase in both male and female RA patients compared with healthy subjects [31]. A Korean cohort including 47,034 RA patients and 235,170 controls also indicated an increased risk of osteoporotic fractures for RA patients across all sexes, age groups, and various anatomic sites, compared with non-RA patients [32]. High disease activity, long disease duration, and joint damage were reported as determinants of reduced BMD in RA patients [15, 33]. Studies on the molecular mechanism also suggested that the pathogenesis of generalized BMD loss and focal erosions had common pathways mediated by osteoclasts, particularly by the receptor activator of nuclear factor kappa B ligand (RANKL) pathway [14, 34, 35].
However, importantly, no conclusions can be drawn on whether RA has a direct influence on OP, or the other way around. Some studies found that focal or generalized bone loss occurred before the diagnosis of RA in some patients, and BMD seems to be predominantly related to demographic factors in those patients without disease-modifying antirheumatic drugs (DMARDs) or corticosteroid treatment [14, 15]. These findings suggested that the reduction of BMD in RA patients might be partially intermediated by other factors, such as antirheumatic treatment. Our MR analysis results found that there is no causal association between RA and OP/BMD, which also suggested a secondary effect of RA on OP/BMD.
Some kinds of classic DMARDs might have effects on the progression of bone formation. For example, methotrexate (MTX) could inhibit the differentiation of osteoblasts and exert direct negative effects on bone metabolism in RA patients [36]. However, Minaur NJ et al. indicated that reduced BMD associated with MTX was due to confounders such as disease activity, and no adverse effect of low-dose MTX on bone formation in RA was detected [37]. Corticosteroids, widely used in treating RA due to their strong suppressive effect on inflammatory activity, can decrease generalized BMD loss. However, as a side-effect, they could also increase BMD loss [38]. Nevertheless, a study including 342 patients with RA found no differences in BMD loss between four common treatment strategies, including high doses of corticosteroids, sequential monotherapy/step-up combination therapy of high doses of MTX, and antitumor necrosis factor-α [15]. Moreover, lack of physical exercise might be another reason for bone loss in patients with RA. Patients with RA usually experience pain, swelling, and immobilization in one or more joints; thus, they are more likely to take less physical activity. Habitual levels of higher impact physical activity were reported to be positively related to lower limb bone strength in older women [39]. A study focused on healthy young men also found that sedentary activities were inversely related to FN-BMD [40].
As far as we know, no MR study on the effect of RA on OP/BMD or OP/BMD on RA has yet been reported. Our study uses several variants summarized from large-scale GWA studies on RA and BMD to date to increase the statistical power to detect causal associations. Compared with the analysis of individual-level data from a small study, a key strength of our research is that effect-size data were drawn from separate large-scale GWASs for exposure and outcome traits so that we can assess the effect sizes more precisely.
However, our study also has certain limitations. First, stratified analyses such as menopausal status would also have been of interest, due to the increased risk of developing OP of postmenopausal women [41], however, since we used summary-level data for two-sample MR analyses, the analyses in specific subgroups were not possible. Otherwise, female sex is also an independent risk factor for RA. However, a Korean large-scale observation study showed an increased risk of osteoporotic fractures for RA patients across all sexes, age groups, and various anatomic sites, suggesting that there might be no sex and age differences in bone loss in patients with RA [32]. Second, compared with other sites, hand BMD is an independent predictor of subsequent radiographic damage since quantitative hand bone loss in RA patients occurs before radiographic joint damage. Therefore, hand BMD may be used as an instrument for the evaluation of bone involvement in RA patients [42]. Unfortunately, we did not obtain available public GWAS data for hand BMD. The forearm BMD summarized data we used in our study included wrist but did not include hand. Nevertheless, some studies have suggested that hand BMD is associated with both the lumbar spine and total hip BMD among postmenopausal women with RA [43], and other RA patients [44]. Thus, using BMD in different sites as phenotypes for OP could reduce the statistical bias in our MR analyses.
In summary, we found no evidence that RA is causally associated with OP/BMD, or the other way around. Therefore, the associations between RA and OP/BMD indicated in previous observational studies are possibly related to secondary effects such as antirheumatic treatment and less physical activity. At present, the clinical treatment of RA is aimed at suppressing inflammatory activity and anti-osteoporosis. Based on our current results, we suggested no indications for anti-osteoporosis treatment in RA patients without other risk factors for OP.
Materials and Methods
Data sources
Summary statistics for RA were obtained from the MRC IEU OpenGWAS database (https://gwas.mrcieu.ac.uk/), which comprises mainly publicly available GWAS summary data, serving as an input source to a variety of analytical methods, such as Mendelian randomization, fine mapping, and colocalization [45, 46]. The GWAS for RA was examined in imputed genotype data from the UK Biobank study, which included 1,523 cases and 461,487 controls from European populations [45]. RA cases were obtained from the UK Biobank study using Hospital Episode Statistics (HES) with ICD-10 codes M06. Since the absence of the original individual measures, gender- or age-specific analyses were difficult to perform. Fortunately, gender and age were adjusted in the original GWAS analysis. GWAS summary statistics for heel, TB-, FA-, FN- and LS-BMD were obtained from GEFOS (http://www.gefos.org/). The GWAS for heel BMD included 265,627 individuals from the European population in the UK Biobank study. GWAS for heel BMD was estimated from quantitative heel ultrasounds, the age, sex, genotyping array, assessment center, and ancestry informative principal components 1 to 20 were included in the fixed model as covariates [47]. The GWAS dataset for TB-BMD contains summary statistics for a GWAS meta-analysis study involving 66,628 European individuals and was adjusted for age, weight, height, sex, genomic principal components, and other study-specific covariates (such as recruiting center) [23]. GWASs for FA-, FN-, LS-BMD were obtained from a meta-analysis released at GEFOS in 2015 [22]. Separately, 8,143 individuals for FA-BMD, 32,735 individuals for FN-BMD, and 28,498 individuals for LS-BMD from European populations were included, and BMD was adjusted for sex, weight, age, and age squared. More details for assessment, quality control, and association analysis were presented in the original studies [22, 23, 45, 47].
Selection of genetic variants
In the first stage, genetic variants associated with BMD were used as instrumental SNPs. To satisfy the 3 assumptions for MR analysis, we selected independent SNPs (r2 < 0.001) that were strongly (p < 5 × 10−8) associated with the exposure. Next, we obtained the association results of the corresponding SNPs with the outcome. If the corresponding SNPs were not available in the outcome data, we used proxy SNPs that were highly correlated (r2 > 0.8) with the corresponding SNPs (if possible). 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. In the second stage, since there were only 3 SNPs with a p-value less than 5 × 10−8 for RA, we broadened the threshold to 1 × 10-5 for selecting RA-associated variants as instrumental SNPs. To ensure that the selected IVs have enough power for detecting the causal effect of exposure on the outcome, we calculated the F statistic of selected IVs with an online tool (https://sb452.shinyapps.io/overlap) [48]. The selected IVs with F statistics >10 are considered powerful enough for the causal effect estimate.
MR analysis
IVW method was conducted as the primary method to estimate the causal effect between exposure and outcome, which was calculated as the effect size of the association between SNP and outcome divided by the effect size of the association between SNP and exposure [49]. When there was no evidence of directional pleiotropy (p for MR-Egger intercept > 0.05) among the selected IVs, the IVW method was considered the most reliable [26]. To ensure the robustness of our results, MR-Egger, weighted median, and weighted mode methods were also performed to estimate the causal effect of exposure on outcome. Detailed information about these MR methods mentioned above can be found in published studies [17, 18]. The MR analyses were performed in the R software (http://www.r-project.org) with the TwoSampleMR package [46]. For the IVs with a p-value < 0.05 for IVW analysis, we then performed MR-PRESSO with the MRPRESSO package [50], which can detect, remove the potential pleiotropic IVs (outliers), and provide the outlier-adjusted estimates. In addition, for pairs of exposure-outcomes with more than 100 IVs, we also performed a recently proposed MR method called Robust Adjusted Profile Score (RAPS) [51], which is unbiased even when there are many (such as hundreds of) weak instruments.
Sensitivity analysis
To further ensure the robustness of our MR estimates, the following sensitivity analyses were performed: First, Cochran’s Q statistics were employed to assess the heterogeneity among the IVs. Second, MR Egger regression was used to examine whether our MR analyses were driven by the directional horizontal pleiotropy. Moreover, to examine whether the casual association was driven by a single SNP, we performed the leave-one-out analysis.
Negative control
To further ensure the validity of the selected IVs, we included myopia as a negative control outcome, since no evidence showed that OP or RA is correlated with myopia. The GWAS data for myopia were derived from the FinnGen biobank (https://www.finngen.fi/en), including 621 myopia cases and 93,606 controls from the European population.
Author Contributions
Ying-Qi Liu: Conceptualization, Writing - Original Draft, Formal analysis Yong Liu: Software, Formal analysis Zhuo-Yuan Chen: Writing - Review and Editing Hui Li: Editing and Supervision Tao Xiao: Supervision.
Conflicts of Interest
The authors declared that they have no conflicts of interest to this work.
Funding
This article was funded by the National Nature Science Foundation of China (grant number 81502332), Research Project of Human Health Commission (grant number B2019162), and the Nature Science Foundation of Hunan (grant number 2019JJ50861).
References
-
1.
McInnes IB, Schett G. The pathogenesis of rheumatoid arthritis. N Engl J Med. 2011; 365:2205–19. https://doi.org/10.1056/NEJMra1004965 [PubMed]
-
2.
Helmick CG, Felson DT, Lawrence RC, Gabriel S, Hirsch R, Kwoh CK, Liang MH, Kremers HM, Mayes MD, Merkel PA, Pillemer SR, Reveille JD, Stone JH, and National Arthritis Data Workgroup. Estimates of the prevalence of arthritis and other rheumatic conditions in the United States. Part I. Arthritis Rheum. 2008; 58:15–25. https://doi.org/10.1002/art.23177 [PubMed]
-
3.
Chen L, Huang Z, Yang B, Cai B, Su Z, Wang L. Association of E26 Transformation Specific Sequence 1 Variants with Rheumatoid Arthritis in Chinese Han Population. PLoS One. 2015; 10:e0134875. https://doi.org/10.1371/journal.pone.0134875 [PubMed]
-
4.
Johnell O, Kanis JA, Oden A, Johansson H, De Laet C, Delmas P, Eisman JA, Fujiwara S, Kroger H, Mellstrom D, Meunier PJ, Melton LJ 3rd, O’Neill T, et al. Predictive value of BMD for hip and other fractures. J Bone Miner Res. 2005; 20:1185–94. https://doi.org/10.1359/JBMR.050304 [PubMed]
-
5.
Lodder MC, de Jong Z, Kostense PJ, Molenaar ET, Staal K, Voskuyl AE, Hazes JM, Dijkmans BA, Lems WF. Bone mineral density in patients with rheumatoid arthritis: relation between disease severity and low bone mineral density. Ann Rheum Dis. 2004; 63:1576–80. https://doi.org/10.1136/ard.2003.016253 [PubMed]
-
6.
Cummings SR, Bates D, Black DM. Clinical use of bone densitometry: scientific review. JAMA. 2002; 288:1889–97. https://doi.org/10.1001/jama.288.15.1889 [PubMed]
-
7.
MacGregor AJ, Snieder H, Rigby AS, Koskenvuo M, Kaprio J, Aho K, Silman AJ. Characterizing the quantitative genetic contribution to rheumatoid arthritis using data from twins. Arthritis Rheum. 2000; 43:30–37. https://doi.org/10.1002/1529-0131(200001)43:1<30::AID-ANR5>3.0.CO;2-B [PubMed]
-
8.
Ioannidis JP, Ng MY, Sham PC, Zintzaras E, Lewis CM, Deng HW, Econs MJ, Karasik D, Devoto M, Kammerer CM, Spector T, Andrew T, Cupples LA, et al. Meta-analysis of genome-wide scans provides evidence for sex- and site-specific regulation of bone mass. J Bone Miner Res. 2007; 22:173–83. https://doi.org/10.1359/jbmr.060806 [PubMed]
-
9.
Richards JB, Zheng HF, Spector TD. Genetics of osteoporosis from genome-wide association studies: advances and challenges. Nat Rev Genet. 2012; 13:576–88. https://doi.org/10.1038/nrg3228 [PubMed]
-
10.
Okada Y, Wu D, Trynka G, Raj T, Terao C, Ikari K, Kochi Y, Ohmura K, Suzuki A, Yoshida S, Graham RR, Manoharan A, Ortmann W, et al, and RACI consortium, and GARNET consortium. Genetics of rheumatoid arthritis contributes to biology and drug discovery. Nature. 2014; 506:376–81. https://doi.org/10.1038/nature12873 [PubMed]
-
11.
Cagnetta V, Patella V. The role of the immune system in the physiopathology of osteoporosis. Clin Cases Miner Bone Metab. 2012; 9:85–88. [PubMed]
-
12.
Corvaisier M, Delneste Y, Jeanvoine H, Preisser L, Blanchard S, Garo E, Hoppe E, Barré B, Audran M, Bouvard B, Saint-André JP, Jeannin P. IL-26 is overexpressed in rheumatoid arthritis and induces proinflammatory cytokine production and Th17 cell generation. PLoS Biol. 2012; 10:e1001395. https://doi.org/10.1371/journal.pbio.1001395 [PubMed]
-
13.
Güler-Yüksel M, Bijsterbosch J, Goekoop-Ruiterman YP, de Vries-Bouwstra JK, Ronday HK, Peeters AJ, de Jonge-Bok JM, Breedveld FC, Dijkmans BA, Allaart CF, Lems WF. Bone mineral density in patients with recently diagnosed, active rheumatoid arthritis. Ann Rheum Dis. 2007; 66:1508–12. https://doi.org/10.1136/ard.2007.070839 [PubMed]
-
14.
Güler-Yüksel M, Allaart CF, Goekoop-Ruiterman YP, de Vries-Bouwstra JK, van Groenendael JH, Mallée C, de Bois MH, Breedveld FC, Dijkmans BA, Lems WF. Changes in hand and generalised bone mineral density in patients with recent-onset rheumatoid arthritis. Ann Rheum Dis. 2009; 68:330–36. https://doi.org/10.1136/ard.2007.086348 [PubMed]
-
15.
Güler-Yüksel M, Bijsterbosch J, Goekoop-Ruiterman YP, de Vries-Bouwstra JK, Hulsmans HM, de Beus WM, Han KH, Breedveld FC, Dijkmans BA, Allaart CF, Lems WF. Changes in bone mineral density in patients with recent onset, active rheumatoid arthritis. Ann Rheum Dis. 2008; 67:823–28. https://doi.org/10.1136/ard.2007.073817 [PubMed]
-
16.
Mohammad A, Lohan D, Bergin D, Mooney S, Newell J, O’Donnell M, Coughlan RJ, Carey JJ. The prevalence of vertebral fracture on vertebral fracture assessment imaging in a large cohort of patients with rheumatoid arthritis. Rheumatology (Oxford). 2014; 53:821–27. https://doi.org/10.1093/rheumatology/ket353 [PubMed]
-
17.
Davey Smith G, Hemani G. Mendelian randomization: genetic anchors for causal inference in epidemiological studies. Hum Mol Genet. 2014; 23:R89–98. https://doi.org/10.1093/hmg/ddu328 [PubMed]
-
18.
Didelez V, Sheehan N. Mendelian randomization as an instrumental variable approach to causal inference. Stat Methods Med Res. 2007; 16:309–30. https://doi.org/10.1177/0962280206077743 [PubMed]
-
19.
Lawlor DA. Commentary: Two-sample Mendelian randomization: opportunities and challenges. Int J Epidemiol. 2016; 45:908–15. https://doi.org/10.1093/ije/dyw127 [PubMed]
-
20.
Kanis JA. Diagnosis of osteoporosis and assessment of fracture risk. Lancet. 2002; 359:1929–36. https://doi.org/10.1016/S0140-6736(02)08761-5 [PubMed]
-
21.
Iguacel I, Miguel-Berges ML, Gómez-Bruton A, Moreno LA, Julián C. Veganism, vegetarianism, bone mineral density, and fracture risk: a systematic review and meta-analysis. Nutr Rev. 2019; 77:1–18. https://doi.org/10.1093/nutrit/nuy045 [PubMed]
-
22.
Zheng HF, Forgetta V, Hsu YH, Estrada K, Rosello-Diez A, Leo PJ, Dahia CL, Park-Min KH, Tobias JH, Kooperberg C, Kleinman A, Styrkarsdottir U, Liu CT, et al, and AOGC Consortium, and UK10K Consortium. Whole-genome sequencing identifies EN1 as a determinant of bone density and fracture. Nature. 2015; 526:112–17. https://doi.org/10.1038/nature14878 [PubMed]
-
23.
Medina-Gomez C, Kemp JP, Trajanoska K, Luan J, Chesi A, Ahluwalia TS, Mook-Kanamori DO, Ham A, Hartwig FP, Evans DS, Joro R, Nedeljkovic I, Zheng HF, et al. Life-Course Genome-wide Association Study Meta-analysis of Total Body BMD and Assessment of Age-Specific Effects. Am J Hum Genet. 2018; 102:88–102. https://doi.org/10.1016/j.ajhg.2017.12.005 [PubMed]
-
24.
Kemp JP, Morris JA, Medina-Gomez C, Forgetta V, Warrington NM, Youlten SE, Zheng J, Gregson CL, Grundberg E, Trajanoska K, Logan JG, Pollard AS, Sparkes PC, et al. Identification of 153 new loci associated with heel bone mineral density and functional involvement of GPC6 in osteoporosis. Nat Genet. 2017; 49:1468–75. https://doi.org/10.1038/ng.3949 [PubMed]
-
25.
Bowden J, Del Greco M F, Minelli C, Davey Smith G, Sheehan NA, Thompson JR. Assessing the suitability of summary data for two-sample Mendelian randomization analyses using MR-Egger regression: the role of the I2 statistic. Int J Epidemiol. 2016; 45:1961–74. https://doi.org/10.1093/ije/dyw220 [PubMed]
-
26.
Holmes MV, Ala-Korpela M, Smith GD. Mendelian randomization in cardiometabolic disease: challenges in evaluating causality. Nat Rev Cardiol. 2017; 14:577–90. https://doi.org/10.1038/nrcardio.2017.78 [PubMed]
-
27.
Kweon SM, Sohn DH, Park JH, Koh JH, Park EK, Lee HN, Kim K, Kim Y, Kim GT, Lee SG. Male patients with rheumatoid arthritis have an increased risk of osteoporosis: Frequency and risk factors. Medicine (Baltimore). 2018; 97:e11122. https://doi.org/10.1097/MD.0000000000011122 [PubMed]
-
28.
Nava-Valdivia CA, Saldaña-Cruz AM, Corona-Sanchez EG, Murillo-Vazquez JD, Moran-Moguel MC, Salazar-Paramo M, Perez-Guerrero EE, Vazquez-Villegas ML, Bonilla-Lara D, Rocha-Muñoz AD, Martín-Marquez BT, Sandoval-Garcia F, Martínez-García EA, et al. Polymorphism rs2073618 of the TNFRSF11B (OPG) Gene and Bone Mineral Density in Mexican Women with Rheumatoid Arthritis. J Immunol Res. 2017; 2017:7680434. https://doi.org/10.1155/2017/7680434 [PubMed]
-
29.
Abu-Shakra M, Zisman D, Balbir-Gurman A, Amital H, Levy Y, Langevitz P, Tishler M, Molad Y, Aamar S, Roser I, Avshovich N, Paran D, Reitblat T, et al. Effect of Tocilizumab on Fatigue and Bone Mineral Density in Patients with Rheumatoid Arthritis. Isr Med Assoc J. 2018; 20:239–44. [PubMed]
-
30.
Di Munno O, Ferro F. The effect of biologic agents on bone homeostasis in chronic inflammatory rheumatic diseases. Clin Exp Rheumatol. 2019; 37:502–07. [PubMed]
-
31.
Tong JJ, Xu SQ, Zong HX, Pan MJ, Teng YZ, Xu JH. Prevalence and risk factors associated with vertebral osteoporotic fractures in patients with rheumatoid arthritis. Clin Rheumatol. 2020; 39:357–64. https://doi.org/10.1007/s10067-019-04787-9 [PubMed]
-
32.
Kim SY, Schneeweiss S, Liu J, Daniel GW, Chang CL, Garneau K, Solomon DH. Risk of osteoporotic fracture in a large population-based cohort of patients with rheumatoid arthritis. Arthritis Res Ther. 2010; 12:R154. https://doi.org/10.1186/ar3107 [PubMed]
-
33.
Deodhar AA, Brabyn J, Pande I, Scott DL, Woolf AD. Hand bone densitometry in rheumatoid arthritis, a five year longitudinal study: an outcome measure and a prognostic marker. Ann Rheum Dis. 2003; 62:767–70. https://doi.org/10.1136/ard.62.8.767 [PubMed]
-
34.
Yeo L, Toellner KM, Salmon M, Filer A, Buckley CD, Raza K, Scheel-Toellner D. Cytokine mRNA profiling identifies B cells as a major source of RANKL in rheumatoid arthritis. Ann Rheum Dis. 2011; 70:2022–28. https://doi.org/10.1136/ard.2011.153312 [PubMed]
-
35.
Gravallese EM, Manning C, Tsay A, Naito A, Pan C, Amento E, Goldring SR. Synovial tissue in rheumatoid arthritis is a source of osteoclast differentiation factor. Arthritis Rheum. 2000; 43:250–58. https://doi.org/10.1002/1529-0131(200002)43:2<250::AID-ANR3>3.0.CO;2-P [PubMed]
-
36.
Scheven BA, van der Veen MJ, Damen CA, Lafeber FP, Van Rijn HJ, Bijlsma JW, Duursma SA. Effects of methotrexate on human osteoblasts in vitro: modulation by 1,25-dihydroxyvitamin D3. J Bone Miner Res. 1995; 10:874–80. https://doi.org/10.1002/jbmr.5650100608 [PubMed]
-
37.
Minaur NJ, Kounali D, Vedi S, Compston JE, Beresford JN, Bhalla AK. Methotrexate in the treatment of rheumatoid arthritis. II. In vivo effects on bone mineral density. Rheumatology (Oxford). 2002; 41:741–49. https://doi.org/10.1093/rheumatology/41.7.741 [PubMed]
-
38.
Habib GS, Haj S. Bone mineral density in patients with early rheumatoid arthritis treated with corticosteroids. Clin Rheumatol. 2005; 24:129–33. https://doi.org/10.1007/s10067-004-0989-1 [PubMed]
-
39.
Hannam K, Deere KC, Hartley A, Al-Sari UA, Clark EM, Fraser WD, Tobias JH. Habitual levels of higher, but not medium or low, impact physical activity are positively related to lower limb bone strength in older women: findings from a population-based study using accelerometers to classify impact magnitude. Osteoporos Int. 2017; 28:2813–22. https://doi.org/10.1007/s00198-016-3863-5 [PubMed]
-
40.
Weeda J, Horan S, Beck B, Weeks BK. Lifetime physical activity, neuromuscular performance and body composition in healthy young men. Int J Sports Med. 2014; 35:900–05. https://doi.org/10.1055/s-0033-1364027 [PubMed]
-
41.
Curtis E, Litwic A, Cooper C, Dennison E. Determinants of Muscle and Bone Aging. J Cell Physiol. 2015; 230:2618–25. https://doi.org/10.1002/jcp.25001 [PubMed]
-
42.
Hoff M, Haugeberg G, Odegård S, Syversen S, Landewé R, van der Heijde D, Kvien TK. Cortical hand bone loss after 1 year in early rheumatoid arthritis predicts radiographic hand joint damage at 5-year and 10-year follow-up. Ann Rheum Dis. 2009; 68:324–29. https://doi.org/10.1136/ard.2007.085985 [PubMed]
-
43.
Desai SP, Gravallese EM, Shadick NA, Glass R, Cui J, Frits M, Chibnik LB, Maher N, Weinblatt ME, Solomon DH. Hand bone mineral density is associated with both total hip and lumbar spine bone mineral density in post-menopausal women with RA. Rheumatology (Oxford). 2010; 49:513–19. https://doi.org/10.1093/rheumatology/kep385 [PubMed]
-
44.
Naumann L, Hermann KG, Huscher D, Lenz K, Burmester GR, Backhaus M, Buttgereit F. Quantification of periarticular demineralization and synovialitis of the hand in rheumatoid arthritis patients. Osteoporos Int. 2012; 23:2671–79. https://doi.org/10.1007/s00198-012-1897-x [PubMed]
-
45.
Elsworth B, Lyon M, Alexander T, Liu Y, Matthews P, Hallett J, Bates P, Palmer T, Haberland V, Smith GD, Zheng J, Haycock P, Gaunt TR, et al. The MRC IEU OpenGWAS data infrastructure. bioRxiv. 2020. https://doi.org/10.1101/2020.08.10.244293
-
46.
Hemani G, Zheng J, Elsworth B, Wade KH, Haberland V, Baird D, Laurin C, Burgess S, Bowden J, Langdon R, Tan VY, Yarmolinsky J, Shihab HA, et al. The MR-Base platform supports systematic causal inference across the human phenome. Elife. 2018; 7:e34408. https://doi.org/10.7554/eLife.34408 [PubMed]
-
47.
Morris JA, Kemp JP, Youlten SE, Laurent L, Logan JG, Chai RC, Vulpescu NA, Forgetta V, Kleinman A, Mohanty ST, Sergio CM, Quinn J, Nguyen-Yamamoto L, et al, and 23andMe Research Team. An atlas of genetic influences on osteoporosis in humans and mice. Nat Genet. 2019; 51:258–66. https://doi.org/10.1038/s41588-018-0302-x [PubMed]
-
48.
Palmer TM, Lawlor DA, Harbord RM, Sheehan NA, Tobias JH, Timpson NJ, Davey Smith G, Sterne JA. Using multiple genetic variants as instrumental variables for modifiable risk factors. Stat Methods Med Res. 2012; 21:223–42. https://doi.org/10.1177/0962280210394459 [PubMed]
-
49.
Burgess S, Butterworth A, Thompson SG. Mendelian randomization analysis with multiple genetic variants using summarized data. Genet Epidemiol. 2013; 37:658–65. https://doi.org/10.1002/gepi.21758 [PubMed]
-
50.
Verbanck M, Chen CY, Neale B, Do R. Detection of widespread horizontal pleiotropy in causal relationships inferred from Mendelian randomization between complex traits and diseases. Nat Genet. 2018; 50:693–98. https://doi.org/10.1038/s41588-018-0099-7 [PubMed]
-
51.
Qingyuan Zhao JW, Bowden J, Small DS. Statistical inference in two-sample summary-data Mendelian randomization using robust adjusted profi. arXiv.