Multi-omics Data Integration for Identifying Osteoporosis Biomarkers and Their Biological Interaction and Causal Mechanisms.

Osteoporosis is characterized by low bone mineral density (BMD). The advancement of high-throughput technologies and integrative approaches provided an opportunity for deciphering the mechanisms underlying osteoporosis. Here, we generated genomic, transcriptomic, methylomic, and metabolomic datasets from 119 subjects with high (n = 61) and low (n = 58) BMDs. By adopting sparse multiple discriminative canonical correlation analysis, we identified an optimal multi-omics biomarker panel with 74 differentially expressed genes (DEGs), 75 differentially methylated CpG sites (DMCs), and 23 differential metabolic products (DMPs). By linking genetic data, we identified 199 targeted BMD-associated expression/methylation/metabolite quantitative trait loci (eQTLs/meQTLs/metaQTLs). The reconstructed networks/pathways showed extensive biomarker interactions, and a substantial proportion of these biomarkers were enriched in RANK/RANKL, MAPK/TGF-β, and WNT/β-catenin pathways and G-protein-coupled receptor, GTP-binding/GTPase, telomere/mitochondrial activities that are essential for bone metabolism. Five biomarkers (FADS2, ADRA2A, FMN1, RABL2A, SPRY1) revealed causal effects on BMD variation. Our study provided an innovative framework and insights into the pathogenesis of osteoporosis.

Introduction

Osteoporosis is a chronic progressive disorder characterized by low bone mineral density (BMD) and deterioration of bone microarchitecture, resulting in increased bone fragility and susceptibility to fracture (Kanis, 2002). The prevalence of osteoporosis in the aging population is over 20% in the United States, and it is becoming an increasingly serious public health problem in the elderly (Wright et al., 2014). Previous genetic studies have indicated that BMD is under strong genetic influence, with estimates of heritability ranging from 0.50 to 0.85 (Ralston and de Crombrugghe, 2006, Ralston and Uitterlinden, 2010). In the past decade, researchers have interrogated a wide variety of biological components (e.g., genetic variation, gene expression, and DNA methylation) and uncovered a number of risk biomarkers for BMD. Nonetheless, most of the prior studies have been focused on identification of biomarkers in a single molecular level through univariate statistical methods (e.g., t test, ANOVA, or linear model) and rarely integrated evidences from multiple omics levels. Consequently, so far, the specific functional roles of these identified molecular biomarkers are largely unknown, and their in vivo biological interaction and causal mechanisms are not explored.

The advancement of high-throughput technologies, such as whole genome sequencing (WGS), RNA-sequencing (RNA-seq), reduced-representation bisulfite sequencing (RRBS), and liquid chromatography-mass spectrometry (LC-MS), has dramatically increased our ability to comprehensively interrogate diverse molecular features at different omics levels. Meanwhile, several statistical integrative approaches have recently been developed for combining the molecular biomarkers identified from separate analyses of each omics (Liu et al., 2013, Gunther et al., 2012, Rohart et al., 2017), which lead to the discovery of crucial biological insights in a holistic manner and substantially enhance our understanding of molecular networks/pathways underlying the development of human complex diseases. In particular, Le Cao and her colleagues (Singh et al., 2019, Rohart et al., 2017) recently proposed a sparse multiple discriminative canonical correlation method that enables feature selection in multi-omics datasets and answers cutting-edge integrative questions in system biology. This method provided several attractive properties: (1) it allows relax assumptions about data distribution and thus is highly flexible to answer topical questions across various biology-related fields; (2) it is computationally efficient to handle large datasets, where the number of biological markers is much larger than the number of samples; (3) it implements dimension reduction by projecting the data into a smaller subspace while capturing the correlation structure and highlighting the largest sources of variation from the data, resulting in a powerful explanation of the biological system under study.

Genetic variation has a substantial impact on multiple genomic contexts and/or molecular/cellular phenotypes in humans (Pierce et al., 2018, Albert and Kruglyak, 2015, Lemire et al., 2015, McVicker et al., 2013, Kraus et al., 2015), such as transcript abundance (Pierce et al., 2018, Albert and Kruglyak, 2015), DNA methylation (Pierce et al., 2018, Lemire et al., 2015), histone modification (McVicker et al., 2013), and metabolites (Kraus et al., 2015). Whole-genome association scans to detect regions that harbor such variants for gene expression (known as expression quantitative trait locus, eQTL), DNA methylation (meQTL), and metabolite product (metaQTL) have been conducted in multiple human cell/tissue types (GTEx Consortium et al., 2017). Interestingly, previous studies have shown that many QTLs may appear to influence multiple molecular phenotypes. For instance, single nucleotide polymorphisms (SNPs) associated with expression of nearby genes were often associated with methylation of nearby CpG sites (Pierce et al., 2018). These common QTLs in multiple phenotypes may suggest a potential shared biological mechanism by which the common causal variant influences both gene expression and DNA methylation en route to eventually influencing phenotypes. Nonetheless, the precise mechanisms underlying these genetic associations remain poorly understood because of short of the study approaches. Thus, analytical approaches for dissecting the complex biological processes are needed to prioritize the plausible functional variants for further studies. Notably, Mendelian randomization (MR) analysis has been widely used to assess potential causal relationships of genetic/environmental risk factors and diseases (Davey Smith and Hemani, 2014). Recently, MR analysis has been adopted to inspect the causality of biomarkers in disease etiology, utilizing multiple independent SNPs identified by QTL analysis (QTL SNPs) as instrumental variables (Taylor et al., 2019, Yao et al., 2018, Chen et al., 2018). As an example, by applying eQTLs as genetic instruments, Chen et al. recently revealed a causal relationship between LINC00339 gene expression and BMD variation (Chen et al., 2018).

