Sage Journals: Discover world-class research

Abstract

Background:

The co-occurrence of osteoporosis (OP) and rheumatoid arthritis (RA) has long been observed; their intrinsic link, however, has not been fully understood.

Objectives:

We aimed to inform the importance of integrated care targeting both diseases by investigating the phenotypic as well as the genetic relationships underlying OP and RA.

Design:

This is a prospective cohort study and genome-wide cross-trait analysis.

Methods:

We evaluated phenotypic associations using longitudinal follow-up data from the UK Biobank (N = 472,050). We investigated genetic relationships by leveraging summary statistics from the largest genome-wide association study in European ancestry conducted for RA (N_case/control = 22,350/74,823), seropositive RA (N_case/control = 17,221/74,823), and a reliable proxy of OP—the heel estimated bone mineral density (eBMD; N = 426,824).

Results:

Observational analysis suggested a bidirectional relationship (OP → RA: hazard ratio (HR) = 1.59, 95% confidence interval (CI) = 1.28–1.97; RA → OP: HR = 2.94, 95% CI = 2.59–3.35). A negative overall genetic correlation was observed for eBMD with RA ( = −0.06, p = 1.37 × 10⁻⁵) and with its seropositive subtype (r_g = −0.06, p = 9.15 × 10⁻⁶). Cross-trait meta-analysis replicated 96 previously reported trait-associated loci and discovered three novel pleiotropic loci (rs72836346, rs2613812, and rs76458888). Transcriptome-wide association study revealed 23 shared genes. Mendelian randomization analysis suggested a putative causal effect of eBMD on RA (odds ratio (OR) = 0.90, 95% CI = 0.84–0.96), but not of RA on eBMD.

Conclusion:

Our work demonstrates a significant biological pleiotropy as well as a putative causal relationship between OP and RA, emphasizing an intrinsic link underlying the pronounced phenotypic association. These findings highlight the possibility of preventing and predicting RA development by monitoring and interfering with bone loss in preclinical high-risk individuals.

Plain language summary

Investigating the relationship between osteoporosis and rheumatoid arthritis

• Observational analysis suggested a bidirectional relationship between OP and RA.

• Genetic analyses observed a significant negative genetic correlation for eBMD with RA and with its seropositive subtype, further substantiated by seven local genomic signals, 99 shared independent loci (3 novel loci) and 139 shared tissue-gene pairs.

• Mendelian randomization analysis suggested a putative causal effect of a decreased eBMD on an increased risk of RA, but not of RA on eBMD.

• These findings highlight the possibility of preventing RA development by monitoring bone health and interfering with BMD decline in preclinical high-risk individuals.

Keywords

genetic correlation Mendelian randomization osteoporosis phenotypic association pleiotropic loci rheumatoid arthritis

Introduction

Osteoporosis (OP) and rheumatoid arthritis (RA) have been important drivers of an increasing global burden of disability over the past 30 years.^1,2 The co-occurrence of these two common musculoskeletal disorders has long been observed with evidence of both phenotypic and genetic connections.^3–9 Findings from previous case–control studies have suggested a bidirectional association between OP (defined as a T-score of bone mineral density (BMD) below −2.5 standard deviations measured by dual-energy X-ray absorptiometry) and RA, including a significant association of a history of OP with RA (odds ratio (OR) = 1.43, 95% confidence interval (CI) = 1.32–1.55) by comparing 28,341 incident RA cases with 283,226 matched controls⁴ and a significant association of a history of RA with OP (OR = 6.58, 95% CI = 2.29–18.91) by comparing 135 incident OP cases with 135 matched controls.⁵ Downstream analysis of genome-wide association studies (GWASs) has further identified a significant genetic correlation between lumbar spine BMD and RA ( $r_{g}$ = −0.13),¹⁰ implicating the phenotypic association partly attributable to shared genetic basis. Multiple independent loci influencing both traits have also been observed (e.g., HIPK1, CTLA4, and ANKRD55).⁷ Mendelian randomization (MR) has confirmed a bidirectional causal association between estimated BMD (eBMD) measured by quantitative heel ultrasound and RA (eBMD → RA: beta = −7.57 × 10⁻⁴, 95% CI = −0.0014 to −0.0001⁶; RA → eBMD: beta = −0.021, 95% CI = −0.030 to −0.011⁷), as well as a putative causal effect of RA on OP (OR = 1.10, 95% CI = 1.06–1.14).⁹

Despite knowledge gained from observational and genetic analyses advancing our understanding of the co-existence of OP and RA, recommendations targeting both disorders in clinical guidelines remain sparse due to several unaddressed gaps.^11,12 First, existing observational studies lacked prospective follow-up, prone to recall bias or reverse causality.^4,5 Second, previous MR analyses leveraged GWAS data from completely overlapping samples,^6,7,9 of small sample sizes,^6,7,9 or with severe case–control imbalance,^6–8 substantially restricting the statistical power to detect a causal effect. Nonetheless, a recent MR reported no association of genetically predicted RA with eBMD (beta = −0.006, 95% CI = −0.015 to 0.003, p = 0.20), casting doubt on the OP–RA causal association.⁸ Third, previous MR analyses have not examined subgroups of RA defined by serological status,^6,8,9 an important aspect of high clinical relevance.¹³ Fourth, almost no cross-phenotype association (CPASSOC) analysis has been implemented to systemically search the whole genome for potential pleiotropic variants beyond a few well-known loci nor examined their biological mechanisms.^6,8,9 To our knowledge, no large-scale genetic analysis has been performed to comprehensively evaluate the magnitude and direction of shared etiology underlying OP and RA.

Therefore, the current study aimed to comprehensively interrogate the relationship between OP and RA, with overarching goal of providing genetic insights into the observed associations. These insights may contribute to advancing precision prevention as well as precision medicine for common musculoskeletal diseases.

Methods

The overall study design is shown in Supplemental Figure 1. We first investigated the bidirectional phenotypic association using more than 470,000 participants available in the UK Biobank (UKB). We next performed a genome-wide cross-trait analysis to quantify the overall and local genetic correlations, identify pleiotropic loci, detect tissue-expression associations, and make causal inferences. This study followed the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement,¹⁴ the STrengthening the REporting of Genetic Association Studies (STREGA) statement,¹⁵ and the Strengthening the Reporting of Observational Studies in Epidemiology using Mendelian Randomization (STROBE-MR) statement.^16,17

Data sources

UKB data

