In the last decade, a series of large and well-powered studies have dramatically increased our appreciation of a multitude of genetic factors that influence skeletal diseases, including osteoporosis. Significant advances of the post-genomic era are expected to translate into enhanced ability to predict who is at risk for disease, and to enable better treatment of those who already have bone disease (1, 2). GWAS and whole genome sequencing (WGS) analyses have transformed the genetic analysis of complex diseases in general and osteoporosis in particular. The results of GWAS are increasingly being used by the pharmaceutical industry as an effective means of prioritizing compounds for development, as well as for repurposing existing medications for new indications (3).

Bone mineral density (BMD) remains the strongest predictor of fracture risk. As BMD also is highly heritable, it is frequently used in GWAS (4). The most significant study to date on the genetics of osteoporosis is a 2019 UK Biobank study involving approximately 420,000 participants (5). This study identified a total of 518 loci associated with estimated heel BMD, of which 301 were new loci. Of note, GWAS is mostly useful to identify common variants (usually defined as variants with a minor allele frequency >1%). On the other hand, genetic mutations are frequently discovered for less-common skeletal diseases; these mutations might be rare. A recent strategy that has already been employed in skeletal research is to use WGS, which is able to identify rare variants with large effect sizes. Such studies have identified several rare mutations in LGR4 (6) and COL1A2 (7) loci that are associated with low BMD. One particularly powerful study combining sequencing and GWAS identified a non-coding variant at EN1 (minor allele frequency = 1.6%) that also has large effects on BMD (8).

Interpreting results from human association studies remains a challenge, especially nominating causal genes for complex traits based on genome-wide significant SNPs (9). To bridge the gap between the genetic association and molecular function, a systematic functional investigation is necessary to interpret GWAS variants and to infer the exact disease-causing genes, or genes they regulate, and the cells in which they act (9). Here we review current practice for functional dissection of human complex traits and propose a roadmap for data integration and target prioritization for the skeletal outcomes.

Rare skeletal disorders span a broad clinical spectrum of bone-related pathologies, occasionally exhibiting extra-skeletal manifestations. Besides being genetically heterogeneous, the severity of these disorders is highly variable, ranging from neonatal lethality to minor complications discovered incidentally during adulthood (10, 11).

In contrast to complex traits having a multifactorial genetic etiology, genetic studies of monogenic diseases focus primarily on identifying the underlying causal rare variant(s) in affected patients, isolated or as members of a multiplex family (12). The first step consists of deep phenotyping of the clinical skeletal manifestations. The differential diagnosis based on the clinical and radiological observations might strongly indicate one candidate gene as an initial hypothesis explaining the underlying genetic causality. However, most cases have an uncertain genetic basis; necessitating hypothesis-free approaches. These include on the one hand homozygosity mapping and linkage analysis in multiplex families resulting in the delineation of a genomic region where the disease-causing gene resides. On the other hand, and potentially in combination with the previous approaches, high-throughput sequencing provides insight into the genetic variation within an individual. The widespread availability of recent -omics technologies has permitted researchers to focus their efforts on this approach by utilizing customized gene panels, whole exome sequencing (WES) or WGS. Nonetheless, all genetic discoveries resulting from traditional approaches such as linkage analysis or high-throughput technologies require translational assessment and annotation using in vitro or ex vivo bone cell work and/or in vivo knockout models to confirm disease association.

Examples of successful genetic findings with functional validation (i). Four consanguineous and distantly related individuals with autosomal recessive osteopetrosis were analyzed using homozygosity mapping (13). A single 1.22 Mb genomic region shared by all affected subjects was identified on chromosome 7, harboring five genes: NFE2L3, HNRNPA2B1, CBX3, SNX10, and SKAP2. Among these genes, SNX10 (sorting nexin 10) readily stands out as an excellent candidate due to its involvement in endosome homeostasis. A missense mutation was identified in all affected patients and was hence taken forward for functional investigation, whereby osteoclasts derived from monocytes of patients revealed gross abnormalities in the endocytic system and resorptive activity, abnormal SNX10 expression and altered subcellular localization of the encoded protein. Subsequently, Snx10 silencing experiments in mice highlighted the essential role of SNX10 in osteoclast vesicle trafficking and osteoclastic bone resorption (14, 15).

(ii) WES was conducted on three sisters. These sisters had a history of atypical femoral fractures after long-term bisphosphonate treatment for their underlying osteoporosis (16). WES analyses identified the presence of a rare missense mutation in GGPS1 (Geranylgeranyl Diphosphate Synthase 1), encoding the GGPPS enzyme, which acts downstream of the point of bisphosphonate action. Functional validation of the exact missense change, together with gene knockdown in osteoblasts and osteoclasts, was essential to confirm causality and to demonstrate the importance of the gene in atypical femoral fracture susceptibility (17). Additionally, other WES-prioritized variants, such as CYP1A1, were found mutated in other atypical femoral fracture cases (17, 18), opening the possibility of digenic or oligogenic inheritance. It might also reflect the idea that clinical variability, observed in many monogenic diseases, can be explained by variants in modifier genes (19). The discovery of such genetic variants opens an application window into personalized medicine (20).

(iii) Well known in the field, is the G171V missense variant in the gene encoding LRP5. The discovery of this variant was the outcome of a linkage study combined with a focused sequencing effort in a large family with several cases characterized by high bone mass (21). A combination of genetic and functional studies soon provided strong support for the involvement of LRP5 (22). Loss-of-function mutations in LRP5 explain the low bone mass in osteoporosis-pseudoglioma syndrome and other missense mutations in the same domain were identified in high bone mass phenotypes (23, 24). Soon afterwards, a wealth of in vitro and in vivo data confirmed the important role of the LRP5 gene in the regulation of bone mass (25), corroborated by GWAS and candidate gene association studies indicating the effect of a few common LRP5 variants on BMD and the risk for osteoporosis in the general population (26). Indeed, LRP5 is an empirical example of a gene that may harbor mutations or polymorphisms contributing to monogenic or complex forms of abnormal bone mass respectively. Identification of the defective gene in monogenic diseases serves as an optimal natural-occurring ‘knockout’ model, with population-based studies enabling gene prioritization and validation, disentangling the underlying pathogenesis in monogenic conditions. Other examples of genes and loci exist, discussed in more detail in this issue’s paper (12). Monogenic diseases are also not only constrained to mutations in the protein-coding regions. It has been shown the homozygous 52-kb deletion in the SOST-MEOX1 intergenic region on 17q12-q21 occurs in van Buchem disease patients (27). This region 35 kb downstream of the SOST gene fosters a long-range enhancer for it. Thus, the patients have reduced SOST transcription which reflects in lower sclerostin levels (28) ( Figure 1 ).