In this work, we performed multi-omics integrative analyses with the largest datasets so far in the bone field to identify osteoporosis biomarkers as well as their biological interaction and causal mechanisms. A simple overview of our workflow is illustrated in Figure 1. Briefly, our approach consisted of four phases. First, we performed individual transcriptomic, methylomic, and metabolomic analyses in 119 Caucasian female subjects with high (n = 61) and low (n = 58) BMDs to identify potential differentially expressed genes (DEGs), differentially methylated CpG sites (DMCs), and differential metabolic products (DMPs) for osteoporosis risk. The basic characteristics of the study subjects were summarized in Table S1. Second, we integrated the identified DEGs, DMCs, and DMPs via a sparse multiple discriminative canonical correlation analysis (SMDCCA) to retrieve prominent osteoporosis biomarkers that not only reliably distinguish the high-BMD and low-BMD groups, but also are highly correlated spanning different biological layers. Third, we used targeted QTL analyses to test the effects of SNPs on prominent osteoporosis biomarkers in each omics, followed by interaction network analyses, as well as functional annotation and classification analyses, to assess the biological importance of the identified biomarkers. At last, by applying MR analyses using the multiple independent QTL SNPs as instrumental variables, we assessed the causality of the functionally classified biomarkers in BMD variation and inspected whether the identified biomarkers are causally related to one another (e.g., by functional regulation) with the purpose of gaining insights into the in vivo molecular functional mechanisms of the etiology of osteoporosis. In aggregate, we identified several osteoporosis biomarkers and reconstructed multi-omics networks/pathways that may mediate variation in risk of osteoporosis in vivo in humans. Our study pioneered an innovative integrative approach, and the results illuminated the advantages of multi-omics integrative analysis and provided valuable insights into the pathogenic mechanisms of osteoporosis.

Figure 1

The Overall Workflow for Identifying Osteoporosis Biomarkers and Their Biological Interaction and Causal Mechanisms

The workflow consisted of four phases. Phase 1: individual transcriptomic, methylomic, and metabolomic analyses. Phase 2: SMDCCA integrative analysis of DEGs, DMCs, and DMPs. Phase 3: targeted QTL analyses followed by interaction network analyses, as well as functional annotation and classification analyses. Phase 4: MR analyses. PBMs, peripheral blood monocytes; DEGs, differentially expressed genes; DMCs, differentially methylated CpG sites; DMPs, differential metabolic products; SMDCCA, sparse multiple discriminative canonical correlation analysis; IV, instrumental variable; MR, Mendelian randomization.

Results

Multi-omics Integration with Supervised SMDCCA

A total of 25,342 genes, 17,462,566 CpG sites, and 4,209 metabolites were measured in 61 subjects with high BMD and 58 subjects with low BMD, of which 18,774 genes, 763,265 CpG sites, and 2,608 known metabolites were selected (Methods) for subsequent analyses. We identified a total of 1,594 DEGs with false discovery rate (FDR) < 0.05, 1,219 DMCs (q < 0.05 with methylation difference large than 10%), and 204 DMPs (p < 0.05) by comparing the high-BMD and low-BMD groups in prevailing single-omics analyses.

By multi-omics integrative analysis on DEGs, DMCs, and DMPs via SMDCCA approach, we identified an optimal multi-omics biomarker panel for discriminating high-BMD and low-BMD groups (classification error rate of 0.1, Figure S1). This optimal multi-omics biomarker panel was composed of three components (Table S2) involving a total of 74 DEGs, 75 DMCs, and 23 DMPs (henceforward termed prominent osteoporosis biomarkers). We observed moderate correlations between DEGs and DMCs/DMPs and a few weak correlations between DMCs and DMPs (Figures 2 and S2). Notably, a substantial proportion of DEGs were found to be correlated with both DMCs and DMPs. For instance, gene expression of HAUS2 has positive correlation with DNA methylation at CpG sites Chr1:5874307 (NPHP4) but is negatively correlated with metabolite threonine. Interestingly, there is a significant interaction between HAUS2 and NPHP4 with a combined interaction score of 0.905 (Szklarczyk et al., 2019). However, the specific functional roles of threonine and mechanisms underlining the connection between these correlations are unclear.

Figure 2

Circos Plot for Prominent Osteoporosis Biomarkers Selected in the Optimal Multi-omics Biomarker Panel

Circos plot displays the different types of osteoporosis biomarkers in the first component on a circle, with links between each omics indicating the positive (brown) or negative (black) correlations with cutoff r = 0.4. The blue line and orange line represent the high-BMD and low-BMD groups, respectively.

