Inflammation and bone mineral density: A Mendelian randomization study

Osteoporosis is a common age-related disorder leading to an increase in osteoporotic fractures and resulting in significant suffering and disability. Inflammation may contribute to osteoporosis, as it does to many other chronic diseases. We examined whether inflammation is etiologically relevant to osteoporosis, assessed from bone mineral density (BMD), as a new potential target of intervention, or whether it is a symptom/biomarker of osteoporosis. We obtained genetic predictors of inflammatory markers from genome-wide association studies and applied them to a large genome wide association study of BMD. Using two-sample Mendelian randomization, we obtained unconfounded estimates of the effect of high-sensitivity C-reactive protein (hsCRP) on BMD at the forearm, femoral neck,﻿ an﻿d lumbar spine. After removing potentially pleiotropic single nucleotide polymorphisms (SNPs) possibly acting via obesity-related traits, hsCRP, based on 16 SNPs from genes including CRP, was not associated with BMD. A causal relation of hsCRP with lower BMD was not evident in this study.

Osteoporosis increases the risk of fracture. Bone mineral density (BMD), a biomarker of osteoporosis, is a risk factor for fracture which contributes to the global burden of disease from falls 1 . Between 1990 and 2010, global deaths and disability-adjusted life years attributable to low BMD increased, probably due to the aging global population 1 . The burden of osteoporosis (BMD 2.5 standard deviations (SDs) below the reference, i.e. the mean in healthy young adult women 2 ) and osteopenia (BMD 1.0-2.5 SDs below the reference) 2 among older adults is projected to rise in the coming decade, for example in the United States and Australia 3,4 . Apart from the effects on population health, this escalating global burden of osteoporotic disorders places an increasing burden on healthcare providers [4][5][6] , because falls and fractures in older people can generate long hospital stays, extensive convalescence and permanent disability. Inflammation plays an important role in many major non-communicable diseases [7][8][9] , and is increasingly realized to be relevant to the bone remodeling 10 that precedes osteoporosis 11 .
In a recent epidemiological study, older men (aged 65 years or above) with three or more inflammatory markers in the top quartile had a higher risk of hip or vertebral fracture, but the specific associations of each inflammatory marker with fracture at different skeletal sites varied 12 . Observational studies have reported associations of several major inflammatory markers, such as C-reactive protein (CRP), interleukin-6 (IL-6) and monocyte chemotactic protein-1 (MCP-1), with osteoporosis and/or osteopenia [13][14][15] , but findings have not always been consistent 16,17 . Moreover, distinguishing between biomarkers and causal factors is difficult in observational studies, making it unclear whether inflammation is etiologically relevant to osteoporosis as a new potential target of intervention, is a symptom or is a correlate of other unmeasured causal factors. The controversy over current precautionary treatments for low BMD, such as calcium 18 , demonstrates the importance of thorough validation of targets of intervention, as well as the importance of the discovery of new targets.
In this situation assessing whether people with genetically higher levels of inflammatory markers have lower BMD, i.e., Mendelian randomization (MR), provides a way forward. Given genotypes are randomly assorted at conception, MR, or instrumental variable analysis with genetic variants, can mimic a randomized controlled trial and is becoming an important tool for assessing causal relations in observational studies 19,20 . To the best of our knowledge, only one previous MR study has assessed the role of inflammation in BMD and fractures. It did not find a causal relation of higher high-sensitivity CRP (hsCRP) with higher fracture risk, suggesting hsCRP may be a marker of unknown confounders that affect fracture risk 21 . We used a two-sample MR study to estimate the unconfounded associations of several inflammatory markers with BMD. Specifically, we obtained genetic predictors of inflammatory markers from genome-wide association studies (GWAS) and applied them to a GWAS of BMD from a large international consortium. Inflammation and obesity have a complex interplay, which a bidirectional MR study clarified as adiposity likely causing higher CRP rather than CRP causing adiposity 22 . As such, we also checked whether the genetic predictors of inflammation also affected adiposity, to exclude potentially pleiotropic effects.