Benefits from recent advances in sequencing technologies: Novel techniques have remarkably facilitated the elucidation of the underlying molecular etiology of many rare skeletal diseases. Indeed, this has enabled the classification of conditions based on the implicated genetic defect and/or the altered metabolic/signaling pathways. This is why, the monogenic mutations can serve as human knockout models and help to uncover the gene function while GWAS findings serve to prioritize genes to scrutinize the cause of the monogenic conditions (12). In the case of osteogenesis imperfecta (OI), genetic discoveries have prompted a nosology revision of the existing classification, with causative genes added as new subtypes of the OI types I-V (29). A total of 17 genetic causes of OI have been described to date, with inheritance patterns ranging from autosomal dominant (e.g. COL1A1, COL1A2, IFITM5), autosomal recessive (e.g. LEPRE1, PPIB, SERPINH1, PLOD2, BMP1, WNT1), and X-linked (PLS3). An impressive 92% of 461 skeletal disorders described by the Nosology Committee of the International Skeletal Dysplasia Society have been solved at the genetic level thanks to high-throughput sequencing, creating new well-defined entities and sub-classifications of previously unknown or ill-defined skeletal disorders (11). Improved nosology based on careful clinical phenotyping coupled with genetic data leads to better patient care, both in terms of diagnosis and treatment (30). The discovery of causative genes and defective proteins has aided in the diagnosis, prognosis, and management of affected individuals, accelerating the development of personalized therapy. A good example is the treatment-option of bone marrow transplantation in patients with malignant recessive osteopetrosis. Unraveling the genetic cause in these patients before treatment decision is essential, as in RANKL deficiency, bone marrow transplantation will not have any beneficial effect (31–33). Finally, the identification of genes involved in monogenic diseases have resulted in novel treatments for osteoporosis, as is the case for denosumab, an anti-RANKL monoclonal antibody (34, 35) and romosozumab, an antibody against the Wnt-inhibitor sclerostin (36).

Fragility fractures represent a very complex phenotype. So far, most genetic studies have focused on BMD rather than bone fracture risk. There is a realization that the genes affecting BMD are not necessarily the same genes that influence fracture risk (37); however, there are no such indications from the fracture GWAS (38). There is the need for new phenotypes, which will enable and support the causal genes validation.

Imaging techniques like QCT, high resolution peripheral QCT, magnetic resonance imaging (MRI), and radiofrequency echographic multi spectrometry (39), together with fracture and BMD traits are considered as measurable “exophenotypes”, while “endophenotypes” - parameters that are more biologically proximal to gene actions - are currently lacking. Here we define the term, endophenotype (a.k.a. intermediate phenotype), as a characteristic able to mark genetic vulnerability independent of the clinical state (40, 41). Therefore, endophenotypes have the potential to identify the genetic dysfunction prior to disease manifestation. Similar to exophenotypes, the endophenotypes might be influenced by many genes, each with a relatively small effect, making endophenotype-identification difficult. Lifestyle factors (e.g. diet, physical activity) as well as other confounders can influence the exophenotypes such as BMD, QCT, fractures and others, and can mask the effects of genetic factors that we aim to assess in functional genomics (42). Hence, the main advantage of endophenotypes vs. exophenotypes is that their correlation with genetic changes is stronger, as they are more proximal to genes. The levels of molecules like proteins, metabolites, miRNA in bone cells, bone tissue, and/or in the blood can be chosen as the endophenotype.

Current status and needs in the field: It is desirable to have endophenotypes that can be used as specific biomarkers of bone cell activities in order to compensate for the shortcomings of BMD. In contrast to BMD, the potential serum/plasma bone biomarkers would ideally be able to reflect bone remodeling (43–46) in a dynamic fashion. Increased bone turnover results in microarchitectural deterioration of bone and has been associated with fracture risk independent of BMD. However, the evidence is currently not robust enough to use any biomarker in the fracture risk prediction tool (47).

Examples of established bone formation biomarkers used as endophenotypes in treatment monitoring are procollagen I N-terminal propeptides (PINP), bone-specific alkaline phosphatase, procollagen type I C-terminal propeptide (PICP) and osteocalcin, while C-terminal telopeptide of type I collagen (CTX), N-terminal telopeptide of type I collagen (NTX), tartrate-resistant acid phosphatase isoenzyme 5b, C-terminal crosslinked telopeptide of type I collagen, (ICTP), and deoxy-pyridinoline serve as resorption biomarkers (48). Bone biomarkers’ specificity for their respective process is convincing (48). However, there are major challenges even with the recommended reference markers CTX and PINP. Namely, especially CTX fluctuates during the day requiring blood samples to be collected from fasting patients in the morning, and both CTX and PINP vary tremendously among different individuals (49). Therefore, new biomarkers are being investigated, including proteins regulating bone resorption (RANKL, OPG), bone formation [sclerostin (50)] or bone non-collagenous proteins [periostin (51)].