Next, we examined the regulatory status for regions containing the 75 DMC prominent osteoporosis biomarkers. Our results showed that 51 (68.0%) DMCs were mapped to the DNaseI hypersensitivity cluster; 28 (37.3%) DMCs were assigned to the H3K27ac chromatin mark, which was often observed near active regulatory elements; and 39 (52.0%) DMCs also showed gene expression signals in multiple tissues from GTEx RNA-seq data (Figure 3A, Table S3). In addition, there were 28 DMCs at sites that were predicted to be fast-evolving (PhyloP score < −1) and 6 DMCs at sites that were predicted to be conserved (PhyloP score >1) (Table S3). Interestingly, a previous study has indicated that the methylation-PhyloP correlation is dependent on the sequence context (Chuang and Chen, 2014), although the exact mechanisms are still unclear. Therefore, the functional properties of these DMC biomarkers are worth further exploration.

Figure 3

Functional Enrichment of Osteoporosis Biomarkers

Venn diagrams showing the functional enrichment of DMCs (A), BMD-associated QTL SNPs (B), and common QTLs (C) in different potential regulatory elements. The detailed results can be found in Tables S3, S5, and S8.

Targeted eQTL, meQTL, and metaQTL Analyses

To assess the genetic effects on the prominent osteoporosis biomarkers (74 DEGs, 75 DMCs, and 23 DMPs), we performed targeted eQTL, meQTL, and metaQTL analyses by utilizing WGS data generated from the same subjects. At significance threshold of FDR = 0.05, we detected 6,778 SNP-DEG pairs, corresponding to 64 DEGs and 4,401 eQTL SNPs (eSNPs). Significant associations were also detected at 2,062 meQTL SNPs (meSNPs) for 42 DMCs and 1,380 metaQTL SNPs (metaSNPs) for 23 DMPs (Table 1). Notably, a total of 199 QTL SNPs, including 67 eSNPs, 126 meSNPs, and 6 metaSNPs, were associated with BMD (p < 5.0×10−8) in previous genome-wide association studies (GWASs) (Kemp et al., 2017, Morris et al., 2019) (Figure 3B, Table S4), supporting their impact on the BMD phenotype. Among the 67 BMD-associated eSNPs, 59 (88.1%) SNPs were predicted to alter transcription factor (TF) binding motifs and 50 (74.6%) SNPs were annotated to enhance histone marks (Figure 3B, Table S5). Similar enrichment for potential regulatory elements were also observed for the 126 meSNPs and 6 metaSNPs associated with BMD (Figure 3B, Table S5). Remarkably, we detected 28 highly correlated (linkage disequilibrium, LD r2 > 0.8) BMD-associated eSNPs for biomarker FADS2. These eSNPs lie in the gene cluster of MYRF, FEN1, FADS1, and FADS2 at the 11q12.2 locus, of which 3 eSNPs were mapped to transcription start sites (TSSs) of FEN1 (rs174538) and FADS2 (rs5792235, rs99780), 1 eSNP rs1535 mapped to FADS2 enhancer, and 1 eSNP rs174562 overlapped with both TSS and enhancer of FADS1/FADS2 in a wide variety of cell/tissue types (Figure 4). Interestingly, these 5 eSNPs also showed the top five highest functionality scores among the 28 eSNPs (Figure S3, Table S6) through 3DSNP prioritization analysis (Lu et al., 2017). This result is supported by a previous functional study, which demonstrated the influence of rs174538 on the expression of FEN1 and its enzyme activity (Yang et al., 2009).

Table 1

Patterns of eQTLs, meQTLs, and metaQTLs

Type of Analysis	Significant SNP-Biomarker Pairs	QTL SNPs	Biomarkers
eQTL	6,778	4,401 (327)	64
meQTL	2,387	2,062 (75)	42
metaQTL	1,559	1,380 (96)	23

Note: The association between SNP and biomarker was tested with a linear regression model in R package Matrix eQTL. The significance threshold was defined as FDR <0.05. QTL SNPs: The SNPs identified by QTL analysis. Biomarkers: DEGs, DMCs, and DMPs. The number in the bracket indicated the independent SNPs after LD pruning (r2 < 0.1).

Figure 4

Chromatin-State Annotation of 28 BMD-Associated eSNPs for Biomarker FADS2

The chromatin state annotation tracks were generated by the 18-state ChromHMM model from the Roadmap Epigenomics Project under the human reference genome assembly GRCh37 (hg19) and visualized in the UCSC Genome Browser. eSNP rs174538, rs5792235, rs99780, and rs174562 were mapped to transcription start site (TSS); eSNP rs1535 and rs174562 were mapped to enhancer in a wide variety of cell/tissue types.

Next, we attempted to identify common QTLs across DEGs, DMCs, and DMPs. A total of 448 potential eSNPs were also associated (p < 5.0×10−5) with methylation at 48 DMCs, which corresponded to 883 SNP-DEG-DMC combinations (Table S7). For example, SNP rs2236373 is associated with gene expression of BMP3 (p = 3.45×10−6) and DNA methylation at CpG site Chr16:75279661 (p = 4.42×10−5) in BCAR1 gene. Interestingly, both biomarkers were reported to be associated with bone-related signaling pathways, such as WNT pathway (Kokabu and Rosen, 2018) and RANK/RANKL pathway (Robinson et al., 2009). Similarly, we identified 451 SNP-DEG-DMP combinations and 148 SNP-DMC-DMP combinations (Table S7). Functional annotation analyses showed that many of these multi-marker QTLs were predicted to alter TF binding motifs and/or mapped to other putative regulatory regions (Figure 3C, Table S8), suggesting that a substantial number of genetic variants may have multi-level effects on functional biomarkers across different omics levels, possibly through regulation via complex functional network/modules across multi-omics layers.