UKB is a population-based cohort with more than 500,000 participants aged 40–69 years recruited in 22 assessment centers across the United Kingdom between 2006 and 2010.¹⁸ For the current study, we identified a subset of 472,050 individuals of White ancestry. We defined OP as self-reported medical conditions code 1309 at baseline, the International Classification of Diseases, Ninth Revision (ICD-9) code 7330 and ICD-10 codes M80, M81, M82, and RA as self-reported medical conditions code 1464 at baseline, ICD-9 code 7140 and ICD-10 codes M05, M06. These combined definitions follow the algorithms developed by the UK Biobank Outcome Adjudication Group and were used to minimize bias toward the null arising from case misclassification.

For the analysis of the association between OP and incident RA, we excluded 6078 participants with a history of RA at baseline, resulting in a final sample of 465,972 individuals. For the analysis of the association between RA and incident OP, we excluded 9416 participants with a history of OP at baseline, yielding a final sample of 462,634 individuals.

GWAS summary statistics

EBMD GWAS

Heel eBMD, the most heritable site among all skeletal sites, has been recognized as the most prominent predictor of OP.^19–21 Given the limited power of the hitherto available GWAS on OP, we selected the GWAS of eBMD in European ancestry. The hitherto largest GWAS of eBMD was conducted by the Genetic Factors for Osteoporosis Consortium in 426,824 UKB participants.¹⁹ eBMD was assessed by quantitative ultrasound speed of sound and broadband ultrasound attenuation using a Sahara Clinical Bone Sonometer (Hologic Corporation, Bedford, MA, USA). A total of 1097 independent eBMD-associated single-nucleotide polymorphisms (SNPs) at genome-wide significance (p < 6.6 × 10⁻⁹, $r^{2}$ ⩽ 0.1 across a 20 Mb window) were identified as instrumental variables (IVs). We extracted the effect size and relevant information of these IVs (Supplemental Tables 1 and 2) as well as retrieved the full set GWAS summary statistics.

RA GWAS

The largest GWAS of RA in European ancestry was conducted by Ishigaki et al.,²² meta-analyzing data from 25 studies comprising 22,350 cases and 74,823 controls. Since approximately two-thirds of individuals with RA are seropositive for rheumatoid factor and/or autoantibodies against citrullinated peptide antigens,²³ to confirm the robustness of findings, we also included the most important type of RA in clinical practice, seropositive RA. The largest GWAS of seropositive RA in European ancestry was also conducted by Ishigaki et al., meta-analyzing data from 25 studies comprising 17,221 cases and 74,823 controls. A total of 57 independent RA-associated SNPs and 61 independent seropositive RA-associated SNPs at genome-wide significance (p < 5 × 10⁻⁸, r² ⩽ 0.1 across a 1 Mb window) were identified as IVs. We extracted the effect size and relevant information of theses IVs (Supplemental Tables 3 and 4) as well as retrieved the full set GWAS summary statistics.

Statistical analysis

Observational analysis

We first assessed the association between OP and the subsequent incidence of RA. Person-years at risk were calculated from the date of attending assessment center until the date of RA diagnosis, death, loss to follow-up, or end of follow-up, whichever occurred first. We fitted a Cox proportional hazards regression model to estimate the hazards ratio (HR) for incident RA associated with prior OP diagnosis with incremental adjustments as follows: model 1 (crude model) without any adjustments; model 2 adjusted for sex, age, assessment center, income, Townsend deprivation index, and the first 10 genetic principal components; model 3 adjusted on top of the model 2 covariates for smoking, drinking, physical activity, sleep duration, and body mass index; and model 4 adjusted on top of the model 3 covariates for diabetes mellitus, hypertension, dyslipidemia, use of antidiabetic drugs (Anatomical Therapeutic Chemical classification code A10; data field 20003), antihypertensive drugs (C02), and HMG-CoA reductase inhibitors (C10AA). In the sensitivity analysis, we excluded participants who had less than 2 years of follow-up or a diagnosis of RA within 2 years after follow-up.^24,25 We repeated all these analyses to assess the association between RA and the subsequent incidence of OP. The same covariate structure was applied, with Model 4 further adjusted for the use of non-steroidal anti-inflammatory and antirheumatic drugs (M01A) and glucocorticoids (R03A). We assessed the proportional hazards assumption using Schoenfeld residuals and identified violations for the bidirectional associations between OP and RA (p < 0.05). Therefore, the HRs should be interpreted as representing the weighted average effect over the follow-up period.²⁶ All analyses were performed using SAS 9.4 software (SAS Institute, Cary, NC, USA). A two-sided p value < 0.05 was considered statistically significant.

Global and local genetic correlation analysis

To assess the overall shared genetics, we estimated the global genetic correlation using the cross-trait linkage disequilibrium (LD) score regression.²⁷ This algorithm relies on the fact that the effect for a given SNP aggregates the effects of all SNPs in LD with that SNP. Genetic correlation ranges from −1 to 1, where 1 represents complete positive genetic correlation. A Bonferroni-corrected p value (p < 0.025 = 0.05/2) was considered statistically significant.

While global genetic correlation represents an average of shared genetic association across the genome, the underlying architecture of correlations at local genomic regions or individual loci can vary. We thus estimated the local genetic correlation in 2353 predefined LD-independent regions (average length 1.6 cM) using SUPERGNOVA.²⁸ A Bonferroni-corrected p value (p < 2.12×10⁻⁵ = 0.05/2353) was considered statistically significant.

Cross-trait meta-analysis

To identify potential pleiotropic loci, we performed a cross-trait meta-analysis using CPASSOC through the statistic S_het.²⁹ CPASSOC integrates GWAS summary statistics from multiple correlated traits to detect variants with shared association, controlling for population structure or cryptic correlation. To obtain independent shared SNPs, we applied PLINK clumping function (parameters: –clump-p1 5e-8 –clump-p2 1e-5 –clump-r2 0.2 –clump-kb 500). Variants with p_CPASSOC < 5 × 10⁻⁸ and $p_{e B M D o r R A G W A S}$ < 1 × 10⁻⁵ were identified as significant pleiotropic SNPs.

Specifically, all significant pleiotropic SNPs can be classified into one of the four categories: (i) a “known” shared SNP, defined as a SNP achieved genome-wide significant association with both single traits; (ii) a “single-trait-driven” shared SNP, defined as a SNP achieved genome-wide significant association with one of the two single traits; (iii) a “LD-tagged” shared SNP, defined as a shared SNP in LD with index SNPs previously identified by single-trait GWAS (r² ⩾ 0.05), despite not driven by any single trait; and (iv) a “novel” shared SNP, defined as a shared SNP neither driven by any single trait nor in LD with index SNPs previously identified by single-trait GWAS (r² < 0.05).