miRNAs as endophenotype markers: Examples of potential new molecular biomarkers include non-coding RNAs (ncRNAs), of which miRNAs currently seem to be more promising (52, 53). They are small, 20–24 nucleotides long, noncoding, single-stranded RNA molecules that act as post-transcriptional regulators of gene expression. A cluster of miRNAs can target a single gene, and every single miRNA can regulate several different protein-coding genes. Their role in bone homeostasis is well established since miRNAs were shown to significantly affect the differentiation, proliferation, and function of both osteoblasts and osteoclasts (54, 55). Besides being intracellular, they are also present in several biological fluids where they are remarkably stable. Several studies have shown differences in circulating miRNA levels between osteoporotic and control subjects, both in primary and secondary osteoporosis [reviewed in (56)]. Based on these studies, it was proposed that circulating miRNAs could serve as a clinical tool for fracture risk prediction giving additional information on bone metabolism not captured by BMD, FRAX^®, or traditional bone turnover markers. A commercial test for fracture risk based on a panel of 19 miRNAs called OsteomiR™ was shown to effectively assess fracture risk (52). A cost utility model showed that its implementation could reduce fracture incidence compared with standard approaches such as monitoring BMD, no monitoring, or FRAX^® calculation alone (57). However, the key circulating miRNAs in osteoporosis are not consistent between studies, and before their implementation in routine clinical practice can become a reality, further studies are required to obtain validated disease-specific signatures (58, 59).

Endophenotypes as functional markers: Described endophenotypes are also relevant for in vitro functional validation of GWAS hits. Ideally, the endophenotypes should be able to reflect the “effects” of a particular genetic variation even with subtle changes in gene expression i.e. mRNA, protein level/activity, or other metabolites levels. In this context, identification and selection of endophenotypes depends on the position of the genetic variation of interest. If it is positioned in regulatory gene regions (e.g. the promoter), the gene expression and mRNA level is a well suited endophenotype. On the other hand, in the case of coding genetic variants, the protein should be qualitatively and quantitatively analyzed. Suggested approaches to functionally evaluate SNPs are presented in Table 1 .

In conclusion, endophenotypes that are needed for identification and evaluation of risk genes are important not just for progress in functional genomic research, but also for better prognosis and prevention of bone disease, which is ultimately the goal of the GEMSTONE consortium ( Figure 2 ).

With the boom of genetic testing, an opportunity for establishing possible causalities for the skeletal traits and diseases based on genomics was discovered. A seminal paper published in 2003 laid out the foundations and conveyed the key message of the method called “Mendelian randomization” (60). The main focus of this method is to establish the relationship between an exposure (a SNP) and an outcome (an endo- or exo-phenotype) via a “proxy”, an instrumental variable, that within a specific degree of certainty cannot be influenced by neither intrinsic nor extrinsic influences, i.e. confounders. Thus, the association between the exposure and outcome can indeed be established to arise from the two alone and not due to other factors. The results of such analyses may predict both the direction as well as the effect-magnitude. In terms of genomics, SNPs were proven to be a good fit for these instrumental variables. With time and technological advancements, large-scale GWAS provided well-powered and reproducible association results for risk prediction of common diseases (61), including the skeletal field (62). Many such studies in the field of skeletal diseases have since been performed, with a notable one scrutinizing the clinical risk factors of fracture (38). In contrast to low BMD, an established “causal” determinant of fracture, three is no evidence to suggest that increasing vitamin D (25-hydroxyvitamin D) levels in “sufficient” individuals will modify fracture risk.

Use of GWAS to identify quantitative trait loci (QTLs): The method of Mendelian randomization is not constrained to the results of GWAS alone. SNPs can also act as QTLs where they correlate with the expression of genes (eQTLs). Other QTLs include protein expression (pQTLs) (63, 64), DNA methylation (mQTLs) (65, 66), and chromatin acetylation and chromatin accessibility QTLs [reviewed by (9)]. eQTLs are abundant, with 48% of common genetic variants estimated to act as eQTLs for at least one gene (9).

Within such Mendelian randomization configuration, the effect of gene expression, as an exposure, can be tested for association against chosen traits, the outcome(s). The challenge that such studies face is that gene expression is highly tissue-specific; thus, expression in one tissue may not fully predict expression in another. This is especially important since it is unclear which cells are the true drivers of a disease (i.e., in which cell type GWAS variants act), as the pathophysiology of complex diseases often implicates interactions of multiple cell types (9). As of yet, there are not many studies that implement gene expression from bulk bone tissue. There are some that leveraged eQTL data obtained from whole blood and tested for effect on estimated BMD (eBMD) (5), and more recently eQTLs obtained from osteoclast-like cells derived from human peripheral blood mononuclear cells were tested for effect on the same eBMD trait (67, 68) as well as in the case of the newly reported osteomorphs (69) using a mouse model.

Co-localization: Co-localization analyses integrate eQTL and GWAS data (70, 71). Within their scope, the position of the (usually) topmost associated SNP(s) on each locus between the two datasets are analyzed, with results indicating whether the same SNP(s) drive both the gene expression at that particular region and the GWAS signal (e.g. where the SNPs effect on the GWAS trait is mediated by the gene expression). The difference with the Mendelian randomization-based approach is that the co-localization does not estimate the effect size and direction, but provides the probability of (a) shared causal variant(s).