Functional Interaction Network of the Prominent Osteoporosis Biomarkers

Through STRING interaction analysis (Szklarczyk et al., 2019), we revealed functional interactions among the 74 DEGs and/or 75 DMCs (corresponding to 61 DMC-annotated genes) that were identified in the optimal multi-omics biomarker panel (Figures 5 and S4). Notably, by integrating 74 DEGs and 61 DMC-annotated genes together, we revealed a complex interaction network (Figure 5), which may indicate the potential regulatory relationships between gene expression and DNA methylation biomarkers in bone metabolism.

Figure 5

A Functional Interaction Network for DEGs and DMC-Annotated Genes

DMCs were assigned to their nearest gene/gene cluster based on the human reference genome assembly GRCh37 (hg19). Connections are based on co-expression and experimental evidence with a STRING v11.0 summary score above 0.4. Each color-filled node denotes a querying gene; edges between nodes indicate protein-protein interactions between protein products of the corresponding genes. Different edge colors represent the types of evidence for the association (Szklarczyk et al., 2019). Note: * represents DEGs.

Furthermore, by integrating the STRING interaction results with the known biological functions (through literature review) of each prominent osteoporosis biomarker (74 DEGs, 61 DMC-annotated genes, and 23 DMPs), we assigned these biomarkers into different signaling pathways or functional activities related to bone metabolism. Finally, we identified 29 DEGs, 38 DMCs-annotated genes, and 8 DMPs that may interactively act upon several critical bone-related signaling pathways/activities in osteoblasts and/or osteoclasts (Figure 6), such as RANK/RANKL pathway, MAPK/TGF-β pathway, WNT/β-catenin pathway, G protein-couple receptor activity, GTP binding/GTPase activator activity, and telomere/mitochondrial activity.

Figure 6

Biomarkers in Bone-Related Signaling Pathways/Activities

The 29 gene expression (DEGs), 38 DNA methylation (DMC-annotated genes), and 8 metabolic (DMPs) biomarkers that may interactively act upon several critical bone-related signaling pathways/activities in osteoblasts and/or osteoclasts. The purple line indicates interaction identified in the STRING interaction network, and the gray line indicates the association reported in the literature. The known osteoporosis biomarkers reported in previous GWASs or functional studies were marked in bold italics. * Several polypeptide members of TGF-β/TGF-β receptor and their coactivators (TGFB1I1, TGFBR3, and TGFBRAP1) were differentially expressed in single-omics analysis. DOPEY2 (CpG_Chr21:37635103) was also identified as a DNA methylation biomarker.

MR Analysis

The finding of functionally classified biomarkers prompted us to investigate their causal biological mechanisms in BMD variation. Focusing on the prominent osteoporosis biomarkers assigned to signaling pathways/activities critical for bone metabolism (Figure 6), we conducted MR analysis and identified five biomarkers (DEG biomarkers FADS2, ADRA2A, FMN1, RABL2A, and DMC biomarker CpG_4:124356866 at SPRY1) that may have putative causal effects on BMD variation (Table 2). Interestingly, gene expression of FADS2 has a robust causal effect on BMD variation based on either median-based method or inverse-variance weighted (IVW) method (p < 0.001). Previous GWAS meta-analysis has shown that the genetic variants in FADS2 were associated with BMD (Kemp et al., 2017). Furthermore, our QTL analysis revealed that these variants have significant effects on FADS2 expression (FDR <0.05, Table S4). Collectively, these results provided convergent and compelling evidence for the significance of genetic regulation of FADS2 expression in bone metabolism.

Table 2

Significant Causal Biomarkers for BMD Variation

Causal Biomarker	Simple Median	Weighted Median	IVW	MR-Egger	Intercept
FADS2	<0.001	<0.001	<0.001	0.23	0.92
ADRA2A	0.014	0.014	0.044	0.216	0.057
FMN1	0.032	0.023	0.005	0.154	0.626
RABL2A	0.123	0.105	0.185	0.041	0.078
CpG_4:124356866 (SPRY1)	0.043	0.078	0.026	0.21	0.026

Note: significant results (p < 0.05) are marked in bold.