To gain detailed functional annotation for pleiotropic loci, we applied Ensembl Variant Effect Predictor to map SNPs to their nearest genes. To investigate whether the same causal variant is responsible for two GWAS signals as opposed to distinct causal variants in proximity, we performed a colocalization analysis using Coloc.³⁰ We extracted summary statistics for variants within 500 kb of each shared index SNP and estimated the posterior probability for H4 (PPH4, the probability that both traits share a causal variant). A locus was considered colocalized if PPH4 was greater than 0.8.

Transcriptome-wide association study analysis

To identify pleiotropic genes whose expression pattern across tissues implicates shared etiology or biological mechanisms, we performed transcriptome-wide association study (TWAS) analyses using FUSION.³¹ This algorithm integrates GWAS summary statistics with precomputed gene expression weights to evaluate the association of each gene with the disease. First, we performed single-trait TWAS using precomputed expression reference weights from 49 postmortem Genotype-Tissue Expression (GTEx) project tissues. The Bonferroni correction within each tissue ( $p_{B o n f e r r o n i}$ < 0.05) was used to identify significant expression–trait associations. Second, we performed joint/conditional tests for loci with multiple associated features to identify independent genes at each locus. Third, we performed a colocalization analysis to test whether GWAS signals and GTEx expression quantitative trait loci (eQTLs) signals were colocalized at the same causal variant. Fourth, we orchestrated single-trait TWAS results to identify shared tissue–gene pairs across traits.

MR analysis

To assess potential causal relationships, we performed a bidirectional two-sample MR analysis using TwoSampleMR. We applied the inverse-variance weighted (IVW) method³² as our primary approach, which assumes that all IVs are valid and provides the greatest statistical power. To evaluate the robustness of the findings and to assess potential violations of the independence and exclusion restriction assumptions, we implemented several sensitivity analyses using the weighted-median,³³ MR-Pleiotropy Residual Sum and Outlier (MR-PRESSO),³⁴ and MR-Egger regression method.³⁵ We additionally repeated the IVW analyses after excluding palindromic IVs or pleiotropic IVs (SNPs associated with potential confounders according to GWAS Catalog). A causal relationship was considered significant if the estimate in IVW was statistically significant (p < 0.0125 = 0.05/4) and showed directional consistency across all sensitivity methods. Finally, we estimated the phenotypic variance explained by IVs (R² $R^{2}$ ) as well as calculated the F-statistics to test the strength of IVs, thereby ensuring that the MR “relevance” assumption was satisfied.

Results

Phenotypic association

Baseline characteristics of UKB participants were presented in Supplemental Tables 5 and 6.

In the analysis for the risk of incident RA associated with OP, participants were followed for 5,646,311 person-years (12.1 ± 2.0 years), during which 210 individuals with OP and 4730 individuals without OP developed RA (Table 1). Individuals with OP had a significantly increased risk of RA compared with individuals without OP (HR = 2.36, 95% CI = 2.05–2.71). After adjusting for sex, age, assessment center, income, Townsend deprivation index, and genetic principal components, the effect attenuated to some extent but remained significant (HR = 1.46, 95% CI = 1.23–1.72). Further adjustments for smoking, drinking, physical activity, sleep duration, and body mass index did not substantially change the effect (HR = 1.56, 95% CI = 1.29–1.89). In the fully adjusted model, the effect stabilized to 1.59 (95% CI = 1.28–1.97). A similar result was observed in the sensitivity analysis (HR = 1.48, 95% CI = 1.17–1.86).

Table 1.

Observational associations between OP and RA.

Exposure status at baseline	Cases/person-years	Primary analysis				Sensitivity analysis
Exposure status at baseline	Cases/person-years	Model 1	Model 2	Model 3	Model 4	Sensitivity analysis
OP → RA
No	4730/5,541,673	1.00 (ref)	1.00 (ref)	1.00 (ref)	1.00 (ref)	1.00 (ref)
Yes	210/104,638	2.36 (2.05, 2.71)	1.46 (1.23, 1.72)	1.56 (1.29, 1.89)	1.59 (1.28, 1.97)	1.48 (1.17, 1.86)
RA → OP
No	11,330/5,516,206	1.00 (ref)	1.00 (ref)	1.00 (ref)	1.00 (ref)	1.00 (ref)
Yes	566/64,153	4.38 (4.03, 4.77)	3.16 (2.86, 3.48)	3.06 (2.73, 3.43)	2.94 (2.59, 3.35)	2.86 (2.50, 3.27)

HRs are provided with 95% CIs. Model 1: unadjusted (crude) model; Model 2: adjusted for sex, age, assessment center, income, Townsend deprivation index, the first 10 genetic principal components; Model 3: model 2 plus adjusted for smoking, drinking, physical activity, sleep duration, body mass index; Model 4: model 3 plus adjusted for diabetes mellitus, hypertension, dyslipidemia, use of antidiabetic drugs, antihypertensive drugs, HMG-CoA reductase inhibitors, anti-inflammatory and antirheumatic products, non-steroids and glucocorticoids (just adjusted in the analysis for the risk of incident OP associated with RA).

CI, confidence intervals; HR, hazard ratio; OP, osteoporosis; RA, rheumatoid arthritis.

In the analysis for the risk of incident OP associated with RA, participants were followed for 5,580,359 person-years (12.1 ± 2.0 years), during which 566 individuals with RA and 11,330 individuals without RA developed OP. Individuals with RA had a significantly increased risk of OP compared with individuals without RA (HR = 4.38, 95% CI = 4.03–4.77). After adjusting for sex, age, assessment center, income, Townsend deprivation index, and genetic principal components, the effect weakened but remained significant (HR = 3.16, 95% CI = 2.86–3.48). Further adjustments for smoking, drinking, physical activity, sleep duration, and body mass index had minimal impact on the effect (HR = 3.06, 95% CI = 2.73–3.43). In the fully adjusted model, the effect stabilized to 2.94 (95% CI = 2.59–3.35). A similar result was observed in the sensitivity analysis (HR = 2.86, 95% CI = 2.50–3.27).

Global and local genetic correlation