Identification of candidate genes is more straightforward for coding variants, which may directly disrupt the structure of a protein (9). However, as early as 2012, it was realized that only a minority of GWAS hits fall within transcribed regions, with most of them mapping to introns (4.9% and 41.2% respectively (72). The leftover majority of GWAS hits thus cannot be easily linked to a candidate causal gene.

Also, in 2012, the ENCODE Project Consortium set out to map and describe functional elements encoded in the human genome across 1,640 data sets involving 147 different cell types, amongst which were also human osteoblasts (73). The mapping was expanded in 2015 by the Roadmap Epigenomics Consortium (74). The consortia assayed the available cells for eight histone modifications. By further integrating five specific histone modification marks (histone H3 lysine 4 trimethylation and monomethylation – H3K4me3 and H3K4me1 respectively; trimethylations of histone H3 lysine 36 – H3K36me3; histone H3 lysine 27 – H3K27me3; histone H3 lysine 9 – H3K9me3) they were able to build the 15-chromatin-state model (and later expand it to the 18-state model by inclusion of the histone H3 lysine 27 acetylation; Figure 3 ).

In short, by combining information on the methylation and acetylation dynamics, they were able to predict local chromatin states that were roughly divided into 8 active and 7 repressed states, now used to study the relationship between histone modification patterns, RNA expression levels, DNA methylation, and DNA accessibility. Their findings additionally showed tissue-specific enhancer regions, epigenomic dynamics during lineage specification, and both similarities and differences thereof between various tissue and cell types (74). Importantly, both consortia freely released their data repositories for use by others. This data can be integrated and tested against for enrichment e.g., by incorporating it into GARFIELD (GWAS Analysis of Regulatory of Functional Information Enrichment with LD correction) or by using the partitioned heritability function of the LDSC (linkage disequilibrium score) regression method (75, 76). This allows researchers to further functionally analyze and annotate their results and potentially discover novel cell and tissue specific genomic patterns.

One of the goals of in silico and in vitro methods is to translate findings to in vivo models. Historically, mouse – and more recently, zebrafish models - have been used to explore variants present in protein-coding regions. Recently it has been shown that despite poor evolutionary conservation in non-coding genome sequences, this model can still be used to compare enhancer activity of putative variants as predicted by in silico findings (77).

With WES being increasingly applied to large population-based settings, the American College of Medical Genetics and Genomics and the Association for Molecular Pathology (ACMG-AMP) have released standards and guidelines for the interpretation of sequence variants (78). According to these guidelines, the variants are classified as 1) benign, 2) likely benign, 3) uncertain significance, 4) likely/expected pathogenic, and 5) (known) highly pathogenic. Databases can follow this classification system, whereas others, such as Human Gene Mutation Database (HGMD), use their adaptation of functional classifications (79, 80). The Combined Annotation Dependent Depletion (CADD) (81) is a tool that uses a machine learning approach for scoring the deleteriousness of both coding as well as non-coding variants. It has been shown that the occurrence of known pathogenic and expected pathogenic variants in healthy populations are higher than expected based on disease prevalence (82–84). Many variants predicted to be potentially pathogenic have a lower than expected penetrance in healthy populations. To this effect, the UK BioBank has recently forayed into the venture of sequencing whole genomes of its participants, identifying rare variants to improve the prediction of monogenic and complex traits. The first tranche of results (n=49,960) was released to be used by the wider scientific community (85), with an addition of exome sequencing data on 150,000 volunteers added to the UK BioBank database (https://www.ukbiobank.ac.uk/2020/10/uk-biobank-makes-available-new-exome-sequencing-data).

As a proof of this concept, evaluation of WES data with clinical information has been done in a deeply phenotyped cohort study (86). They discovered 26 variant carriers, but only three of them experienced a clinical event related to the identified variant. When they consulted two main databases, ClinVar and HGMD, for clinical interpretation of the variants, they found a high degree of disagreement between the two databases. Moreover, the clinical classifications within ClinVar in different releases over five years (2014-2018) evidenced a trend of changing the clinical interpretation of formerly expected pathogenic variants towards class 1 (benign), 2 (likely benign) or 3 (uncertain significance). As shown, the definition of known pathogenic variants is ambiguous between databases, yet also differs between different versions of the same database. Moreover, potentially pathogenic variants do not always have a clinical impact. This presents challenges that researchers, as well as clinicians, face and must address while interpreting their findings.

Traditional GWAS using microchip analysis has been limited to less than half of the genome, since the major part consists of various repeated sequences and retrotransposons in which accurate localization of genomic variants is not possible. This may, however, change with a trend towards WGS, enabling very long reads spanning repeated sequences with the newest DNA sequencing tools. The most abundant retrotransposon, Long Interspersed Nuclear Element 1 (LINE1, L1), is 6 kb long and constitutes 20% of the genome (103). Of relevance to bone metabolism, a recent study showed that blocking L1 activity hinders differentiation of bone marrow mesenchymal stromal cells (BMSC) into osteoblasts, while transfection of BMSC from osteoporotic women with the L1 transcript stimulates osteoblast differentiation and bone production. In line with these results, the L1 copy number was increased in bone from healthy postmenopausal women as compared to osteoporotic women {Mangiavacchi, 2019 #1773}.

A decrease in CpG methylation status of repetitive sequences has been demonstrated during osteogenic differentiation of BMSCs suggesting that their methylation status is important for the induction of osteogenic differentiation (104). In the blood of postmenopausal osteoporotic women, the methylation status of Alu short intersperse elements (SINEs) has been associated with lower BMD, suggesting a positive correlation between Alu hypomethylation and age-related phenotype such as loss of BMD (105). In rare cases, when retrotransposition occurs in germ cells, unhindered by the host restriction mechanisms, novel mutations arise in the host genome that can be vertically transmitted to the next generation. Mutations induced by repetitive sequences were also reported in genes involved in musculoskeletal development. Mouse mutant chagun with skeletal dysplasia and male infertility has been demonstrated to have LINE-1 insertion in Poc1a gene (106) [reviewed in (107)]. Mutation in Poc1a (encoding protein of the centriole 1a) caused centrosome dysfunction which led to disorganized epiphyseal growth plates of their long bones [reviewed in (107)]. Seen from an evolutionary perspective, the youngest and most active retrotransposons are Alu SINEs comprising 11% of the human genome and unique to primates. Human endogenous retroviruses (HERVs) are a group of repetitive sequences that comprises 8% of the human genome (103). Regarding MSK physiology and pathology, HERV-W expression has been shown to increase in synovial fluids of human patients with osteoarthritis (108) and to play a role in human osteoclast fusion (109, 110).

Mobile genetic elements also shape the regulatory landscape of the human genome; for example, HERVs alone provide 320,000 transcription binding sites, Alu elements provide numerous splicing donor sites, whereas LINE-1s contribute to the generation of retrogenes and probably cause somatic mosaicism. Under normal conditions the majority of repetitive sequences are methylated, thus not transcribed. Under stress conditions and in pathological states (e.g. viral infections, inflammation, cancers) the expression status of mobile elements is altered and repetitive sequences are transcribed which can change the transcriptional regulation of genes and could lead to genetic instability (111). Even though mobile genetic elements form half of the human genome, their role in transcriptional regulation in skeletal diseases awaits further elucidation.

SNPs and other variants in non-coding regions can potentially influence the binding affinity of transcription factors and consequently change the regulation of bone homeostasis. Recently, several comprehensive functional studies of SNPs in intergenic regions have combined bioinformatics data analysis followed by functional validation in vitro and in animal models (112–115). A targeted search for novel potentially functional SNPs in enhancers that associate with bone metabolism was performed in five independent cohorts including 5,905 patients (116). In this study, correlation of SNPs with gene expression and biological processes resulted in 15 novel SNPs in enhancer regions (116). Analysis of transcriptional binding sites in the vicinity of osteoporosis associated SNPs revealed that they could affect the binding affinity of common transcription factors (NFATC2, MEF2C, SOX9, RUNX2, ESR2, FOXA1 and STAT3) which may be affected by SNPs and are involved in bone metabolism (117). High-throughput assays to speed the identification of functional variants have been developed in the last decade. Massive parallel reporter assays allow for the testing thousands of candidate regulatory sequences through cloning to a reporter gene followed by deep sequencing and analysis of transcription activation. This technique has recently been used to test 1605 SNPs residing in haplotypes implicated in osteoarthritis (118), highlighting its value to accelerate SNPs functional tests and genetic prioritization.

SNPs in the enhancer regions can change transcription factor binding sites (TFBS) and thus influence transcriptional regulation in bone metabolism. A study identified 5081 osteoporosis-related SNPs residing in enhancers (119). Transcription factor enrichment analyses identified EZH2 TFBS as a common binding site typical for osteoporosis associated enhancer SNPs (119). Comprehensive analysis combining integrative functional genomics and experimental validation methods was reported for functional assessment of osteoporosis risk locus on 1p36 (112). First, the authors prioritized a particular SNP (rs6426749) with functional genomics. They then confirmed with dual luciferase assays and CRISPR/Cas9 silencing that this SNP acts as an allele-specific enhancer regulating the expression of a lncRNA (LINC00339). The downregulation of LINC00339 increases the expression of an important regulator of skeletal development, CDC42 (112).

Moreover, Zhu et al. have performed a comprehensive analysis to explain associations between SNPs in a potential RANKL enhancer region located 100 kb upstream of the RANKL gene and risk for osteoporosis (113). They employed eQTL, high-throughput chromosome conformation capture (Hi-C), epigenetic annotation, and functional assays to show that several SNPs residing in non-coding regions exclusively correlated with RANKL expression. This study revealed that RANKL transcriptional regulation is mediated by a long-range super-enhancer.

The importance of long-range enhancers was also demonstrated for SOST. As already mentioned, in van Buchem disease, patients carry a homozygous 52-kb noncoding deletion that is essential for the transcriptional activation of SOST in the bone (27). Deletion of specific long-range regulatory element Ecr5 in mice caused the elevated bone formation and higher bone mass implying that the ECR5 long distant region is responsible for transcriptional activation of Sost in the adult skeleton (114). In an integrative study, Carey et al. searched for BMD associated SNPs that are enriched in lineage-specific pathways during osteoclast differentiation (115). An overlay between BMD GWASs and active enhancers in the myeloid compartment revealed that the PU.1 transcription factor network is important for osteoclast differentiation.

The above exemplifies that (a) identification of enhancers is important, and (b) the combination of GWASs with experiments on model organisms helps to decipher pathways for skeletal conditions.

Over the past few decades tremendous advances in -omics technologies (transcriptomics, epigenomics, proteomics and metabolomics) have greatly expanded our knowledge into the cellular and molecular diversity, and pathological mechanisms underlying many diseases, including those affecting mostly the skeleton like osteoporosis and other skeletal conditions (120). While each -omic technology possesses the potential to capture a snapshot in a cell’s lifetime or disease state, individually they lack the power and capacity to capture holistic and spatiotemporal changes that occur at both the cell and tissue level during disease pathogenesis and progression. Therefore, there is growing momentum towards concerted multi-omic studies in the effort to integrate and unify data from different -omics platforms and thus better encapsulate all of the multilevel molecular and functional pathways that underpin a particular disease. However, the unification of -omics data presents challenges in combining and interpreting multilevel data sets, which are inherently large, complex, and call for significant computational grunt coupled with high-end bioinformatics. There are advantages of single-omics based approaches described herein, each that have contributed significantly to our current understanding of bone cell function and skeletal disease.

Since bone is a relatively inaccessible tissue, few -omics data resources are available; more specifically, data related to chromatin/DNA structure are missing. A summary of -omics data resources originating from human bone tissue is presented in Table 2 .

To best exploit the power of the various -omics data, genetic alterations must be combined in order to understand the interaction between features like SNPs, chromatin structure, DNA methylation, coding transcripts, non-coding transcripts, metabolomics, and bone status/bone metabolism, ideally supplemented with functional studies in cells and model organisms. GWAS follow-up studies are necessary to interpret GWAS results and to infer the exact disease-causal variants, the genes they regulate, and the cells in which they act (9).

With Hi-C, the chromatin loops and topologically associating domains (TADs) can be mapped (236). The hierarchical organization of chromatin can be further detailed with Assay of Transposase Accessible Chromatin sequencing (ATAC-seq), which maps nucleosome-free DNA available for transcription (237, 238). Furthermore, DNA methylation, as well as genome variants, can influence binding of gene regulatory proteins, thereby regulating gene expression levels. The DNA methylation pattern is also associated with the three-dimensional structure of DNA (239). DNA accessibility peaks indicate regions available for transcription factor (TF) binding to histone modifications (e.g. H3K4me1, H3K4me3, H3K27ac, and H3K27me3) (240). In particular, H3K4me3 peaks highlight gene promoters while H3K27ac peaks mark active enhancer and promoter regions (241). Thus, to promote the understanding of the underlying molecular mechanism of bone metabolism, several -omics analyses should be performed on the same sample in cells from patients with osteoporosis and controls. Unfortunately, such comprehensive studies are still missing, and at best, analyses of two or three different -omics layers have been combined in the same study. Current searches are largely limited by the availability of comprehensive reference functional data sets and the emerging set of analytical tools for multi-omic analysis.

The various studies often have different designs and purposes, and therefore are not directly comparable. Sclerostin is a central inhibitor of the Wnt signaling pathway, and various parts of the Wnt signaling system have been associated with bone status in most types of -omics analyses including GWAS, e.g., β-Catenin and DKK1 in proteomics (187); SOST, DKK1, WIF1, CTNNB1, and WNT5B in transcriptomics (193, 196); FZD10, TBL1X, CSNK1E, WNT8A, CSNK1A1L, SFRP4, and SOST in DNA methylation studies (221, 223). Also, TGF-β signaling genes (187, 196, 198) and regulators of osteoclast function (187, 197, 199, 242) have been identified in more than one type of -omics analysis. Furthermore, a study of DNaseI hypersensitive sites during osteoblast differentiation identified changes in chromatin and expression of several genes within the Wnt and TGF-β signaling pathways (232).

Some GWAS studies have included eQTL results, but often transcript levels of genes near the resulting variants are neither associated with the allele frequency nor BMD (243). A recent study indicates that Hi-C type methods are well suited to identify the effectors of causal genomic variants (230). In this study, the chromatin capture technique was combined with ATAC-seq to map 46 BMD GWAS loci to 81 gene promoters in human MSC-derived osteoblasts. Consequently, several novel genes physically interacting at the three-dimensional (3D) genome level with the causal variants of BMD were identified.

From Table 2 it follows that gene expression is available for a wider set of cells and tissues than other types of -omics data (9). However, to get a comprehensive understanding of the genomic changes underlying bone diseases, it is necessary to identify eQTLs with different effect sizes at different stages of cell differentiation (dynamic eQTLs). This applies to cells of MSC progeny as well as those of monocyte lineage.

To date, there is not one single online resource that has collected data from all the available -omics analyses done on human bone tissue. At best, the resources are scattered throughout several such outlets. Sequence Read Archive (13), wherein high-throughput sequencing data is curated, is one of those (244). Even though the database itself is rather large, advanced search functionality built into it allows for easier navigation through the contents, which can be used to discover human bone tissue derived data. SRA is complemented by Gene Expression Omnibus (GEO) that may also host these same datasets, yet often enough, unique data can be found there as well (245). Since GEO is integrated into the National Center for Biotechnology Information (NCBI), like SRA is, its built-in search functionality allows for a similar navigation of the contents. ProteomeXchange is a portal dedicated to protein expression datasets (246), listed as one of the primary information resources by the Human Proteome Organization. Although human bone is a rare tissue in proteomic studies, bone-derived proteomic datasets are available on ProteomeXchange, although not in the same abundance as transcriptomic studies in other comparable resources.

Perhaps the closest approximation of “one-size-fits-all” collection of -omics datasets may be the OmicsDI platform, which acts as an integrational portal for proteomics, genomics, metabolomics, and transcriptomics datasets (247). The built-in search function in the portal has certain limitations e.g. improper implementation of the “NOT” operator in order to filter out undesirable results. The platform does though include a feature RESTful API for a possibility of implementing the search functionality within another website, or automating and curating the search results via a scripting language of choice, which may overcome the aforementioned limitations.

SkeletalVis is a portal devoted to exploration, visualization, linking of- and meta-analyzing skeletal transcriptomic data (248). Publicly available data resources (such as SRA, GEO, and ProteomeXchange mentioned before) are mined, undergo a QC procedure, and re– analyzed. The strengths of the platform include inter-experiment comparison options using signed Jaccard index, an approach that is also species-agnostic, visualizing datasets’ (dis)similarity using the t-distributed stochastic neighboring embedding, as well as other possible downstream analyses, whilst presenting the results in a user-friendly web interface.

Another user-friendly tool is a correlation browser to identify highly correlated transcripts in trans-iliac bone biopsies from 84 postmenopausal women (193). The correlation browser enables targeted searches among >260 million transcript correlations. This tool (http://app.uio.no/med/klinmed/correlation-browser/iliac-v2.0/) enables e.g. identification of candidate targets of transcription factors. It has been expanded to include mature miRNAs, thus also enabling identification of candidate interacting mRNAs/miRNAs (unpublished).

The flourishing of analytic in silico tools and software is remarkable, and increases the speed at which data can be processed and analyzed. However, with this abundance of possibilities, caution is warranted, as no single tool is comprehensive and none is infallible. It is imperative to understand the principles behind bioinformatic tools and to sensibly choose the most suitable one(s) for the purposes of the end user’s project(s).

Cellular models that accurately resemble/reflect the morphology and physiology of their originating/native tissue are pivotal tools to study bone biology and disease. However, the isolation of homogenous and functional primary bone cell populations remains technically challenging as most cellular residents are bound tightly to bone surfaces (i.e. osteoclasts and osteoblasts) or deeply entrenched in mineralized tissue (i.e. osteocytes), thus requiring specialized isolation methods. An amalgamated cellular repository: i) composed of a wide variety of primary and transformed skeletal cells; ii) derived from relevant skeletal tissue types, and iii) of multiple species of origin would therefore facilitate rapid and transparent inter-institutional exchange of cellular resources and methods. To this end, a survey of cellular resources within reach of this consortium (from 17 research teams) is illustrated in Figure 4 . The scheme encapsulates both primary and immortalized cell lines, the majority being human or mouse origin (N=79 and N=55, respectively) but also includes rabbit, rat, and monkey. Multipotential MSCs derived from primary tissue sources are best represented (N=46 human) and (N=39 mouse) followed by tumor-derived osteoblast cell lines (N=18) and those of myeloid lineage/PBMCs (N=4 human, N=16 mouse). Osteoblasts are by far the most represented cell type (N=71) followed by osteoclasts derived from either human (N=13) or mouse origin (N=16), respectively. Not surprisingly, osteocytes are comparatively underrepresented (N=11), with only four derived from primary human sources. This cellular resource is not limited to bone cells, but also extends to neighboring and intersecting tissues/cell types.