Furthermore, we assessed whether DEGs, DMCs, and DMPs within the same signaling pathway/activities are causally related to one another. As an example, we selected two biomarkers from RANK/RANKL pathway (Figure 6), namely, DEG biomarker ADCY3 and DMC biomarker NFATC1 (CpG_18:77225621). ADCY3, inducted by RANKL, encodes a membrane-associated enzyme that catalyzes the formation of secondary messenger cAMP in response to G protein-coupled receptor activity. NFATC1, a ubiquitous TF in many cell types, is a well-known master regulator of both osteoblastogenesis and osteoclastogenesis (Fromigue et al., 2010, Kim and Kim, 2014) and has been demonstrated to play crucial roles in regulating bone homeostasis and bone mass (Winslow et al., 2006, Lee et al., 2009). Previous studies have shown that the expression of NFATC1 was markedly elevated in ADCY3-silenced cells since the elevation of intracellular cAMP culminates the PKA-mediated phosphorylation and subsequently inhibits gene expression of NFATC1 (Yoon et al., 2011, Sheridan et al., 2002). Nevertheless, no interaction was reported in STRING interaction network. Notably, in our study, MR analyses revealed a significant causal effect of gene expression of ADCY3 on DNA methylation of CpG_18:77225621 at NFATC1 based on either median-based method or IVW method with p < 0.001 (Table 3). Interestingly, NFATC1 was identified not only through a DMC (CpG_Chr18:77225621) (Figure 7A), but also as a DEG (p = 4.37×10−2), although it did not quite (but close to) reach the adjusted threshold for integrative analysis (Figure 7B). Remarkably, functional annotation analysis revealed that DNA methylation biomarker CpG_Chr18:77225621 was linked to regulatory regions, such as TSSs, transcription regions, enhancer histone marks, and DNaseI hypersensitivity cluster (Figure 7C, Table S3). GTEx RNA-seq data also showed extensive gene expression signals of NFATC1 across multiple human cell/tissue types (Figure 7D). Taken together, these results suggested that DNA methylation changes in NFATC1 may regulate its gene expression activity and ultimately regulates bone metabolism.

Table 3

Significant Causal Effects of Gene Expression of ADCY3 and/or ADRA2A on NFATC1 (CpG_18:77225621) DNA Methylation

Causal Biomarker	Simple Median	Weighted Median	IVW	MR-Egger	Intercept
ADCY3	<0.001	<0.001	<0.001	0.525	0.451
ADRA2A	<0.001	<0.001	0.002	0.009	0.068
ADCY3 + ADRA2A	<0.001	<0.001	<0.001	0.002	0.097

Note: significant results (p < 0.05) are marked in bold.

Figure 7

Biological Importance of Biomarker NFATC1

(A and B) Boxplot for NFATC1 DNA methylation (A) and gene expression (B) levels in high-BMD and low-BMD groups. The vertical axis in (A) represents DNA methylation level (M-value); the adjusted q-value was determined by R package methylKit. The vertical axis in (B) represents the gene expression level [log2(TMM)]; p values for moderated statistics were determined by R package limma.

(C) The chromatin-state annotation for DNA methylation biomarker CpG_Chr18:77225621 at NFATC1. The chromatin state annotation tracks were generated by the 18-state ChromHMM model from the Roadmap Epigenomics Project under the human reference genome assembly GRCh37 (hg19) and visualized in the UCSC Genome Browser.

(D) The gene expression signals of NFATC1 across multiple human cell/tissue types from GTEx RNA-seq data.

In addition, DEG biomarker ADRA2A, a key player in G protein-coupled receptor activity, was highly associated with ADCY3 (evidence score 0.919) in STRING interaction network. Interestingly, our MR analysis also showed a significant causal impact of ADRA2A gene expression on BMD variation (Table 2). Therefore, we were interested in testing whether ADRA2A by itself or together with ADCY3 also has a causal effect on NFATC1 DNA methylation. Indeed, we identified a significant causal effect of ADRA2A gene expression or combined effects of ADCY3 and ADRA2A expression on DNA methylation of NFATC1 (Table 3). In contrast, we also conducted reversed causality analysis (e.g., test the causal effect of DNA methylation of NFATC1 on ADCY3 and/or ADRA2A gene expression) and found no significant results (data not shown).

Discussion

Recent development in high-throughput profiling technologies and integrative analysis of multi-omics data offered advanced and powerful approaches to dissect complex biological problems. In this study, we pioneered an innovative approach by synthesizing the state-of-the-art methods recently developed and performed multi-omics analyses integrating gene expression and DNA methylation data in bone-related cells, as well as serum metabolomics data. We initially identified 74 DEGs, 75 DMCs (in 61 genes), and 23 DMPs for BMD variation (Table S2). There are a set of osteoporosis biomarkers in addition to biomarkers that are known to be associated with BMD. We then investigated the effects of genetic variants on these prominent osteoporosis biomarkers via targeted QTL analysis in each omics and identified hundreds of potential QTL SNPs shared by different omics, which may suggest their common biological mechanisms in pathogenesis of osteoporosis. Furthermore, we reconstructed the integrative networks/pathways in bone metabolism via STRING interaction analysis combined with functional annotation and classification analyses and revealed that a substantial proportion of these biomarkers not only interacted with each other, but also were enriched in several well-known signaling pathways or functional activities (Figure 6) that are crucial for osteoblastogenesis and osteoclastogenesis, such as RANK/RANKL pathway, MAPK/TGF-β pathway, and WNT/β-catenin pathway, as well as G protein-coupled receptor activity, GTP binding/GTPase activator activity, and telomere/mitochondrial activity. By considering the perplexing relationship between functionally classified biomarkers, we implemented MR analysis to investigate their potential causality. Our MR results further provided supporting evidence that several gene expression and DNA methylation biomarkers have causal effects on the final BMD variation or were causally related to one another. In aggregate, our multi-omics data integration captured the complexity of the prominent interplay among different omics and pointed out a list of candidate biomarkers that may help refine biological hypotheses and propose biological validations for future studies. Furthermore, the integration framework taken here can be adopted to other complex traits/disorders and further extended to incorporate additional types of omics data (e.g., proteomics, lipidomics, and metagenomics) to enhance our understanding of the pathogenesis of human diseases.