Genetic predictors of inflammatory markers. SNPs predicting inflammatory markers, including
hsCRP, IL-6, erythrocyte sedimentation rate (ESR) and MCP-1, were obtained from SNPs of genome-wide significance (p < 5 × 10 −8 ) from the largest available published GWAS. Correlations between these SNPs (linkage disequilibrium (LD)) were obtained from SNP Annotation and Proxy Search (SNAP) (http://www.broadinstitute. org/mpg/snap/ldsearchpw.php) using the 1000 genomes reference panel. Highly correlated SNPs (r 2 ≥ 0.8) for each exposure were discarded to retain SNPs with the smallest p-values. A SNP highly correlated with the original SNP was used as proxy when the original SNP was not available in GEFOS. When the remaining SNPs for an exposure were correlated an LD matrix was built and included in the MR analysis to account for that correlation. MR studies assume that the SNPs (instrumental variables) are associated with the outcome only via the exposure 23 , so we performed sensitivity analysis excluding SNPs with potentially pleiotropic effects. Ensembl (http:// www.ensembl.org) integrates genome annotation with available biological data, and provides the phenotypes associated with a given SNP. We manually searched Ensembl for the phenotypes associated with each SNP. SNPs with phenotypes that may affect BMD other than via the relevant exposure were considered as having potentially pleiotropic effects.

Genetic predictors of BMD. Genetic associations with BMD were contributed by the GEnetic Factors for
OSteoporosis (GEFOS) Consortium (http://www.gefos.org). The GEFOS-seq project investigated the associations of SNP with BMD in a general population of European descent (N = 32,965) 24 . Genetic associations with BMD for SNPs with a minor allele frequency ≥ 0.5% in an additive model adjusted for age, age-squared, sex, and weight 24 were made publically available in 2015 25 . BMD was measured using Dual-energy X-ray Absorptiometry at three skeletal sites, the forearm (distal 1/3 of radius), lumbar spine (L1-4) and femoral neck 24 . BMD within each study was standardized to account for systematic differences between Dual-energy X-ray Absorptiometry machines 24 .
Statistical Analysis. SNPs predicting the exposures (i.e. inflammatory markers) were used as instrumental variables in two-sample MR to obtain unconfounded estimates of the effects of hsCRP, IL-6, ESR, and MCP-1 on BMD (SD units) at the forearm, femoral neck, and lumbar spine. To obtain an overall estimate for each exposure, we combined SNP-specific Wald estimates using inverse variance weighting with the standard error calculated as the ratio of the standard error of SNP on outcome to the estimate of SNP on exposure. We used fixed effects when only three or fewer SNPs were available, and used random effects when four or more SNPs were available. Where SNPs for a particular exposure were correlated, we used generalized weighted linear regression and a matrix of their correlations. In order to assess the robustness of these findings, we used a weighted median (WM) estimate, which is consistent when up to 50% of the SNPs are invalid instrumental variables 26 . We also used MR-Egger to test for potentially pleiotropic effects as it may generate correct estimates even when all SNPs are invalid instruments as long as the assumption of instrument strength independent of direct effect (InSIDE) is satisfied. Average directional pleiotropy across genetic variants was assessed from the p-value of the intercept term from MR-Egger 26 . Analysis using MR-Egger has a lower false positive rate but a higher false negative rate than IVW 27 . We also similarly checked the associations of the SNPs predicting inflammation with obesity-related traits from the largest available GWAS among people of European descent in the Genetic Investigation of ANthropometric Traits (GIANT) consortium databases [28][29][30] , and Early Growth Genetics (EGG) consortium database 31 . SNPs with potentially pleiotropic effects were excluded in an additional analysis. We used a p-value of 0.0125 as significance threshold, because we performed MR analysis for four inflammatory markers, i.e. hsCRP, IL-6, ESR, and MCP-1.
All statistical analyses were conducted using R version 3.3.0 (R Foundation for Statistical Computing, Vienna, Austria) and MR analyses were performed using the MendelianRandomization package. This study only used publicly available summary data and hence no ethical approval from an Institutional Review Board was required.
A GWAS by Naitza et al. provides genetic predictors of IL-6, ESR, and MCP-1 (in SD units) based on 4292, 3596, and 4295 people respectively aged 14-102 years from the Lanusei Valley of Sardinia, Italy 35,36 . Four SNPs from ABO were used to predict IL-6 (3 SNPs in the analysis of each skeletal site). Two proxy SNPs were used, i.e. rs657152 (ABO) as a proxy for rs643434 (ABO) and rs544873 as a proxy for rs687289 (ABO). However, highly pleiotropic effects of ABO cast doubt on whether these SNPs solely affect BMD via IL-6. We repeated the analysis for IL-6 using SNPs from two other sources. Three SNPs in IL6R, predicting natural log-transformed IL-6, were obtained from a previous MR study by the IL-6R MR Analysis Consortium in 4489 Europeans 37 . rs61812598 (IL6R) was used as a proxy for rs7529229 in the analysis of BMD at the femoral neck and lumbar spine. Another three SNPs predicting IL-6 (from KPNB1, SERPINE2, TC2N) were obtained from a GWAS by Ahola-Olli et al. in SCieNTiFiC RepoRts | 7: 8666 | DOI:10.1038/s41598-017-09080-w 1664 Finns 38 . Four SNPs from CR1 were used to predict ESR, of which 3 were included in the analysis of BMD at the femoral neck and lumbar spine respectively. Twenty-three SNPs were used to predict MCP-1 (from ACKR1, AIM2, CADM3, DUSP23, FCER1A, OR10J1, OR10J2P, OR10J7P, OR10J9P, and intergenic variant), among which 22 were included in the analysis of BMD at the femoral neck and lumbar spine respectively. rs79433881 (OR10J1) was used as a proxy for rs11265186 (OR10J9P) in the analysis of BMD at the femoral neck, and rs4290055 was used as a proxy for rs11265186 (OR10J9P) in the analysis of BMD at the lumbar spine.
We matched the chosen SNPs with GWAS of obesity-related traits and found none of the SNPs achieved genome-wide significance in GIANT or EGG ( Table 1 in the Supplementary material). Nevertheless, some of these SNPs (rs1260326 (GCKR), rs13233571 (BCL7B), rs2847281 (PTPN2), and rs4129267 (IL6R) predicting hsCRP and rs630014 (ABO), rs651007 (ABO), and rs687289 (ABO) predicting IL-6) were associated with obesity-related traits ( Table 2 in the Supplementary material). None of the individual SNPs predicting MCP-1 were associated with obesity-related traits, but the IVW estimate showed potential associations of MCP-1 with higher BMI in men and higher risk of obesity class 3 ( Table 3 in the Supplementary material). rs12075 (DARC/CADM3) that predicted MCP-1 also has potential pleiotropic effects on obesity-related traits based on a comprehensive search of Ensembl.
Funnel plots for hsCRP, IL-6, ESR, and MCP-1 show some asymmetry (Figures 1 to 4 in the Supplementary material). Figures 5 to 8 in the Supplementary material show the scatter plots for hsCRP, IL-6, ESR, and MCP-1. We performed analyses including all available SNPs and sensitivity analyses excluding SNPs with potentially pleiotropic effects such as SNPs in the pleiotropic genes GCKR and ABO. Table 1 shows the estimates of the effects on BMD. hsCRP was not associated with BMD in any model. Figure 1 shows the SNP-specific associations and MR estimates for hsCRP. Considering the limitations of the SNPs predicting IL-6, ESR, and MCP-1, estimates for these exposure are only presented in Table 4 and Figures 9 to 11 in the Supplementary material. MCP-1 was associated with lower BMD at the forearm, but with higher BMD at the femoral neck and lumbar spine. However, the association of MCP-1 with BMD at the forearm was no longer evident after removing rs12075 (DARC/CADM3) in a sensitivity analysis ( Table 4 in the Supplementary material). Intercepts from MR-Egger showed little directional pleiotropy, but the direction of effect estimates was not always consistent for both IVW and WM.

Discussion
For the first time, we assessed the associations of four major inflammatory markers (i.e. hsCRP, IL-6, ESR, and MCP-1) with BMD using two-sample MR to provides unconfounded estimates because genotypes are randomly assorted at conception and usually unaffected by the confounders that typically bedevil epidemiological studies, such as socio-economic position and health status. Use of two different samples for the exposures and outcomes also minimizes susceptibility to confounding because of the different underlying data structures in each sample 39 . We only used SNPs strongly associated with inflammatory markers from the largest available GWAS to reduce the risk of false positives, and we used LD matrices to account for any correlation between SNPs predicting each exposure. We also checked for potentially pleiotropic effects and removed them in sensitivity analyses. We did not observe associations of higher hsCRP with lower BMD.
Although we checked the assumptions of MR rigorously, some limitations exist. First, genetic predictors of IL-6, ESR, and MCP-1 were based on an isolated population in the Lanusei Valley of Sardinia, Italy 36 . These results may be susceptible to population stratification and might not be applicable to other populations. In addition, no  given that genetic instruments may only explain a small amount of variation of the phenotype 40 . We were unable to calculate the power for this study because the GWAS of Naitza et al. and of Prins et al. did not report the explanation power of SNPs on inflammatory markers. However, one advantage of two-sample MR is increased statistical power 41 . We also estimated the F-statistics based on the effect allele frequency, estimated SNP effect on exposure, and sample size 42 . F-statistics range from 32 to 4419, except for the three SNPs from Ahola-Olli et al., which had F-statistics ranging from 3 to 8 ( Table 5 in the Supplementary material). An F-statistic higher than the conventional threshold of 10 is usually taken as implying little weak instrument bias. Fifth, associations of higher genetically predicted MCP-1 with lower BMD were only observed at the femoral neck and lumbar spine. Given that BMD is likely to be correlated at different skeletal sites, this may suggest chance findings. Sixth, local and systemic inflammation may impact bone growth via different pathways 43 , which we could not distinguish. Seventh, GEFOS is based on the general population, but osteoporosis is more prevalent among the elderly. However, bone loss may start in middle age and increase the risk of osteoporosis and fracture in later life. Lastly, while the prevalence of osteoporosis is higher in women, sex-specific analysis was not possible in this study, although causal relations are usually consistent. Previous observational studies have reported inconsistent associations of inflammatory markers with BMD 13,16,17,44 . Inflammation may affect BMD via the bone remodeling cycle regulated by osteoblasts and osteoclasts 45 affecting osteocalcin 46 . Previous studies have observed inverse associations of hsCRP with serum osteocalcin 47,48 , which may play a role in osteoclastogenesis 49 . However, an in vitro study suggested that CRP inhibits both osteoblast and osteoclast differentiation 50 . Hence, the impact of these inflammatory markers on BMD is unclear. SNPs predicting hsCRP with potentially pleiotropic effects via obesity-related traits were identified in our study, which may imply that the previously reported associations of inflammation with BMD are not causal. Also, higher follicle-stimulating hormone is associated with both bone loss 51 and higher levels of inflammatory markers including CRP and MCP-1 52 , implying a possible confounder of the association of inflammation with bone loss.
A causal relationship of hsCRP with lower BMD was not evident in this study. A replication in the UK Biobank would help clarify whether these null findings are due to lack of power or indicate no effect and that other means of preventing osteoporosis should be sought.