Metabolomic Pathways to Osteoporosis in Middle‐Aged Women: A Genome‐Metabolome‐Wide Mendelian Randomization Study

ABSTRACT The metabolic state of the body can be a major determinant of bone health. We used a Mendelian randomization approach to identify metabolites causally associated with bone mass to better understand the biological mechanisms of osteoporosis. We tested bone phenotypes (femoral neck, total hip, and lumbar spine bone mineral density [BMD]) for association with 280 fasting blood metabolites in 6055 women from TwinsUK cohort with genomewide genotyping scans. Causal associations between metabolites and bone phenotypes were further assessed in a bidirectional Mendelian randomization study using genetic markers/scores as instrumental variables. Significant associations were replicated in 624 participants from the Hong Kong Osteoporosis Study (HKOS). Fifteen metabolites showed direct associations with bone phenotypes after adjusting for covariates and multiple testing. Using genetic instruments, four of these metabolites were found to be causally associated with hip or spine BMD. These included androsterone sulfate, epiandrosterone sulfate, 5alpha‐androstan‐3beta17beta‐diol disulfate (encoded by CYP3A5), and 4‐androsten‐3beta17beta‐diol disulfate (encoded by SULT2A1). In the HKOS population, all four metabolites showed significant associations with hip and spine BMD in the expected directions. No causal reverse association between BMD and any of the metabolites were found. In the first metabolome‐genomewide Mendelian randomization study of human bone mineral density, we identified four novel biomarkers causally associated with BMD. Our findings reveal novel biological pathways involved in the pathogenesis of osteoporosis. © 2017 American Society for Bone and Mineral Research.


Introduction
I t is now established that genetic factors, environmental factors, and their interactions play a major role in osteoporosis pathogenesis. (1) In addition, metabolic pathways play an important role in age-related bone loss. Several circulating proteins and metabolites are known biomarkers of bone turnover and can help to identify the state of bone formation/resorption for diagnostic or therapeutic purposes. (2) Observation of bone loss accompanying several endocrinological, inflammatory, gastrointestinal, renal, and nutritional disorders, (3) often known as secondary osteoporosis, also suggests that the metabolic state of the body can be a major determinant of bone health. Identification of metabolites that are causally associated with This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. reduced bone mass will increase our power to understand the biological mechanisms of osteoporosis and to develop drug targets for bone health.
Causality, however, cannot be inferred from observational studies, (4) and randomized controlled trials are considered the gold standard for establishing causal relationships. Recently, Mendelian randomization (MR), an instrumental variable analysis method, has been proposed as an alternative for establishing causal relationships. (5) For inference of causal relationships between an exposure (eg, a metabolite) and an outcome (eg, bone mineral density [BMD]), a susceptibility variant (or the genetic risk score) of the exposure is used as instrument. The association of this instrument with the exposure is not affected by any confounding factors (because the genetic sequence is randomly assigned during conception), and it is associated with the outcome only via its effects on the exposure. (6) MR has been successful in inferring causality for several clinically important traits (7) and for elucidating the complex etiology of common disease. (8) In this study, we used MR on a genome-metabolome-wide scale. We first looked for circulating metabolites associated with bone density measures in a large cohort of twins with metabolomic data. We then searched for the genetic markers of those metabolites associated with BMD and used them as instrumental variables to identify metabolites causally associated with bone health. For evaluating bidirectional effects, we also used genetic risk scores for BMD measures as instrumental variables. Lastly, we replicated our findings in an independent cohort from Hong Kong.

Materials and Methods
The flowchart of the study design is presented in Fig. 1.

Discovery cohort
The TwinsUK study, started in 1992, is a nationwide registry of volunteer twins in the UK, with about 12,000 registered twins (83% female predominantly middle-aged and older). (9) The cohort has extensive demographic, physiological, behavioral, and lifestyle data available and is one of the most deeply phenotyped and genotyped cohorts in the world. The study has been an active partner of many large collaborative projects including Genetic Factors for Osteoporosis (GEFOS) consortium. (10) The study has been approved by St. Thomas' Hospital Research Ethics Committee, and all twins have provided informed written consent.
Dual-energy X-ray absorptiometry (DXA) was used to measure BMD at the lumbar spine (L 1 to L 4 ) and hip regions (Hologic Discovery W devices, Hologic, Bedford, MA, USA). In total, 15,491 hip and spine DXA scans were performed in 7056 twins during 17 years of follow-up. (11) All BMD measurements were performed by trained technicians using a standardized protocol of measurement. Daily quality-control scans were performed using the spine phantom. For this study, we included all female twins with metabolomic data and at least one BMD measurement in hip/spine regions. DXA scans from the same date of blood sampling for metabolomic analysis were used in 4650 twins, and BMD measurements from the closest date (if multiple measurements) were selected for 453 twins without DXA in the same date of blood sampling (up to 2 years differences were allowed). This resulted in a study population of 5103 participants (mean age 53.8 AE 13.8 years).