Collectively, the shared cellular models will serve as a powerful resource towards accelerating the functional validation of new molecular targets potentially implicated with skeletal biology and disease. However, although bone cell cultures are an easy and useful tool for the discovery and/or functional validation of variants associated with bone cell differentiation, they do not fully reflect the in vivo situation. Complex interactions between cells and the surrounding matrix are often missing. Despite these drawbacks, cell cultures have several benefits including analysis on a specific cell type, excluding the influence of endocrine factors and complex tissue interactions, ease of gene manipulation, and they are cost effective. It can be seen from Figure 4 that roughly half of the cellular resources used are immortalized/transformed cells. These are especially helpful because the cells in culture are rather uniform and can more readily be genetically modified, in contrast to many primary cell cultures. However, immortalized cells should be used with caution since the cell cycle is artificially altered through the transformation that may potentially affect cellular signaling. Where possible, confirmation in primary cell cultures could be of benefit.

We find that the most important step needed to make progress on the functional validation of GWAS findings using cell cultures, is that experts throughout Europe and the rest of the World share cells, protocols, and expertise with each other (249). Resource from Figure 4 is a first step in this direction, but the GEMSTONE consortium will further substantiate this collaboration within the network, thus we strongly encourage other research groups to join us in this effort.

To date, most of the data generated in laboratory settings are either 2D in vitro culture systems or animal models. Cell cultures often involve human cell types, whereby one or more cell lines are (co)cultured, which nonetheless lack some of the complexity as observed in real life. Hence, there is a pressing need for models that better reflect human bone metabolism, in which certain aspects can be included, such as a 3D environment, shear stress and chemotaxis. Approaches that have been taken to overcome these issues can be divided into 3 categories: 1) 3D printed bone scaffolds, 2) Bioprinting of scaffolds (containing cells), and 3) Microfluidic models (Organ-on-Chip).