We observed a significant negative global genetic correlation for eBMD with RA ( $r_{g}$ = −0.06, p = 1.37 × 10⁻⁵) and with its seropositive subtype ( $r_{g}$ = −0.06, p = 9.15 × 10⁻⁶).

Partitioning the whole genome into 2353 LD-independent regions, we identified four genomic regions (1p22.3-p22.2, 15q26.1, 20p12.3, and 21q22.13) presenting a significant local genetic correlation for eBMD and RA (Figure 1(a)–(d)). The locus PKN2 on 1p22.3-p22.2, C20orf196 and CASC20 on 20p12.3, and CLDN14 on 21q22.13 were known to be associated with eBMD.¹⁹ Notably, 15q26.1 that harbors a previously reported eBMD-associated gene IQGAP1 remained significant for eBMD and seropositive RA.¹⁹ Three additional genomic regions (2p23.2, 13q32.2-q32.3, and 15q23) showed a significant local signal specific to eBMD and seropositive RA. The loci PLB1 and PPP1CB on 2p23.2, and DOCK9 on 13q32.2-q32.3 were known to be associated with eBMD.¹⁹ Within region 15q23, TLE3 and LOC101929151 were known to be associated with eBMD,¹⁹ while PCAT29 and LINC00593 were known to be associated with seropositive RA.²²

Figure 1.

Genome-wide local genetic correlation between eBMD and RA. In the quantile–quantile plots (a, c) and Manhattan plots (b, d), red points represent genomic regions that contribute significant local genetic correlation as estimated by SUPERGNOVA, with significance defined by the multiple-testing threshold (p < 0.05/2353).

Cross-trait meta-analysis

Through performing pairwise CPASSOC, we identified a total of 99 independent pleiotropic loci ( $p_{C P A S S O C}$ < 5 × 10⁻⁸ and $p_{e B M D o r R A G W A S}$ < 1 × 10⁻⁵), including 49 loci shared between eBMD and RA, and 50 loci shared between eBMD and seropositive RA (Figure 2 and Supplemental Table 7). Of note, 32 loci were shared across all traits.

Figure 2.

Cross-trait meta-analysis between eBMD and RA. (a) Circular Manhattan plot between eBMD and RA. The outermost circle shows the cross-trait meta-analysis results between eBMD and RA; from the periphery to the center, each circle shows the GWAS results on RA and eBMD, respectively. The light blue indicates variants with genome-wide significance (P^eBMD $p_{e B M D}$ < 6.6 × 10⁻⁹ or P_RA < 5 × 10⁻⁸), whereas the dark blue indicates variants with p value above the genome-wide significance threshold. SNPs are divided into four different categories according to their single-trait and cross-trait characteristics: (i) a “known” shared SNP, defined as a SNP achieved genome-wide significant association with both single traits; (ii) a “single-trait-driven” shared SNP, defined as a SNP achieved genome-wide significant association with one of the two single traits; (iii) a “LD-tagged” shared SNP, defined as a shared SNP in LD with index SNPs previously identified by single-trait GWAS (r² $r^{2}$ ⩾ 0.05), despite not driven by single trait; and (iv) a “novel” shared SNP, defined as a shared SNP neither driven by any single trait nor in LD with index SNPs previously identified by single-trait GWAS (r² $r^{2}$ < 0.05). These four types of SNPs are presented in black, brown, green, and red, respectively, and their rsIDs are listed. Same below. (b) Circular Manhattan plot between eBMD and seropositive RA. (c) Bar plot of significant pleiotropic loci between eBMD and RA. (d) Bar plot of significant pleiotropic loci between eBMD and seropositive RA.

After excluding nine “known” shared SNPs (significantly associated with both traits in previous GWASs), 34 “single-trait-driven” shared SNPs (significantly associated with one of the two traits in previous GWASs), and four “LD-tagged” shared SNPs (in LD (r² ⩾ 0.05) with index SNPs identified by previous single-trait GWAS), we determined two novel shared SNPs for eBMD and RA (Table 2). The most significant novel shared SNP was rs72836346 ( $p_{C P A S S O C}$ = 7.26 × 10⁻¹²), located near BCL2L11, a gene related with T cells apoptosis, which may play a critical role in the treatment of RA.³⁶ The other novel shared SNP was rs2613812 ( $p_{C P A S S O C}$ = 8.22 × 10⁻⁹), located near PTPRS, a gene regulating immune response and chronic inflammation by fibroblast-like synoviocytes.³⁷

Table 2.

Novel SNPs identified in cross-trait meta-analysis between eBMD and RA.

SNP	CHR	BP	A1	A2	Beta		P-eBMD	P-RA	P-CPASSOC	Mapped genes^a	PPH4^b
SNP	CHR	BP	A1	A2	eBMD	RA	P-eBMD	P-RA	P-CPASSOC	Mapped genes^a	PPH4^b
eBMD and RA
rs72836346	2	111,876,613	C	G	−0.020	0.155	8.30 × 10⁻⁰⁸	8.55 × 10⁻⁰⁸	7.26 × 10⁻¹²	BCL2L11, AC096670.3, ACOXL	<0.01
rs2613812	19	5,174,882	G	A	−0.009	0.075	8.40 × 10⁻⁰⁶	5.56 × 10⁻⁰⁶	8.22 × 10⁻⁰⁹	PTPRS, CTC-482H14.5	0.94
eBMD and seropositive RA
rs76458888	6	34,483,749	T	C	0.019	0.127	8.20 × 10⁻⁰⁸	9.92 × 10⁻⁰⁶	3.01 × 10⁻¹⁰	PACSIN1	<0.01

Ensemble variant effect predictor.

A locus was considered colocalized if the PPH4 was greater than 0.8.

A1, effect allele; A2, non-effect allele; Beta, effect allele beta coefficient; BP, physical position of SNP (base-pairs); CHR, chromosome; eBMD, heel estimated bone mineral density; PPH4, posterior probability of H4; RA, rheumatoid arthritis; SNP, single nucleotide polymorphisms.

As for eBMD and seropositive RA, we determined 1 novel shared SNP in addition to 11 “known” shared SNPs, 32 “single-trait-driven” shared SNPs, and 6 “LD-tagged” shared SNPs. The novel shared SNPs rs76458888 ( $p_{C P A S S O C}$ = 3.01 × 10⁻¹⁰) was located near PACSIN1, a gene that functions in immune modulation and inflammation via its regulation of Toll-like receptor 7 (TLR-7) and TLR9-mediated type I interferon (IFN) response in plasmacytoid dendritic cells.³⁸

