The analysis of causal relationships between blood lipid levels and BMD

Bone mineral density (BMD) and lipid levels are two of the most extensively studied risk factors for common diseases of aging, such as cardiovascular disease (CVD) and osteoporosis (OP). These two risk factors are also correlated with each other, but little is known about the molecular mechanisms behind this correlation. Recent studies revealed that circulating levels of several metabolites involved in the biosynthesis of androsterone correlate significantly with BMD and have the capacity to affect cholesterol and lipids levels. A main aim of the present study was to investigate the hypothesis that androsterone-related metabolites could provide a link between CVD and OP, as a common cause of lipid levels and BMD. The present study employed data from the NIHR BRC TwinsUK BioResource, comprising 1909 and 1994 monozygotic and dizygotic twin pairs, respectively, to address the causal relationships among BMD and lipids, and their associated metabolites, using reciprocal causation twin modelling, as well as Mendelian randomization (MR) using large publicly-available GWAS datasets on lipids and BMD, in conjunction with TwinsUK metabolite data. While results involving the twin modelling and MR analyses with metabolites were unable to establish a causal link between metabolite levels and either lipids or BMD, MR analyses of BMD and lipids suggest that lipid levels have a causal impact on BMD, which is consistent with findings from clinical trials of lipid-lowering drugs, which have also increased BMD.