In addition to the osteoporosis biomarkers that have been discussed in the results, several other prominent osteoporosis biomarkers that participated in well-known signaling pathways or functional activities of bone metabolism (Figure 6) should also be concerned and highlighted. Briefly, there are 18 genes, including 5 DEGs (ADCY3, DNASE1, HAUS2, GATA1, and FMN1) and 13 DMC-annotated genes (MN1, SPRY1, KCNQ1, NFATC1, ITPKB, BCAR1, NPHP4, MFHAS1, GAS6, PDE9A, RTEL1-TNFRSF6B, MFN2-TNFRSF8-TNFRSF1B, and ESPN), annotated in RANK/RANKL pathway. MN1 acts as a transcriptional activator of the osteoclastogenic cytokine RANKL and plays a crucial role in the formation of the membranous bones in the skull during mammalian development (Zhang et al., 2009). Disruption of MN1 in calvarial osteoblasts resulted in altered morphology, decreased growth rate, impaired motility, and attenuated 1,25(OH)2D3/VDR-mediated transcription, as well as reduced alkaline phosphatase activity and mineralized nodule formation (Zhang et al., 2009). SPRY1 encodes a growth factor regulator for marrow progenitor cells and promotes osteoblast differentiation at the expense of adipocytes (Urs et al., 2012). A recent transgenic mouse model revealed that miR-21, a regulator of osteoclastogenesis, can affect RANK/RANKL signaling pathway by targeting SPRY1 (Hu et al., 2017). Genetic variants in SPRY1 have also been associated with osteoporosis in Korean women (Jin et al., 2013).

We identified 2 DEGs (PIM1, BMP3), 6 DMC-annotated genes (MAPK11, PMEPA1, GSDMD, ADAMTSL2, ADAMTS17, and POFUT2-COL18A1), and 1 metabolite biomarker LysoPC (16:0) in MAPK/TGF-β pathway. PIM1 is a member of the serine/threonine kinase family and can significantly decreases MAP3K5 kinase activity and inhibits MAP3K5-mediated phosphorylation of JNK and JNK/p38MAPK, which subsequently reduces caspase-3 activation and cell apoptosis (Gu et al., 2009). Importantly, PIM1 can also regulate RANKL-induced osteoclastogenesis via NF-κB activation and NFATC1 induction (Kim et al., 2010). BMP3 encodes a secreted ligand of the TGF-β superfamily of proteins. It is one of the most abundant bone morphogenetic proteins in demineralized bone matrix. BMP3 suppresses osteoblastogenesis and negatively regulates bone density by modulating TGF-β receptor availability to other ligands (Wu et al., 2016). Remarkably, we indeed observed the negative regulation of BMP3 gene expression on BMD levels (p = 2.59×10−5, Figure S5A). Several polypeptide members of TGF-β/TGF-β receptor and their coactivators (TGFB1I1, p = 2.80×10−3; TGFBR3, p = 1.26×10−5; TGFBRAP1, p = 7.77×10−4) were also differentially expressed between the high-BMD and low-BMD groups (Figures S5B–S5D). MAPK11, also known as p38-β, is one of the four p38 MAPKs that play a crucial role in osteoblast differentiation and bone development and maintenance (Hu et al., 2003). A recent study showed that MAPK11 can also enhance osteoclastogenesis and bone resorption (He et al., 2014). PMEPA1 encodes a transmembrane protein that contains a Smad interacting motif (SIM). PMEPA1 has a significant role in osteoclastogenesis (Funakubo et al., 2018) and can also act as a TGF-β signaling regulator in osteoblast proliferation (Fournier et al., 2015). ADAMTSL2 is directly involved in TGF-β bioavailability and plays a key role in osteoblast and skeletal development (Le Goff et al., 2008).

In addition, there are 2 DEGs (RBP1, LIMD1) and 4 DMC-annotated genes (GNAS, ILKAP-PER2, SMYD3, and UBE3C-DNAJB6) in WNT/β-catenin pathway. RBP1 can act as a RUNX2 coactivator and promotes osteoblastic differentiation (Monroe et al., 2010). LIMD1 encodes a scaffold protein that has been implicated in the regulation of osteoclastogenesis through an interaction with the p62/sequestosome protein (Luderer et al., 2008). LIMD1 protein can also influence osteoblast differentiation and function; as such, Limd1(−/−) calvarial osteoblasts displayed increased mineralization and accelerated differentiation (Luderer et al., 2008). Furthermore, there is a significant increase in nuclear beta-catenin staining in differentiating Limd1(−/−) calvarial osteoblasts (Luderer et al., 2008), suggesting that LIMD1 is a negative regulator of canonical WNT signaling in osteoblasts. The GNAS gene is a complex imprinted locus that produces multiple transcripts (such as Gsα, XLAS, NESP55) through the use of alternative promoters and alternative splicing. A recent study by Ramaswamy et al. (Ramaswamy et al., 2017) demonstrated that Gnas inactivation in mice negatively affects cortical bone quality and strength, with mutation of the paternal allele causing more severe effects than maternal mutations. These effects of Gsa deletion on bone maintenance were exerted through enhanced osteoclast differentiation and increased bone resorption, mediated by Gsa signaling via cAMP/PKA and WNT/β-catenin pathways (Ramaswamy et al., 2017). SMYD3 encodes a histone methyltransferase that functions in RNA polymerase II complexes by an interaction with a specific RNA helicase (Hamamoto et al., 2004) and controls a WNT-responsive epigenetic switch (Wang et al., 2018).