Among these pleiotropic loci, six loci (rs2016492, rs6908626, rs42032, rs2613812, rs34536443, and rs5754387) shared by eBMD and RA, and four loci (rs6908626, rs42032, rs34536443, and rs5754387) shared by eBMD and seropositive RA were colocalized at the same causal variant (PPH4 > 0.8; Supplemental Table 8). Of note, the novel shared SNP rs2613812 for eBMD and RA showed evidence of colocalization (PPH4 = 0.94).

Transcriptome-wide association study

Through performing TWAS in 49 GTEx tissues, we identified a total of 139 significant tissue–gene pairs, including 23 genes across 41 tissues (Figure 3 and Supplemental Table 9).

Figure 3.

Shared TWAS significant genes between eBMD and RA across 49 GTEx tissues (version 8). Red squares represent shared significant genes for both eBMD–RA and eBMD–seropositive RA trait-pairs; yellow squares represent shared significant genes for eBMD–RA trait-pairs; blue squares represent shared significant genes for eBMD–seropositive RA trait-pairs.

Specifically, 76 significant tissue–gene pairs were found for eBMD and RA, including 19 genes across 40 tissues. Among these tissue–gene pairs, a majority (76.3% = 58/76) remained significant for eBMD and seropositive RA. Of the 19 genes, we found 1 long non-coding RNA (KB-1440D3.14), 13 genes (CCDC116, IRF5, TMEM258, FADS1, YDJC, PEAK1, PRR14, UBE2L3, HLA-J, CD5, CDK6, FAM213A, and AFF3) reported by previous eBMD and/or RA GWAS, and 5 novel genes (PPP1R14B, CCDC189, NEMP2, HIC2, and GTF2IRD2) associated with bone metabolism and/or inflammation response.^39,40 Of note, IRF5, UBE2L3, and CDK6 were also identified by cross-trait meta-analyses.

As for eBMD and seropositive RA, five significant tissue–gene pairs were further identified, including RNF39 in aorta artery, DYDC2 and PEAK1 in basal ganglia, BAD in esophagus mucosa, and CCDC88B in atrial appendage. Of note, PEAK1 was also found by TWAS as shared for eBMD and RA.

Mendelian randomization

We finally performed a bidirectional two-sample MR to test for the causal relationships (Figure 4). Genetically predicted higher level of eBMD significantly decreased the risk of RA (OR = 0.90, 95% CI = 0.84–0.96, p = 6.60 × 10⁻⁴). This causal relationship was further supported by sensitivity analyses and was not affected by horizontal pleiotropy. When we restricted the outcome to seropositive RA, the effect remained consistent (OR = 0.88, 95% CI = 0.82–0.94, p = 2.90 × 10⁻⁴), corroborating the robustness of findings.

Figure 4.

Bidirectional MR analysis between eBMD and RA. The boxes denote point estimate of the causal effects, and the error bars denote the 95% CIs. Inverse-variance weighted was adopted as primary analysis. Weighted median, MR-PRESSO, and MR-Egger were adopted as sensitivity analyses.

On the contrary, neither genetically predicted RA (beta = −0.002, 95% CI = −0.020 to 0.015, p = 0.78) nor genetically predicted seropositive RA (beta = −0.002, 95% CI = −0.017 to 0.014, p = 0.84) seem to affect the level of eBMD.

With the current sample size of exposure, assuming 24.94% (eBMD), 2.97% (RA), and 3.23% (seropositive RA) of phenotypic variance explained by IVs based on the data we used, the mean F-statistics of IVs we calculated were all larger than 10 (ranging from 65.34 to 134.60), indicating strong instruments (Supplemental Table 10).

Discussion

To our knowledge, this is the most comprehensive observational and genetic analyses that systematically investigated the phenotypic and genetic relationships between OP and RA. Our observational analyses found a bidirectional association between OP and RA, which was corroborated by genetic findings. We found a significant negative genetic correlation for eBMD with RA and with its seropositive subtype. Such a shared genetic basis was further substantiated by 99 shared independent loci and 139 shared tissue–gene pairs. Additionally, we found evidence supporting a causal effect of eBMD on RA, but not of RA on eBMD. These findings advance our understanding of the complicated relationship underlying these two common musculoskeletal disorders, and provide important implications for their prevention and treatment.

Our findings are primarily consistent with results from existing studies and greatly extend previous work in several crucial aspects. First, there is a bidirectional association between OP and RA. Although previous case–control studies^4,5 also found a bidirectional association, these studies were prone to control selection bias and recall bias. Our study reinforced previous findings with a substantially improved accuracy by applying a prospective cohort design with a sufficient follow-up period and a larger sample size. Second, eBMD and RA share a genetic basis at both global and local genomic levels. The negative overall genetic correlation we observed (eBMD–RA, = −0.06) was consistent with the result of an earlier analysis (r_g = −0.059).⁷ Different from and extending those findings, we further partitioned the whole genome into LD-independent regions and identified seven significant local signals (1p22.3-p22.2, 2p23.2, 13q32.2-q32.3, 15q23, 15q26.1, 20p12.3, and 21q22.13), of which 15q26.1 remained significant for both eBMD–RA and eBMD–seropositive RA. The gene IQGAP1 on 15q26.1 plays a crucial role in mediating Wnt/β-catenin signaling,⁴¹ a pathway related to the pathogenesis and progression of both OP and RA.⁴² Third, lower eBMD is causally associated with an increased risk of RA. Our MR analyses provided consistent evidence for a causal role of eBMD in RA⁶ and substantially improved the statistical power of prior MR analyses^6–9 by leveraging the most updated GWAS of eBMD and RA in European ancestry and incorporating a larger set of SNPs to improve the proportion of variance explained by IVs. When restricting the outcome to seropositive RA, we found a consistent causal association with even more pronounced magnitude and significance, corroborating the robustness of findings. Furthermore, our bidirectional MR analyses detected neither a causal effect of RA nor its serological subtype on eBMD, suggesting that the observed phenotypic association may be driven by residual confounding or reverse causation. Collectively, our MR findings support BMD decline as a preclinical phenomenon rather than as a clinical outcome in the pathogenesis of RA.