3D printed bone scaffolds have been employed with numerous compositions, surface modifications, coatings, biomechanical properties, and porosities (250, 251). Initially, by implementing MSCs/osteoblasts to study osteogenesis, the complexity has gradually been increased by including e.g. endothelial cells/vasculature leading to vascularized bone tissue-engineered constructs (252). Taking this one step further, efforts in the leukaemia/bone tumor research field have led to models (partially) mimicking the human bone marrow microenvironment, enabling the study of complex pathology through simultaneous interactions between multiple cell types and their extracellular environment in vitro (253). More recently, 3D bioprinting tools have become available, which allow for generating a 3D structure with the cell type(s) of interest included in the printing process (254, 255). 3D scaffolds have been widely used as in vivo bone regeneration models, but their translation value for human bone biology is yet unclear.

Over the last decade, simple microfluidic set ups have evolved into multi compartment-based chips, often coined Organ-on-Chip (OoC), in which relevant physiological aspects can be studied, including shear stress and chemotaxis. Many OoC models have studied shear stress for endothelial cell function, but evidence is growing that also MSCs, osteoblasts and osteocytes perform better under fluid flow as evidenced by increased proliferation and altered marker genes expression (256–258). This suggests that OoC models may better reflect biology than conventional cultures. Similar to 3D printed scaffolds, OoC approaches have also led to employing more complex microenvironments, for example breast cancer-derived bone metastases or the so-called Bone Marrow-on-a-Chip (259, 260). The small format of OoCs may also allow for personalized medicine initiatives and for compound screening, as small amounts of cell numbers and materials suffice for cell culturing (261, 262).