Replication cohort
The Hong Kong Osteoporosis Study (HKOS) is a prospective longitudinal cohort study of osteoporosis and fracture that started in 1995. The cohort participants were communitydwelling Southern Chinese men and women of Han descent recruited from public road shows and health fairs held in various districts of Hong Kong from 1995 to 2010. (12) A total of 9449 participants were recruited. In 2015, a full-scale in-person followup study was conducted. DXA has been used to measure BMD at the lumbar spine (L 1 to L 4 ) and hip regions (Hologic QDR 2000plus and 4500plus systems at baseline; Hologic Discovery A and Horizon A at follow-up study). All BMD measurements were performed by trained technicians using a standardized protocol. Daily quality-control scans were performed using the spine phantom and across the DXA machines to make sure of consistent measurements.

Metabolomic profiling
Non-targeted metabolite detection and quantification was conducted by the metabolomics provider (Metabolon, Inc., Durham, NC, USA) on fasting blood samples of 6055 twins, as previously described. (13) This platform incorporates two separate ultra-high-performance liquid chromatography/tandem mass spectrometry injections and one gas chromatography/ mass spectrometry injection per sample, and can detect metabolites in the range of low nanograms per milliliter. (14) In previous studies, the platform findings have been shown to be highly stable (15) and reproducible (with median relative standard deviation of 5% for internal standards added to samples before mass spectrometry). (14) The same metabolomic profiling was also performed on 340 and 291 serum samples from different participants in the HKOS baseline and follow-up studies, respectively. Blood samples were taken in the same day of DXA scans for all patients.
In this study, we analyzed 280 structurally named biochemicals (known metabolites) categorized into the following broad categories: amino acids, lipids, carbohydrates, vitamins, nucleotides, peptides, xenobiotics, and steroids.

Genomewide association study (GWAS)
Genomewide genotyping was performed for 5710 participants of the TwinsUK study using two chips (HumanHap300 BeadChip and HumanHap610 QuadChip, Illumina, San Diego, CA, USA). The genotype data have been imputed (using HapMap II reference panel containing $2.5 million autosomal singlenucleotide polymorphisms [SNPs]) and used in several international consortia for different phenotypes. (10) In this study, we performed genomewide association studies (GWAS) for all metabolites including those significantly associated with bone phenotypes. (16) The software MERLIN was used to account for twin structure. Results from the largest GWAS meta-analysis for BMD measures (17) were used to calculate genetic risk scores (GRS) for hip and spine BMD.

Statistical analysis
All models were built using Stata (version 14, StataCorp LLC, College Station, TX, USA) and SIMCA software (version 13, Umetrics, Umeå, Sweden). For quality control of metabolomic data, we had to consider that mass spectrometry was performed over several days for the cohort samples, and measured concentrations were variable according to run days. We, therefore, normalized the metabolite data by dividing each metabolite concentration by its median in the respective run day and then inverse normalized the data because the metabolite concentrations were not normally distributed. Moreover, to avoid spurious false-positive associations due to small sample size, we excluded metabolomic traits with more than 20% missing values. We imputed individual metabolites not detected in a sample using the minimum measures in their run day. The same QC procedure was performed in HKOS.

Association analysis in TwinsUK
To assess the predictive power of combined metabolites, we used orthogonal-partial least squares (O-PLS) regression. We used a sevenfold cross-validation technique to validate the O-PLS models. (18) Goodness-of-fit measures were reported for all models in training (coefficient of variation labeled as R 2 ) and cross-validation (coefficient of variation labeled as Q 2 ) sets. For univariate metabolite-BMD trait associations, we used mixedeffects random-intercept models considering family relatedness and zygosity as random effects variables, and age, height, weight, and duration of hormone-replacement therapy (HRT; as a continuous variable in units of years) as confounders.