In addition to the significant genetic correlation and causal relationship identified, our cross-trait meta-analyses replicated 96 previously reported loci and discovered three novel pleiotropic loci, suggesting that the observed OP–RA phenotypic link can be largely attributable to biological pleiotropy. We highlight three novel pleiotropic loci (BCL2L11, PTPRS, and PACSIN1). First, BCL2L11 (index SNP rs72836346) encodes the BCL2-like protein 11 (BIM), a protein acting as an apoptotic activator.⁴³ In individuals with RA, T cells overexpress BIM, increasing the susceptibility to apoptosis and accelerating the process of immune aging.³⁶ BIM also mediates the osteoblast apoptosis in the pathogenesis of OP.⁴⁴ Second, PTPRS (index SNP rs2613812) encodes the receptor protein tyrosine phosphatase sigma, reciprocally interacting with chondroitin sulfate or heparan sulfate containing extracellular proteoglycans in the proteoglycan switch. In individuals with RA, fibroblast-like synoviocytes become activated and contribute to joint inflammation and cartilage and bone destruction, while the function of fibroblast-like synoviocytes can be regulated by the proteoglycan switch, which may be a potential target for RA therapy.³⁷ Third, PACSIN1 (index SNP rs76458888), a gene that specifically expressed in human and mouse plasmacytoid dendritic cells, has been linked to human systemic lupus erythematosus by regulating the production of TLR-7 or TLR-9-mediated IFN.³⁸ Several RA mouse models have also found that plasmacytoid dendritic cells and IFN play a role in joint inflammation and bone damage, serving as a potential therapeutic target for individuals with RA.⁴⁵

By integrating data from GWAS and GTEx tissue expression, TWAS suggested shared mechanistic hypotheses between eBMD and RA on a tissue–gene pair level. The three loci (IRF5, UBE2L3, and CDK6) identified in both CPASSOC and TWAS analyses implicate common biological mechanisms in musculoskeletal disorders, involving cell cycle progression,⁴⁶ inflammatory responses,^47,48 and bone remodeling.^49,50 In addition to tissues of skin, skeletal muscle, and blood (well-recognized as relevant to eBMD and/or RA^20,51), our TWAS reported shared regulatory features in the reproductive (e.g., ovary, testis) and digestive systems (e.g., colon transverse), suggesting the possibility of shared pathways extending to a wider range of organs. The underlying mechanisms may partially attribute to sex steroids and gut microbiome. Sex steroids are essential for skeletal development, where estrogen deficiency in postmenopausal women or androgen deficiency in aging men contributes to the development of OP.^52,53 Gut microbiome may affect the pathogenesis of OP and RA by regulating inflammatory responses relevant to the receptor activator of nuclear factor-κB ligand signaling pathway.^54,55 Follow-up experimental studies are needed to validate the role of these hypothesized mechanisms.

Our findings deliver translational implications for public health and clinical practice. First, bone loss as indicated by a decline in BMD or OP may precede the occurrence of RA. While the American College of Rheumatology (2022) only recommends screening for BMD level in individuals with RA initiating or receiving >3 months of glucocorticoid treatment,⁵⁶ we provide novel etiological evidence supporting the routine screening for BMD level in individuals with new-onset RA. Second, OP and RA are inherently linked through biological pleiotropy. The identification of specific pleiotropic variants and pathways regulating common pathological elements may help discover therapeutic targets that would benefit both the prevention and treatment of common musculoskeletal comorbidities. We anticipate that aggregating large-scale GWASs to identify shared genetic underpinnings may guide the development of novel drugs or drug repurposing in the future.

Limitations

We acknowledge several limitations of our study. First, our primary findings were restricted to individuals of European ancestry, constraining the generalizability to other ethnic populations. Second, given the limited power of the hitherto available GWAS on OP, we selected eBMD as a proxy for OP. Despite BMD acknowledged as the best predictor for OP and fracture,^19,21 further GWAS on OP with larger sample sizes is needed to identify the direct effect. Third, we evaluated tissue-expression associations based on 49 tissue types available in GTEx, which may restrict understanding of gene regulatory mechanisms due to limited tissue availability. Although the most relevant tissues and cell types to BMD (bone or bone cells) were not included, the GTEx project demonstrated that many eQTLs can be shared across tissues.⁵⁷ Additionally, the effects in many non-bone tissues and cell types (e.g., skin, skeletal muscle, blood, ovary, testis, and colon transverse) on BMD and bone have been well-known.^20,51

Conclusion

To conclude, our study confirms a putative causal effect of OP on RA through observational and genetic analyses using data from a large population-based prospective cohort and the largest GWAS conducted in European ancestry for eBMD and RA. The OP–RA association is intrinsic, potentially attributed to biological pleiotropy. Findings of the current study clarify genetic etiology underlying the phenotypic link between OP and RA, and provide novel insights into the prevention and the treatment of these two common musculoskeletal diseases.

Supplemental Material

sj-docx-1-tab-10.1177_1759720X251412846 – Supplemental material for Investigating the relationship between osteoporosis and rheumatoid arthritis: prospective cohort study and genetic analyses

Supplemental material, sj-docx-1-tab-10.1177_1759720X251412846 for Investigating the relationship between osteoporosis and rheumatoid arthritis: prospective cohort study and genetic analyses by Lingli Qiu, Wenqiang Zhang, Zhixin Tan, Xuan Wu, Yutong Wang, Mingshuang Tang, Lin Chen, Yanqiu Zou, Yunjie Liu, Bowen Lei, Xiaofeng Ma, Di Zhang, Wenzhi Wang, Yiping Jia, Qiurong He, Lei Sun, Lu Wang, Jian Xu, Yao Chen, Mengyu Fan, Jiayuan Li, Ben Zhang, Xiaoling Wen and Xia Jiang in Therapeutic Advances in Musculoskeletal Disease

Supplemental Material

sj-docx-2-tab-10.1177_1759720X251412846 – Supplemental material for Investigating the relationship between osteoporosis and rheumatoid arthritis: prospective cohort study and genetic analyses

Supplemental material, sj-docx-2-tab-10.1177_1759720X251412846 for Investigating the relationship between osteoporosis and rheumatoid arthritis: prospective cohort study and genetic analyses by Lingli Qiu, Wenqiang Zhang, Zhixin Tan, Xuan Wu, Yutong Wang, Mingshuang Tang, Lin Chen, Yanqiu Zou, Yunjie Liu, Bowen Lei, Xiaofeng Ma, Di Zhang, Wenzhi Wang, Yiping Jia, Qiurong He, Lei Sun, Lu Wang, Jian Xu, Yao Chen, Mengyu Fan, Jiayuan Li, Ben Zhang, Xiaoling Wen and Xia Jiang in Therapeutic Advances in Musculoskeletal Disease