G protein-coupled receptor activity includes 2 DEGs (ADRA2A, EFR3B) and 4 DMC-annotated genes (HRH4, ACKR3, GPR78, and GPR124). GTP binding/GTPase activator activity includes 7 DEGs (METTL7A, STX1A, ARHGAP26, RABL2A, RABL3, SRL, and DBF4) and 6 DMC-annotated genes (DOCK2, BAHCC1, MICALL2-INTS1, RAB35, A4GALT-ARFGAP3, and TUBB6). For example, ADRA2A, a member of the G protein-coupled receptor superfamily, is involved in neuro-endocrine regulation of bone resorption (Mlakar et al., 2015). ARHGAP26 encodes a GTPase-activating protein, and a mutation in this gene has recently been determined to be associated with BMD (Kemp et al., 2017). Another interesting gene is DBF4, which plays a central role in DNA replication and cell proliferation through nitrogen-containing bisphosphonate-induced cytotoxicity (Bivi et al., 2009). Nitrogen-containing bisphosphonate can potently inhibit the prenylation and function of GTP-binding proteins required for osteoclast formation and now is firmly established as first-line therapy for osteoporosis (Grey and Reid, 2006). TUBB6 plays a key role in GTP binding and has been associated with BMD variation (Daswani et al., 2015). In addition, there are 8 DEGs (TERF1, MRPL10, NSUN4, C11orf83, PDSS2, PACS2, NDUFV3, and ALDH3B1), 6 DMC-annotated genes (PPARGC1A, ILKAP-PER2, SIRT6, RTEL-TNFRSF6B, VARS2-GTF2H4, and MFN2-TNFRSF8-TNFRSF1B), and 3 metabolites biomarkers ({4-((E)-2-(2,3,5-trihydroxyphenyl)ethenyl)phenyl}oxidanesulfonic acid, Acetyl-T2 toxin, and LysoPC (16:0)) associated with telomere/mitochondrial activity. TERF1 encodes a component of the telomere nucleoprotein complex, which can regulate telomere elongation and plays a key role in aging-related disease (Blasco, 2005). Previous study has shown that defects in telomere maintenance molecules impair osteoblast differentiation and promote osteoporosis (Pignolo et al., 2008). PPARGC1A encodes a transcriptional coactivator that mediates mitochondrial biogenesis and energy metabolism (Liang and Ward, 2006). This coactivator interacts with PPARG, which permits the interaction of this coactivator with multiple transcription factors (Vega et al., 2000). PPARGC1A can control skeletal stem cell fate and bone-fat balance in osteoporosis and skeletal aging (Yu et al., 2018). Interestingly, WNT signaling can activate PPARGC1A expression and upregulate mitochondrial biogenesis; this upregulation contributes to the osteoblastic differentiation (An et al., 2010). PER2 encodes a transcriptional repressor that forms a core component of the circadian clock. It directly and specifically represses PPARG proadipogenic activity by blocking PPARG recruitment to target promoters and thereby inhibiting transcriptional activation. PER2 is required for fatty acid and lipid metabolism and is involved as well in the regulation of circulating insulin levels (Grimaldi et al., 2010). SIRT6 encodes a member of the sirtuin family of NAD-dependent enzymes that are implicated in cellular stress resistance, genomic stability, aging, and energy homeostasis. A recent study demonstrated that the way SIRT6 regulated osteoclast was predominantly through osteoblast paracrine manner, rather than osteoclast-autonomous behavior, which provided a valuable insight into the pathogenesis of osteoporosis due to SIRT6 mutation (Zhang et al., 2018). Metabolite Acetyl-T2 toxin is the class of organic compounds known as trichothecenes. Trichothecene inhibition of protein synthesis in the mitochondria allows reactive oxygen species (ROS) to build up in the cell, which inevitably leads to oxidative stress and induction of the programmed cell death pathway, apoptosis (Fang et al., 2012). LysoPC (16:0) is a lysophosphatidylcholine. Previous study showed that LysoPC-induced p38 MAPK signaling pathway can control monocyte migration (Tan et al., 2009) and a novel lysophosphatidylcholine derivative (SCOH) can inhibit osteoclast differentiation and bone resorption (Kwak et al., 2004). Moreover, lysophosphatidylcholine can produce mitochondrial ROS generation, increase intracellular free calcium concentration, activate active adenylate cyclase (e.g., gene expression biomarker ADCY3), and enhance glucose-dependent insulin secretion via an orphan G protein-coupled receptor (Watanabe et al., 2006, Chaudhuri et al., 2003, Soga et al., 2005). LysoPC(16:0) isolated from rats plasma was also proved to be related to osteoporosis (Liu et al., 2012).