Mendelian randomization analysis
For the metabolomic GWAS analysis, we only considered metabolites passing the Bonferroni corrected threshold of p < 1.79 Â 10 À4 (0.05 divided by 280, for the number of known metabolites) for association with at least one of the BMD phenotypes. Genomic loci passing a stringent threshold for association with these metabolites (p < 1 Â 10 À10 for the sentinel SNP in a locus) as well as hip and spine GRS scores were considered for the bidirectional instrumental variable analysis. (16) In the direct metabolite-to-BMD analysis, metabolite levels were used as endogenous variables; age, height, weight, and duration of HRT were included as exogenous variables; and allele frequency (additive model) was used as the instrument. Bidirectional Mendelian randomization (MR) (19,20) also investigates for any evidence of reverse causality (ie, BMD causally influencing metabolite levels) by use of GRS scores as the instrument and BMD measures as the endogenous variables for metabolite level outcomes. We used the generalized method of moments (GMM) with cluster-robust heteroskedastic-consistent variance estimates. (21) Under-, weak-, and over-identification limitations were checked. Bootstrap resampling of the data (10,000 repetitions) was used for the calculation of robust estimates of standard errors from instrumental variables analysis.

Replication in the HKOS
Multivariable robust regression models were used to evaluate metabolites found to be causally associated with BMD in the TwinsUK cohort. We reported all the association between these metabolites and BMD, adjusting for age, sex, height, and weight Journal of Bone and Mineral Research from the same visit as DXA and metabolic measurements. HRT users were excluded in the metabolomic study. Analyses were conducted separately in the baseline and follow-up studies and then meta-analyzed using inverse variance method.

Results
The demographic characteristics of the study populations are presented in Table 1.

Metabolomic signature of osteoporosis
Principal components analysis (PCA) showed 38 significant components in our panel with the cumulative R 2 of 55% in the training set and 35.5% in the cross-validation set (Q 2 ). O-PLS modeling confirmed one component as predictive and two as orthogonal for prediction of BMD in different regions. Models showed higher predictive power of metabolites for total hip BMD (R 2 ¼ 21.5% and Q 2 ¼ 17.7%) and femoral neck BMD (R 2 ¼ 18.9% and Q 2 ¼ 15.6%) compared with lumbar spine (R 2 ¼ 15.0% and Q 2 ¼ 10.5%). A scatter plot of observed versus predicted values for total hip BMD based on the O-PLS model is presented in Fig. 2.

Mendelian randomization
Results of the bidirectional MR analyses for all metabolites are shown in the Supplemental Material. In the direct instrumental variables analysis, four sulfated adrenal androgen levels (androsterone sulfate, epiandrosterone sulfate, 5alpha-androstan-3beta,17beta-diol disulfate, and 4-androsten-3beta,17beta-diol disulfate 1) showed causal associations with BMD traits (Table 3). We could not see any evidence for a reverse relationship (ie, BMD causing any metabolite level changes). The first three metabolites were instrumented on SNP rs4646450 (chromosome 7q22.1), an intron variant of the CYP3A5 gene (cytochrome P450, family 3, subfamily A, polypeptide 5). This gene encodes a member of the cytochrome P450 superfamily of enzymes. These enzymes catalyze many reactions involved in drug metabolism and synthesis of cholesterol, steroids (eg, testosterone, progesterone, and androstenedione), and other lipids. The other metabolite was instrumented on SNP rs16981893 (chromosome 19q13.33). This SNP is downstream of the SULT2A1 gene (sulfotransferase family, cytosolic, 2A, dehydroepiandrosterone-preferring, member 1). Sulfotransferases aid in the metabolism and excretion of drugs and endogenous compounds. The product of SULT2A1 catalyzes the sulfation of steroids and bile acids in the liver and adrenal glands.

Replication in the HKOS
To further evaluate the generalizability of these findings, we replicated these causally associated metabolites in the HKOS. All four metabolites were significantly replicated in the HKOS ( Table 4). The variance explained (adjusted R 2 ) by these metabolites were 2.7% and 1.2% for BMD at the total hip and lumbar spine, respectively, after adjustment for age, sex, height, and weight. The corresponding numbers for the TwinsUK population were 0.83% and 0.80%, respectively.