Supplemental Material

sj-docx-3-tab-10.1177_1759720X251412846 – Supplemental material for Investigating the relationship between osteoporosis and rheumatoid arthritis: prospective cohort study and genetic analyses

Supplemental material, sj-docx-3-tab-10.1177_1759720X251412846 for Investigating the relationship between osteoporosis and rheumatoid arthritis: prospective cohort study and genetic analyses by Lingli Qiu, Wenqiang Zhang, Zhixin Tan, Xuan Wu, Yutong Wang, Mingshuang Tang, Lin Chen, Yanqiu Zou, Yunjie Liu, Bowen Lei, Xiaofeng Ma, Di Zhang, Wenzhi Wang, Yiping Jia, Qiurong He, Lei Sun, Lu Wang, Jian Xu, Yao Chen, Mengyu Fan, Jiayuan Li, Ben Zhang, Xiaoling Wen and Xia Jiang in Therapeutic Advances in Musculoskeletal Disease

Supplemental Material

sj-docx-4-tab-10.1177_1759720X251412846 – Supplemental material for Investigating the relationship between osteoporosis and rheumatoid arthritis: prospective cohort study and genetic analyses

Supplemental material, sj-docx-4-tab-10.1177_1759720X251412846 for Investigating the relationship between osteoporosis and rheumatoid arthritis: prospective cohort study and genetic analyses by Lingli Qiu, Wenqiang Zhang, Zhixin Tan, Xuan Wu, Yutong Wang, Mingshuang Tang, Lin Chen, Yanqiu Zou, Yunjie Liu, Bowen Lei, Xiaofeng Ma, Di Zhang, Wenzhi Wang, Yiping Jia, Qiurong He, Lei Sun, Lu Wang, Jian Xu, Yao Chen, Mengyu Fan, Jiayuan Li, Ben Zhang, Xiaoling Wen and Xia Jiang in Therapeutic Advances in Musculoskeletal Disease

Footnotes

Acknowledgements

We are grateful to all investigators who shared genome-wide summary statistics.

Declarations

ORCID iD

Lingli Qiu

Supplemental material

Supplemental material for this article is available online.

References

Vos

Lim

Abbafati

, et al. Global burden of 369 diseases and injuries in 204 countries and territories, 1990–2019: a systematic analysis for the Global Burden of Disease Study 2019. Lancet 2020; 396: 1204–1222.

Black

Cross

Haile

, et al. Global, regional, and national burden of rheumatoid arthritis, 1990–2020, and projections to 2050: a systematic analysis of the Global Burden of Disease Study 2021. Lancet Rheumatol 2023; 5: e594–e610.

Dougados

Soubrier

Antunez

, et al. Prevalence of comorbidities in rheumatoid arthritis and evaluation of their monitoring: results of an international, cross-sectional study (COMORA). Ann Rheum Dis 2014; 73: 62–68.

Meer

Thrastardottir

Wang

, et al. Risk factors for diagnosis of psoriatic arthritis, psoriasis, rheumatoid arthritis, and ankylosing spondylitis: a set of parallel case–control studies. J Rheumatol 2022; 49: 53–59.

Askari

Lotfi

Azimi

, et al. Risk factors of osteoporosis in females: a hospital-based case–control study, Yazd, Iran. Iran J Public Health 2022; 51: 1371–1380.

Liu

Chen

, et al. Rheumatoid arthritis and osteoporosis: a bi-directional Mendelian randomization study. Aging 2021; 13: 14109–14130.

Yang

Cao

, et al. Rheumatoid arthritis and osteoporosis: shared genetic effect, pleiotropy and causality. Hum Mol Genet 2021; 30: 1932–1940.

Kasher

Williams

FMK

Freidin

, et al. Understanding the complex genetic architecture connecting rheumatoid arthritis, osteoporosis and inflammation: discovering causal pathways. Hum Mol Genet 2022; 31: 2810–2819.

Deng

Wong

MCS

. Association between rheumatoid arthritis and osteoporosis in Japanese populations: a Mendelian randomization study. Arthritis Rheumatol 2023; 75: 1334–1343.

10.

, et al. A powerful approach to estimating annotation-stratified genetic covariance via GWAS summary statistics. Am J Hum Genet 2017; 101: 939–964.

11.

England

Smith

Baker

, et al. 2022 American College of Rheumatology Guideline for exercise, rehabilitation, diet, and additional integrative interventions for rheumatoid arthritis. Arthritis Rheumatol 2023; 75: 1299–1311.

12.

Morin

Feldman

Funnell

, et al. Clinical practice guideline for management of osteoporosis and fracture prevention in Canada: 2023 update. CMAJ 2023; 195: E1333–E1348.

13.

Smolen

Aletaha

McInnes

. Rheumatoid arthritis. Lancet 2016; 388: 2023–2038.

14.

Von Elm

Altman

Egger

, et al. Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: guidelines for reporting observational studies. BMJ 2007; 335: 806–808.

15.

Little

Higgins

Ioannidis

, et al. STrengthening the REporting of Genetic Association Studies (STREGA): an extension of the STROBE statement. PLoS Med 2009; 6: e22.

16.

Skrivankova

Richmond

Woolf

BAR

, et al. Strengthening the reporting of observational studies in epidemiology using Mendelian randomisation (STROBE-MR): explanation and elaboration. BMJ 2021; 375: n2233.

17.

Skrivankova

Richmond

Woolf

BAR

, et al. Strengthening the Reporting of Observational Studies in epidemiology using Mendelian randomization: the STROBE-MR statement. JAMA 2021; 326: 1614–1621.

18.

Sudlow

Gallacher

Allen

, et al. UK biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age. PLoS Med 2015; 12: e1001779.

19.

Morris

Kemp

Youlten

, et al. An atlas of genetic influences on osteoporosis in humans and mice. Nat Genet 2019; 51: 258–266.

20.

Al-Barghouthi

Rosenow

, et al. Transcriptome-wide association study and eQTL colocalization identify potentially causal genes responsible for human bone mineral density GWAS associations. Elife 2022; 11: e77285.

21.

Zhu

Bai

Zheng

. Twelve years of GWAS discoveries for osteoporosis and related traits: advances, challenges and applications. Bone Res 2021; 9: 23.

22.

Ishigaki

Sakaue

Terao

, et al. Multi-ancestry genome-wide association analyses identify novel genetic mechanisms in rheumatoid arthritis. Nat Genet 2022; 54: 1640–1651.