Thus, despite the technical challenges ahead of us, 3D (bio)printing and microfluidics are at the forefront of a new era that may enable us to better recapitulate the physiology of bone tissue. The outcomes from various bone-related GWAS and clues from monogenic disease states has yielded a valuable list of target genes to scrutinize in a 3D environment with all the relevant physiological cues. With the use of primary cells or induced pluripotent stem (iPS) cells, the toolbox expands to generate a ‘bone-on-a-chip’ that relates to the disease or condition of interest. Ultimately, this may lead to improved therapeutic opportunities for bone metabolism and tissue engineering.

A cornerstone of science is the ability to replicate results. In the early 2010s the “replication crisis” was formulated to drive attention to the problem of inability to reproduce many findings (355). Part of this issue is also an increasing volume of scientific research being published (356). As part of the scientific method involves creating a conjecture, which relies on observations and prior knowledge, the mentioned crisis presents a problem that is two-fold: first, even though the information age makes it ever so simple to allocate desired material, the direct cause is also the exponential growth of the amount of data itself (357). Second, even when one is able to find wanted knowledge, the replication crisis should make researchers always second-guess the published results on their journey towards hypothesis creation. Thus, it would be prudent to build upon a service such as SkeletalVis and the “Musculoskeletal Knowledge Portal” [MSK-KP (358)], that allow not only -omics data resources to be available, but would also offer an overview, a “curation”, of results, whilst enriching them with other sources of relevant information. It is only through community effort that such advances are possible.

International Federation of MSK Research Societies (IFMRS) realized that there is a growing urgency for reproducible research using integrated -omics, similar to all disciplines in science. Along with the creators of relevant databases, data archives and knowledge-databases, IFMRS strived to identify common issues in database development, curation and management, to create data portals allowing reproducibility of singular -omics and integrated omics data (358).