Discussion
To the best of our knowledge, this is the first study that uses two levels of "omics" data to detect the causal associations between blood-circulating metabolites and BMD variation. We confirmed the causal role of adrenal-secreted sulfated steroids in osteoporosis and identified two novel genetic loci (CYP3A5 and SULT2A1) to be important in the regulation of metabolites causing osteoporosis. Successful replication of the causally  associated metabolites in another ethnic group also demonstrated a high generalizability of our findings.
The contribution of metabolites into BMD variation is not yet established. Previous metabolomic studies have mainly focused on low BMD associations using low-resolution techniques instead of looking into BMD variations across the whole range in large population cohorts. (22,23) Using a high-throughput metabolomic platform measuring 280 known metabolites in a large population, we showed that circulating metabolites explained about 35% to 55% of BMD variation in all skeletal sites measured. This can have important clinical applications because circulating metabolite levels can be used as markers of bone health, predictors of risk of osteoporotic fractures, and markers of treatment effects or potential intervention targets (diet and medications).
Using the MR approach, we found that sulfated adrenal androgen plays a causal role in bone metabolism. These metabolites are related to androgen metabolism, and androgen has long been known as an anabolic agent in bone metabolism. However, whether the elevated levels of these metabolites represent the higher levels of endogenous androgen in the participants or contribute individually to bone metabolism is unknown and needs further assessment. Our previous study showed that one of these metabolites, epiandrosterone sulfate, was significantly associated with chronic widespread musculoskeletal pain. (24) The SNP rs4646450 has shown significant association with femoral neck BMD in the GEFOS consortium meta-analysis. (17) Importantly, all the causally associated metabolites in the TwinsUK study showed direct association with BMD in the HKOS population, and collectively explained 2.7% variance of BMD variation at the total hip. This raises hope for potential feasible metabolic interventions, either diagnostic or therapeutic, in clinical management of osteoporosis.
It is important to note that we cannot exclude a potential causal association of the metabolites associated with BMD in our study but not included in the MR analysis. This is mainly related to the lack of reliable genetic instruments for running MR analyses for these metabolites. For instance, prolyl-hydroxyproline showed significant association with BMD in all three measurement sites, but we could not find any genetic variant for it in our GWAS. This metabolite is a known marker of bone collagen degradation, and its urinary excretion has been previously linked to postmenopausal osteoporosis. (2,25) Caffeine also showed direct associations with BMD variation in our study, but this should be interpreted cautiously because the levels of this xenobiotic may only reflect short-term exposure to coffee in participants and not be related to the risk of osteoporosis associated with long-term coffee intake. The effect of caffeine on bone has been controversial throughout literature. (26) Our study shows that combining two sources of "omics" studies (genomics and metabolomics) would be beneficial in uncovering the pathogenesis of osteoporosis and other complex health conditions. Although genetic factors are wellknown markers of osteoporosis, combining 63 susceptibility variants in the largest genomewide meta-analysis only explained 5.8% of the variance of femoral neck BMD variation. (17) Although these genetic studies can bring several important clinical applications, (1) the genetic variations per se are not modifiable and so the findings cannot be directly translated into clinical practice. Using these genetic markers as instruments for  discovery of causally associated metabolites can reveal more insight on the pathological pathways to diseases and potential interventions.
Our study has several strengths. This is the first study that used two levels of "omics" data to detect the causal associations between blood-circulating metabolites and BMD variation. The results were further replicated in another ethnic group, highlighting the high reliability and generalizability of the findings. Both TwinsUK and HKOS are well-characterized cohorts for osteoporosis. The metabolomic profiling in TwinsUK and HKOS were performed on the same platform by the same provider with a highly stringent quality control, making the metabolomic data comparable between cohorts.
Our study has some limitations. There is a potential for horizontal pleiotropy in our study (ie, variants used as instruments can directly influence other metabolites, which then influence BMD). We could not employ new extensions of the MR method (namely multivariable MR (27) and MR Egger regression) (28) because of power limitations. Metabolites were measured in a relative manner and so their absolute concentrations cannot be inferred. Moreover, the difference in levels of metabolites between British and Chinese participants could not be assessed (given separate normalized values in each population). Patients receiving HRT were excluded from the metabolomics and genomics studies in the HKOS cohort at the design stage. This may cause some differences in the target populations of the cohorts, although use of HRT is not common in the Chinese population. False-negative causal associations are quite likely attributable to lack of validated genetic instruments for several metabolites. Future collaborative studies with larger sample sizes are required to evaluate the role of other BMD-associated metabolites in bone metabolism.
Our findings have important clinical implications. Currently, there is a crisis in the treatment of osteoporosis worldwide. (29) Patients are generally afraid of initiating or using antiosteoporosis treatments because of the fear of severe but rare side effects. Thus, the development of novel therapeutic agents for improving bone health is important. The causal metabolic pathway identified in this study may provide insight on therapeutic agent development. Moreover, the associated metabolites may be used as a biomarker of osteoporosis or potentially as a dietary supplement to improve BMD. Follow-up studies are now underway to evaluate the importance of these findings.
In conclusion, circulating metabolites measured using a high-throughput metabolomic technique explained more than one-third of the variance of BMD in a large female population. Fifteen unique metabolites showed direct association with BMD variation, four of which were found to be causally associated using genetic instruments. The direct association between these metabolites and BMD were replicated in an independent Chinese population. Combining both genomic and metabolomic data provided a unique platform for elucidating the role of new pathways in the pathogenesis and diagnosis of osteoporosis.

Disclosures
RPM is an employee of Metabolon, Inc. All other authors state that they have no conflicts of interest.