23.

Malmström

Catrina

Klareskog

. The immunopathogenesis of seropositive rheumatoid arthritis: from triggering to targeting. Nat Rev Immunol 2017; 17: 60–75.

24.

Rodriguez-Gomez

Gray

, et al. Osteoporosis and its association with cardiovascular disease, respiratory disease, and cancer: findings from the UK biobank prospective cohort study. Mayo Clin Proc 2022; 97: 110–121.

25.

Zhang

, et al. Impact of gallstone disease on the risk of stroke and coronary artery disease: evidence from prospective observational studies and genetic analyses. BMC Med. 2023; 21: 353.

26.

Stensrud

Hernan

. Why test for proportional hazards? JAMA 2020; 323: 1401–1402.

27.

Bulik-Sullivan

Finucane

Anttila

, et al. An atlas of genetic correlations across human diseases and traits. Nat Genet 2015; 47: 1236–1241.

28.

Zhang

, et al. SUPERGNOVA: local genetic correlation analysis reveals heterogeneous etiologic sharing of complex traits. Genome Biol 2021; 22: 262.

29.

Zhu

Feng

Tayo

, et al. Meta-analysis of correlated traits via summary statistics from GWASs with an application in hypertension. Am J Hum Genet 2015; 96: 21–36.

30.

Giambartolomei

Vukcevic

Schadt

, et al. Bayesian test for colocalisation between pairs of genetic association studies using summary statistics. PLOS Genet 2014; 10: e1004383.

31.

Gusev

Shi

, et al. Integrative approaches for large-scale transcriptome-wide association studies. Nat Genet 2016; 48: 245–252.

32.

Burgess

Scott

Timpson

, et al. Using published data in Mendelian randomization: a blueprint for efficient identification of causal risk factors. Eur J Epidemiol 2015; 30: 543–552.

33.

Bowden

Davey Smith

Haycock

, et al. Consistent estimation in Mendelian randomization with some invalid instruments using a weighted median estimator. Genet Epidemiol 2016; 40: 304–314.

34.

Verbanck

Chen

C-Y

Neale

, et al. Detection of widespread horizontal pleiotropy in causal relationships inferred from Mendelian randomization between complex traits and diseases. Nat Genet 2018; 50: 693–698.

35.

Bowden

Davey Smith

Burgess

. Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression. Int J Epidemiol 2015; 44: 512–525.

36.

Shao

Goronzy

Weyand

. DNA-dependent protein kinase catalytic subunit mediates T-cell loss in rheumatoid arthritis. EMBO Mol Med 2010; 2: 415–427.

37.

Doody

Stanford

Sacchetti

, et al. Targeting phosphatase-dependent proteoglycan switch for rheumatoid arthritis therapy. Sci Transl Med 2015; 7: 288ra276.

38.

Xie

Zhou

Athanasopoulos

, et al. De novo PACSIN1 gene variant found in childhood lupus and a role for PACSIN1/TRAF4 complex in Toll-like receptor 7 activation. Arthritis Rheumatol 2023; 75: 1058–1071.

39.

MacArthur

Bowler

Cerezo

, et al. The new NHGRI-EBI catalog of published genome-wide association studies (GWAS Catalog). Nucleic Acids Res 2017; 45: D896–D901.

40.

Kamat

Blackshaw

Young

, et al. PhenoScanner V2: an expanded tool for searching human genotype–phenotype associations. Bioinformatics 2019; 35: 4851–4853.

41.

Carmon

Gong

, et al. RSPO-LGR4 functions via IQGAP1 to potentiate Wnt signaling. Proc Natl Acad Sci U S A 2014; 111: E1221–E1229.

42.

Nusse

Clevers

. Wnt/β-catenin signaling, disease, and emerging therapeutic modalities. Cell 2017; 169: 985–999.

43.

Luo

Rubinsztein

. BCL2L11/BIM: a novel molecular link between autophagy and apoptosis. Autophagy 2013; 9: 104–105.

44.

Liang

Russell

Hulley

. Bim, Bak, and Bax regulate osteoblast survival. J Bone Miner Res 2008; 23: 610–620.

45.

Nehmar

Alsaleh

Voisin

, et al. Therapeutic modulation of plasmacytoid dendritic cells in experimental arthritis. Arthritis Rheumatol 2017; 69: 2124–2135.

46.

Liu

Konagaya

Chung

, et al. Altered G1 signaling order and commitment point in cells proliferating without CDK4/6 activity. Nat Commun 2020; 11: 5305.

47.

Mishra

Crespo-Puig

McCarthy

, et al. IL-1β turnover by the UBE2L3 ubiquitin conjugating enzyme and HECT E3 ligases limits inflammation. Nat Commun 2023; 14: 4385.

48.

Lawrence

Natoli

. Transcriptional regulation of macrophage polarization: enabling diversity with identity. Nat Rev Immunol 2011; 11: 750–761.

49.

Dawidowicz

Allanore

Guedj

, et al. The interferon regulatory factor 5 gene confers susceptibility to rheumatoid arthritis and influences its erosive phenotype. Ann Rheum Dis 2011; 70: 117–121.

50.

Llorente

García-Castañeda

Valero

, et al. Osteoporosis in rheumatoid arthritis: dangerous liaisons. Front Med 2020; 7: 601618.

51.

Wang

Yin

, et al. Novel insight into the aetiology of rheumatoid arthritis gained by a cross-tissue transcriptome-wide association study. RMD Open 2022; 8: e002529.

52.

Compston

. Sex steroids and bone. Physiol Rev 2001; 81: 419–447.

53.

Almeida

Laurent

Dubois

, et al. Estrogens and androgens in skeletal physiology and pathophysiology. Physiol Rev 2017; 97: 135–187.

54.

Mao

Zhou

, et al. Gut microbiome and osteoporosis: a review. Bone Joint Res 2020; 9: 524–530.

55.

Zhao

Wei

Zhu

, et al. Gut microbiota and rheumatoid arthritis: from pathogenesis to novel therapeutic opportunities. Front Immunol 2022; 13: 1007165.

56.

Humphrey

Russell

Danila

, et al. 2022 American College of Rheumatology guideline for the prevention and treatment of glucocorticoid-induced osteoporosis. Arthritis Rheumatol 2023; 75: 2088–2102.

57.

Aguet

Brown

Castel

, et al. Genetic effects on gene expression across human tissues. Nature 2017; 550: 204–213.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.75 MB

0.03 MB

0.04 MB