Notably, there are two adjacent gene clusters RTEL1-TNFRSF6B and MFN2-TNFRSF8-TNFRSF1B involved in multiple signaling pathways or functional activities. RTEL1 encodes a DNA helicase that interacts with proteins in the shelterin complex and plays a key role in the stability, protection, and elongation of telomeres (Deng et al., 2013). Interestingly, telomere deficiency can impair osteoblast differentiation and promote osteoporosis (Pignolo et al., 2008). TNFRSF6B can suppress RANKL-induced osteoclastogenesis via down-regulating NFATC1 and enhancing cell apoptosis (Cheng et al., 2013). MFN2 encodes a mitochondrial membrane protein that participates in the maintenance and operation of the mitochondrial network (Bach et al., 2003), which has been linked to osteoclast activity and bone metabolism via an iron-related fundamental pathway (Ishii et al., 2009). A recent study reported that MFN2 can facilitate osteoclastogenesis by regulating the calcium-calcineurin-NFATC1 axis as well (Szklarczyk et al., 2019). TNFRSF8 and TNFRSF1B both are members of the TNF-receptor superfamily. Genetic variants in TNFRSF1B gene have been associated with femoral neck BMD (Albagha et al., 2002) and bone structure (Mullin et al., 2008).

We also identified eight unclassified osteoporosis biomarkers, including 2 DEGs (FADS2, CPM), 1 DMC-annotated gene GREM2, and 5 metabolite biomarkers (pipecolic acid, threonine, methylmalonic acid, N-lactoyl-tryptophan, and nicotinic acid), which have been reported in previous association studies. Genetic variants in FADS2, CPM, and GREM2 have been associated with BMD variation in previous GWAS (Kemp et al., 2017). Pipecolic acid is a normal human metabolite present in human blood; it has been associated with both total hip and lumbar spine BMD phenotypes in TwinsUK population (Moayyeri et al., 2018). Interestingly, we also identified two function-related metabolites (N-lactoyl-tryptophan, nicotinic acid). N-lactoyl-tryptophan is lactoyl derivative of tryptophan. Tryptophan acts as the precursor of nicotinic acid (also known as vitamin-B3). Previous study has indicated that tryptophan plays an essential role in osteoblastic differentiation (Pallu et al., 2012). Recently, Michalowska et al. investigated the influence of tryptophan and its metabolites on bone remodeling and observed significant changes of tryptophan levels in bone metabolic diseases (Michalowska et al., 2015). Nicotinic acid occurs naturally in food. Several studies have examined the effect of nicotinic acid on bone metabolism. For example, there is a positive correlation between dietary intake of nicotinic acid and BMD in premenopausal Japanese women (Sasaki and Yanagibori, 2001). The other study reported a significant inverse association of dietary nicotinic acid intake with hip BMD, but there was no significant association with total body BMD (Carbone et al., 2019). These conflicting findings deserve further serious investigation to better understand the effect of supplementation of nicotinic acid on bone biology. Methylmalonic acid is a dicarboxylic acid that is a c-methylated derivative of malonate. Methylmalonic acid was found to induce osteoclastogenesis in a dose-dependent manner, and vitamin-B12 deficiency may lead to decreased bone mass by increased osteoclast formation due to increased methylmalonic acid level (Vaes et al., 2009). Threonine is an essential amino acid that is used in the biosynthesis of proteins. A prior study has demonstrated that threonine can modulate the growth and the differentiation of osteoblasts cultured in vitro and confirmed the relationship between osteoporotic hip fracture and inadequate protein intake (Conconi et al., 2001). Notably, a recent metabolomic study has shown that threonine is associated with BMD and can improve the power for osteoporosis classification in males (Wang et al., 2019).

In summary, we conducted an innovative multi-omics integrative analysis and identified a set of osteoporosis biomarkers as well as biological pathways/networks that may contribute to BMD variation. Our results revealed valuable insights into the pathogenesis of osteoporosis and aided in generating hypotheses for future functional studies.

Limitations of the Study

Several limitations of this study should be noted. First, our sample size is relatively small. However, our study is the largest so far in the bone field for multi-omics analyses and we applied an extreme phenotype sampling strategy with stringent inclusion and exclusion criteria, which is known to provide enhanced statistical power for association analysis compared with studies using comparable numbers of randomly sampled subjects (Bjornland et al., 2018). Moreover, data from different omics levels can provide complementary and inherent validation information with each other, and thus, integrating multi-omics data can partially compensate for the relatively small sample sizes (Hasin et al., 2017). But clearly, our results need to be validated in future studies with large sample size. Second, we used a relatively homogeneous cell type, peripheral blood monocytes (PBMs), as model cells for gene expression and DNA methylation analysis for osteoporosis. PBMs can act as osteoclast precursors, secrete cytokines essential for osteoclast differentiation and function, and represent a major target cell of sex hormones for bone metabolism (Komano et al., 2006). Notably, several transcriptomic and proteomic studies in PBMs have revealed significant insights into the pathogenic mechanisms of osteoporosis (Leung et al., 2011, Kotani et al., 2013, Zhou et al., 2015). On the other hand, we acknowledge that the ideal model cells for osteoporosis study are primary bone cells (e.g., osteoblasts, osteoclasts, and osteocytes). With the continuous development of high-throughput multi-omics profiling technologies, particularly for single-cell sequencing (Macaulay et al., 2017), we will be able to apply multi-omics analysis on human primary bone cells in the near future. Third, the identified biomarkers and results of causal analysis exclusively depend on computational modeling; hence, further experimental validation work should be conducted to confirm the biological significance and causality of these osteoporosis biomarkers. Nonetheless, we want to emphasize that traditional validation/further exploration using in vitro cells or in vivo mice models may be useful in some cases but may not completely reflect human in vivo functional mechanisms in other cases.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.