IFMRS recently created MSK-KP, designed to facilitate functional studies of the genetic factors underpinning skeletal diseases, thus supporting the prioritization of genes and pathways by experimentalists. The ultimate vision of the MSK-KP was to consolidate -omics datasets from human and model organisms into a central repository that can be accessed by researchers to better understand the biological mechanisms underlying musculoskeletal disease and apply this knowledge to identify and develop new disease interventions. To realize this vision, the team from the Broad Institute (Cambridge, Massachusetts) was recruited, who had been instrumental in designing the Knowledge Portals for other complex diseases. Rather than simply serving as a repository for skeletal datasets generated by individual laboratories and large consortia, MSK-KP will provide a much-needed bridge between the statistical genetic, wet lab and clinical communities (358).

This effort requires making -omics data a resource on a web server that is publicly available to the scientific community, via an intuitive and flexible web interface that enables non-specialist users to mine and interpret -omics data easily (359). Thus, the summary of the results of existing and ongoing GWAS and PheWAS analyses is already deposited at the portal. At present, transcriptomics and epigenomics data are populated there. At the next stage, proteomics and metabolomics datasets will be added; in the future, bone-related lipidomics, microbiomics, spatial transcriptomics, and phenomics data will follow.

MSK-KP group will continue identifying, obtaining, curating and integrating various -omics datasets from the international MSK research community, and encourage data sharing through community collaborative spirit (358). Together with hosting the data from cellular experiments and animal models, as well as compound screens, the portal is supposed to integrate such data with bioinformatics resources. Part of the solution of appropriate integrating -omics datasets is a requirement for open sharing of scripts and codes for such analyses. Data provided to the MSK-KP should adhere to Findability, Accessibility, Interoperability, and Reusability (FAIR) principles (360). Besides archiving the data, IFMRS strives to provide guidelines on the integration of -omics datasets for the development of standardized analytical pipelines. This can move the skeletal genetics into the “post GWAS” era.

This consensus statement aims to create a roadmap for using functional genomics to interrogate signals from human genetic studies for osteoporosis and other skeletal conditions. While the number of -omics-based studies in bone biology has exploded over the past decade, high throughput in vitro systems and rapid phenotyping of model organisms remain equally important in order to accelerate the functional validation of identified targets. Our effort nurtures such endeavors by expanding the wealth of knowledge and resources, represented by keen individuals and assets at their disposal, and promoting their exchange amongst participating institutions.

The unification of -omics data will create a wealth of new information towards the study of skeletal diseases (359). Besides increasing sample size, algorithms for standardization of sampling and sample processing, platforms for quantification, and data analyses should be agreed upon. Bone cell-type-specific resources might allow high-throughput massive parallel reporter assays that test the variants altering the activity of putative regulatory elements (361). Advances in computational analyses together with novel gene editing techniques that enable epigenetic manipulations at particular enhancer sites will further unveil relationships between non-coding SNPs and disease development.

Functional genomics studies of osteoporosis have limitations, mostly due to inaccessibility of bone tissue. In particular, not many resources are dedicated to osteocytes; since longitudinal studies of (especially human) specimens are limited, studies of bone loss rate are underrepresented. There are also limited ATAC-Seq and Hi-C datasets derived from bone cells, and it is still challenging to determine protein levels in bone in an in vivo setting, especially for diagnostic purposes. Furthermore, even though mobile genetic elements form half of the human genome, their role in transcriptional regulation in MSK diseases awaits further elucidation. While still vastly underrepresented compared to other -omics technologies, metabolomic studies of resident bone cells will undoubtedly uncover hitherto unappreciated but important metabolites. All this dictates our joint interest for the near future.

As we start to unravel the proteomes of bone cells and build towards more complete protein atlases, future approaches should be aimed at reducing sample complexity. Reducing ‘noise’ will thus enable precision mapping of proteomes at single organelle resolution, such as mitochondria or lysosomes. Similarly, rapidly evolving multiplex technologies such as imaging mass cytometry (e.g. Hyperion) and imaging mass spectrometry (e.g. Cell DIVE) coupled with spatial transcriptomic platforms will undoubtedly afford new and exciting opportunities to systematically map the molecular, cellular and spatial organization of musculoskeletal tissues at unprecedented resolution.

There is a need to continue close communication with the animal modeling groups to assure that in vivo studies cross-fertilize with cellular and translational ones, allowing the planning for functional validation of candidate variants to be taken concomitantly with genetic prioritization. Alongside the plethora of tools available for functional studies in cell culture and mice, zebrafish provide alternative and advantageous solutions for rapid functional screening of coding and non-coding sequences, and translation of genomic findings to therapeutics. The combination of functional expertise made accessible through our collaborative group, allows us to discuss data access, curation, and sharing in a collaborative spirit towards non-overlapping efforts, directing funding resources towards the boosting of skeletal genetic discovery.

Furthermore, modeling of human disease does not only require the standardization of both cellular and animal phenotyping (42), but oftentimes also rethinking of the disease definition. Employing better defined heritable traits (endophenotypes) would benefit both the etiological understanding of the disease in a particular patient, improve the targeted therapeutic approach with fewer side effects, and provide more effective treatments. Moreover, gender-balanced data are needed as most of the -omics studies to date were focused on women.

As we move closer towards an integrative multi-omics and holistic approach to skeletal diseases, consortia such as GEMSTONE become a fundamental tenet of modern genomics. In the long term, these efforts will allow meaningful and physiologically relevant data to be extrapolated, will allow identification of molecular diagnostic biomarkers and translation of findings into new therapeutic targets with higher effectiveness and fewer adverse effects, which will contribute to a higher quality of care for human skeletal diseases.

MR, FR and DK initiated and organised the manuscript. EK, VP, BBa, NAL, IS and KS generated the figures for the manuscript. EK, VP, BBa, SR, JM and DK generated the tables for the manuscript. All authors contributed to the writing and approved the final manuscript.

Funding was obtained from the GEMSTONE COST Action (CA18139).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.