Introduction
Bone mineral density (BMD) and serum total cholesterol (TC) are among two of the most extensively studied clinically-oriented phenotypes, associated with the two of the most common polygenic age-related pathologies, osteoporosis (OP) and cardiovascular disease (CVD). There is an abundant literature showing correlation between these biomarkers and possible mechanisms underlying it [1][2][3][4], although the views and the data remain controversial [5,6].
In particular, it is unclear whether BMD and lipids are related due to pleiotropic genetic influences or through direct causal influence of one on the other. One of the most convincing studies employed conditional false-discovery rate and other approaches on large genome-wide a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Lipid measures in the TwinsUK sample included four traits: total cholesterol (TC), highdensity lipoprotein cholesterol (HDL-C), triglycerides (TG) and low-density lipoprotein cholesterol (LDL-C) calculated using the formula: LDL-C = ¾ (TC-HDL-C) [17]. Subjects were followed up and assessed at multiple time points, over a period of up to 17 years.
DEXA was used to measure BMD at the lumbar spine (L1 to L4) and hip regions (following manufacturer's recommendations (QDR 4500W system, Hologic Inc, Bedford, MA) and described elsewhere [18]. In total, 15,491 DEXA scans were performed in 7056 twins during 17 years of follow-up; control scans were performed using the spine phantom.

Data analysis
The overall analysis strategy, along with data employed, is presented in Fig 1. All the analyses were carried out in the R statistical environment. DEXA scans and lipid profiles were obtained at multiple ages, therefore we used the lme4 package (https://cran.rproject.org/web/packages/lme4/lme4.pdf) to perform linear mixed-effects modelling to produce a single phenotype per individual by regressing out the age effect and using the intercept from those regressions as the individual score for each phenotype. This has the advantage of using all available data while adjusting for age.
Both standard univariate and multivariate twin modelling was performed using the umx package [19] to estimate the additive genetic (A) and shared in common by a twin pair (C) and nonshared (E) environmental variances and covariances. Direction of causation in twin modelling was performed using the OpenMX package [20][21][22]. We fitted the model illustrated in Fig 2, depicted for a single member of a twin pair. Models were fitted to the raw data and fit functions minimized with the CSOLNP optimizer [23]. The model contained three metabolites, two lipid measures, and two BMD measures (lumbar spine and left hip), with the omission of one metabolite and two other lipid measures as explained at the beginning of the results section. The model had a single latent factor for the metabolites, a latent factor for the lipids, and a latent factor for BMD. Parameters were fitted to allow for each latent factor to causally influence other latent factors in both possible directions (6 reciprocal causation parameters). The variances of each of the three latent factors were partitioned into genetic, shared environmental, and unique environmental components (9 parameters). To identify these three latent factors, the first loading on the observed variables for each factor was fixed to unity, leaving two free parameters loading on the metabolites, and one each on TG and left hip BMD (4 measurement model parameters). Finally, each observed variable had genetic, shared environmental, and non-shared environmental phenotype-specific variances estimated (3 components x 7 phenotypes = 21 parameters). All models equated means to be equal across twin pairs and zygosity (7 mean parameters). The total number of parameters estimated in the full model was therefore 47. Variables were re-scaled to have approximate unit variance to aid in optimization. Results of a standardized solution are presented, with all observed and latent variables standardized to unit variance.
Mendelian randomization analyses were performed using the MR-Base platform [24]. Four independent SNPs were chosen as instruments for the metabolites of interest. For ATS, the two independent SNPs which were most significantly associated in TwinsUK (rs7809615 on chromosome 7 and rs182420 on chromosome 19, significantly associated at p<10 −106 and 2x10 -8 , respectively) were selected as instruments. These two SNPs also had large and significant effects on DHEA-S and rs7809615 had a significant effect on EAS, but no other independent SNPs were associated with these two metabolites, so no instruments were available for them. Additionally, the two independent SNPs most associated with Δ 4 -dione in TwinsUK (rs2762353 on chromosome 6 and rs4149056 on chromosome 12, significantly associated at p<3x10 -9 and 2x10 -11 , respectively) were also selected as instruments. The large meta-GWAS of lipids [25] and GWAS for BMD established by heel ultrasound in the UK Biobank [26] were used as outcomes for this series of MR analyses (described below and in Fig 1). We also performed bidirectional MR with SNPs for lipids as instruments and BMD as outcome and the

Twin modelling
Pearson correlations between phenotypes of interest (4 metabolites, 4 lipid measures, and 2 BMD measures) were estimated (Table 1). ATS and AES were highly correlated, as were LDL-C and TC. As expected, HDL-C was negatively correlated with the other lipid measures, but also with the BMD measures. TG was positively correlated with BMD. TC showed no correlation with BMD, while LDL-C was weakly correlated with one of the BMD measures.
All metabolites were significantly correlated with TG and TC, while only Δ 4 -dione was correlated with HDL-C and only DHEA-S and EAS were correlated with LDL-C. To avoid collinearity due to strong correlation of ATS with EAS (r = 0.94) and LDL-C with TC (r = 0.93), we omitted EAS and LDL-C from the standard multivariate genetic model (8 phenotypes analyzed). Also, we used only TG and HDL-C for the reciprocal causation modelling (7 phenotypes analyzed).
A multivariate (multiphenotype) full genetic model was first fitted to the twin data to explore the underlying genetic and environmental architecture. For estimating the genetic and environmental covariance matrices among the 8 measures, 36 genetic, 36 shared environmental, and 36 nonshared environmental variances and covariances were directly estimated (108 parameters). In addition, means were estimated for each variable, equated across twin pairs and zygosity (8 parameters). Tables 2-4 present the genetic, shared environmental, and nonshared environmental phenotypically-standardized covariance matrices, respectively. Summing these three covariance matrices yields the expected phenotypic correlation matrix among the 8 measures. All phenotypes were highly heritable, ranging from .50 for TC to .82 for lumbar spine BMD (diagonal elements in Table 2), and all the heritability estimates were statistically significant (p < .001). Shared environmental influences were generally low, ranging from -.01 for triglycerides to .21 for DHEA-S (diagonal elements in Table 3). While negative variances make no theoretical sense, in the standard twin model, they do imply that the DZ correlation is less than half the MZ correlation, which would suggest dominance variance rather than shared environmental variance, but this negative parameter was not statistically significant. Finally, nonshared environmental variance (diagonal elements in Table 4), while small, is somewhat higher for the metabolite measures and lipids, on average, than for BMD.
Examining the sources of variance underlying the correlations among measures, first genetic variation within domains, all metabolites were significantly genetically inter-correlated (off-diagonal elements in Table 2). All but one (TC with HDL-C) genetic covariance among lipids was significant, and the genetic covariance between BMD measures was substantial and significant. Examining the genetic components of cross-domain correlation, of all the genetic covariances between metabolites and lipids, only ATS with HDL-C was significant, though modest (.07). Five of 6 genetic covariances between metabolites and BMD measures were significant. Finally, of the lipids, only HDL-C was significantly genetically correlated with both BMD measures, while TG was correlated with hip BMD.
Shared environment also contributed significantly to all correlations between the metabolites (off-diagonal elements in Table 3), but did not significantly contribute to BMD variance or covariance, nor to lipid covariances. Some of the covariation between metabolites and lipids was due to shared environmental influences, with 4 out of 9 covariances statistically significant, although all rather small. Finally, the shared environment did not contribute to the correlations between metabolites and BMD measures.
Nonshared environment contributes substantially to the covariances between metabolites, likely reflecting day-to-day co-variation of metabolites involved in the same biochemical pathway (off-diagonal elements in Table 4). Nonshared environment also contributed significantly to the covariances among lipids and among BMD measures, though the magnitude of the covariances among lipids were smaller than for BMD. Nonshared environment generally did not contribute to correlations between metabolites and lipids, between metabolites and BMD, nor between lipids and BMD, though a few of the covariances were statistically significant.
Next, we fitted a direction of causation twin model to the data, with common factors underlying the covariation among metabolites, lipids, and BMD measures, with standardized parameter estimates from the full model presented in Fig 2. We tested all parameters (with the exception of measure-specific E, which contains measurement error) against the full model (Table 5). We first tested the common factor A, C, and E variances. All three common factors (metabolites, lipids, and BMD) had significant genetic and nonshared environmental variance components, but only the metabolites had a significant shared family environmental component. Looking at the measure-specific variances, DHEA-S, TG, and hip BMD did not have specific genetic variance components (the genetic common factor was sufficient to explain the genetic variance present), but the other four measures did. Only HDL-C had a significant measure-specific shared environmental component.
To test the direction of causation, each causal path was dropped and tested against the full model that allows for reciprocal causation. No evidence for causal relationships was found after Bonferroni correction (p = 0.05/6 = .008 for 6 possible causal paths) except for the causal path from BMD to the metabolite common factor. Looking at the underlying correlations among latent factors, ultimately they were likely too low to establish direction of causation, though not too low to establish correlation. As can be seen when testing both causal directions simultaneously, the BMD common factor was found to be phenotypically correlated with just the metabolites, with the other two tests being not statistically significant.

Mendelian randomization
We used Mendelian randomization to test the hypothesis that changes in the levels of ATS and Δ 4 -dione would cause changes in the lipid and BMD levels ( Table 6). We present single SNP tests as well as inverse variance weighted tests of the two SNPs together. Reverse causation couldn't be explored due to small sample size of the metabolite dataset. We employed a conservative Bonferroni multiple testing correction for 42 tests (2 metabolites x 3 instruments (two SNPs and joint test) x 7 traits) giving a p-value threshold of < .0012. We found evidence that ATS levels causally influence TG levels for one of the SNP instruments only (p < 0.0001). There is also suggestive evidence for ATS levels influencing TC and right and total BMD from one SNP instrument only. There is also suggestive evidence that Δ 4 -dione  We also explored the question of whether lipid levels cause BMD, BMD causes lipid levels, or both (Table 7). We used lipid level and BMD SNPs from the same public datasets as instruments. We calculated several MR test statistics, though present results for the new weighted mode estimator, which has been shown to be robust even when the majority of instruments violate assumptions of MR [27].
There is evidence of LDL-C levels causing BMD. Even with a conservative Bonferroni correction for 8 tests, 4 lipid levels causing total BMD and total BMD causing 4 lipid measures (threshold of .05/8 = .0063), LDL-C significantly causes total heel BMD, with p-values also nominally significant for right and left heel measures separately. There was no evidence of directional horizontal pleiotropy for LDL-C, with the MR-Egger intercepts found not to be statistically significant. In addition, with the exception of HDL exposure on BMD, all other MR-Egger intercepts were also not nominally statistically significant. All tests of trait heterogeneity, however, were statistically significant. Fig 3 contains plots of the MR analysis of the LDL-C exposure on total heel BMD outcome. Panel a presents the forest plot of individual SNP effects, along with 95% confidence intervals. Its symmetry suggests no pleiotropic effects. While the MR-Egger regression wasn't statistically significant, the inverse variance weighted (IVW) analysis did yield an overall significant effect of LDL-C. Panel b compares the slopes from the various MR methods employed and shows the weighted mode estimator that we presented above has a slope very similar to the other methods employed. Panel c presents the funnel plot showing the relationship between the causal effect of LDL-C on BMD estimated by each SNP against the inverse of the standard  Causal relationships between blood lipid levels and BMD

Discussion
There are many studies on subjects of different ethnic backgrounds that have demonstrated a significant association between serum lipid profile, specifically TC and LDL-C, and BMD, both in the general population and in alcoholic people [28]. These observations raise the question of possible underlying mechanisms that explain the association between serum lipids and bone metabolism. One such mechanism could be related to metabolism of endogenous steroid hormones. TC is a major structural component of membranes and a substrate for the biosynthesis of other steroids, including ATS and DHEA-S [9]. On the other hand, the important role of endogenous and exogenous steroid hormones in osteoporosis etiology is well known.
Recently, a significant correlation of the metabolites selected in this study with hip and spine BMD was reported [8]. The major aim of the present study was therefore twofold: 1) to test the hypothesis that three physiological facets (lipids, BMD, and steroid hormones-related metabolites) have a common genetic background, and 2) to examine possible causal relationships. Using MR, we found some evidence of ATS levels influencing TG and TC, as well as BMD.
Also, it appears Δ 4 -dione may be causing BMD and TG levels. Given lack of statistical significance for most of these association, or significance for one SNP instrument and not the other, the causal relationship between these metabolites and lipids and BMD has not been convincingly demonstrated. In addition, our results are not completely consistent with a recent analysis of TwinsUK data, which found evidence for both these metabolite levels influencing BMD, and this result was also replicated in a Hong Kong sample [8]. However, our approach used BMD outcome measures from the UK Biobank, rather than from the same individuals on whom metabolites were assayed. Furthermore, with only two SNP instruments available for each of the metabolites of interest, the study is not optimally powered nor can we check for violations of MR assumptions, as can be done when employing multiple instruments. A number of studies looked at the correlation between serum HDL and BMD [29], but no consistent relationship was found, with some studies finding positive correlations and others negative correlations. Addressing issues of causation, one study found lower TC and LDL levels, along with an expected higher BMD in postmenopausal women taking hormone replacement therapy (HRT) compared to those who were not [30]. This suggests either BMD levels influence cholesterol levels or that there is a common cause of both BMD and cholesterol, such as estrogen levels. Our MR results showed no evidence of BMD causing lipid levels. In a recent meta-analysis, statin use was linked to increased BMD, though not to reduced fracture [31], suggesting a causal relationship from cholesterol to BMD, or a common cause of the two. Consistent with this, our MR results demonstrated evidence of LDL-C levels influencing BMD.
Direction of causation modelling in twin studies has been implemented successfully in a number of medical and behavioral phenotypes [11,[32][33][34][35][36]. While the approach relies on very different data and methods, it can address the same questions of causation as RCT or MR studies. The major limiting factor in the twin study approach to causation is that it requires the traits of interest to differ in heritability and shared environment, otherwise the approach lacks power. Our twin modeling demonstrated significant genetic overlap among the analyzed traits, but reciprocal causation modeling of twin data did not yield definitive results. This is likely explained by the little difference in the heritability estimates between the most heritable traits (BMD) and the least heritable traits (metabolites). Shared environmental variance was also minimal in all the measures examined, further limiting the power of our study design. A further limitation of the twin modelling is that because TwinsUK is a community-based sample, detailed clinically-relevant information is unavailable, including such important covariates as use of lipid-lowering or BMD-enhancing drugs, or use of environmental interventions, potentially further attenuating power.
The animal model literature suggests shared genes predispose to both HDL-C and BMD [29]. Our twin analyses also found a significant genetic correlation between HDL-C and BMD. However, using LD-score regression and results from the same large public datasets that we analyzed [37], no significant genetic correlation between HDL levels and BMD levels was shown, while the genetic correlation was -.14 for LDL and lumbar spine BMD (p = .0376) and -.145 for LDL with neck BMD (p = .0109). There was also a genetic correlation of -.137 for TC and lumbar spine BMD (p = .0135) and -.138 for TC with neck BMD (p = .0028), but no statistically significant genetic correlation for TG levels. Our twin analyses failed to find a significant genetic correlation between TC (highly correlated with LDL-C) and BMD, however.
In summary, our study confirmed previous findings on the relationship between lipid levels and BMD, as well as their correlation with the metabolites we examined. The twin modelling failed to establish the direction of causation. However, the MR analyses suggest that causation goes from LDL-C level to BMD, consistent with previous studies. Further studies, particularly of the metabolic factors, would benefit from a larger sample of twins both for reciprocal causation modelling and also to test BMD and lipid genetic variants effects on metabolite levels.