Identification and characterization of BMD Genes

Bone mineral density data on 3,823 genes were recorded and analyzed (Fig 1A). Within this gene set, we identified 200 (5%) genes that caused a statistically significant (p<0.0001) decrease or increase in BMD when compared to wt control animals (Fig 2A and S1 and S2 Tables. Together, the pool of 200 BMD genes was then analyzed further. A majority, 123 genes, produced a low BMD phenotype, while 64 genes caused a high BMD phenotype (Fig 2B). In addition, 13 genes were associated with either a low or high BMD phenotype, depending on sex (Figs 2B and S1). Thus, the 200 BMD genes were derived from 213 mutant mouse lines, of which 161, 47, and 5 were homozygotes, heterozygotes, and hemizygotes, respectively (S1 Table). The identity of the 200 BMD genes is shown in Table 1 and S1 Table. With the exception of 4930591A17Rik, biological functions of all BMD genes had been previously reported to some extent. To determine if novel BMD loci were among the 200 genes, we discriminated between genes with a known (class 1) and unknown (class 2) function in bone biology based on previously published data. A total of 59 and 141 genes were assigned to classes 1 (blue) and 2 (yellow), respectively (Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Mouse phenotyping identified 200 BMD genes.

(A) Of the 3,823 genes deleted in mice, 200 or 5% caused a change in BMD phenotype. (B) Directional changes in BMD. Low and high BMD was found after deletion of 123 and 64 genes, respectively. Another 13 genes caused low or high BMD dependent on sex. (C) Percentage change in BMD compared to wt animals. In both male (blue bars) and female (pink bars) mice, most BMD genes reduced BMD between 5% and 10%. Changes between -2% and +2% were not plotted as these were not considered by PhenStat as significant phenotypes. (D) and (E) A limited weight difference compared to wt controls in combination with a greater than 10% decrease or increase in BMD change prioritized (gray shaded areas) genes for further investigation. Genes with a greater than 10% decrease (D) or increase (E) in BMD were plotted separately, and genes affecting female or male mice represented individually. Priority genes in the gray shaded areas were colored yellow or blue according to class 1 or 2, respectively. For clarity, outside the gray sectors, only genes associated with greater than 20% decrease (D) or increases (E) in BMD were name tagged.

https://doi.org/10.1371/journal.pgen.1009190.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Identity and classification of the 200 genes controlling BMD.

https://doi.org/10.1371/journal.pgen.1009190.t001

To quantitatively assess the change in BMD caused by the 200 genes, we plotted gene numbers against percentage of change in BMD (Fig 2C). While decreases in BMD ranged from 2% to 49%, and increases from 2% to 52%, we observed that the majority of genes elicited a reduction in BMD between 5% and 10%. PhenStat, the statistical package used for the analysis of the high-throughput phenotype data, did not consider BMD changes between -2% and +2% to be significant phenotypes. A comparison between male and female animals showed an asymmetrical data distribution with about twice as many genes causing a 5% to 10% increased BMD in males than in females (Fig 2C). Regardless of sex, BMD decreased greater than 20% only upon deletion for the class 1 genes Lepr, Bbs5, Cyp27b1, and Ghrh, while it increased by more than 20% only after loss of the class 1 gene Lrrk1 (Fig 2D and 2E and S1 Table). To further prioritize the BMD genes, we set criteria by asking which class 2 gene deletions resulted in a greater than 10% decrease or increase in BMD in conjunction with a limited weight change. Weight changes can confound BMD assessments and we set an arbitrary limit to approximately 10% of the body weight of the animal. The 6 low BMD genes Dhx40, Fancl, Nhlh2, Prdm14, Slc9a4, and Spata22 (Fig 2D), and the 3 high BMD genes Cd5l, Dppa1, and Vat1l (Fig 2E) met our criteria. Deletion of Nhlh2 and Dppa1 yielded the most pronounced decrease and increase in BMD, respectively.

To survey for a potential role of the 200 BMD genes in human BMD, GWAS data based on the UK Biobank release was interrogated. We selected the UK Biobank because we believe the size of the data set is an advantage that offsets the difference in BMD acquisition technique, which was based on ultrasound and not the X-ray mean used in this study. A distance range between 0 kilo bases (kB) and 250 kb was chosen. Of the 200 genes, 11, 38, and 52 genes matched the human GWAS derived genes at a 0 kB, 100 kB, and 250 kB distance range, respectively (S3 Table). To take a closer look at genes with potential human relevance, we then utilized the GEFOS data set, which although smaller than the UK Biobank release, utilized X-ray-based BMD acquisition. The six class 2 genes AP4E1, KIAA0825, ACSF2, PKP4, TTC28, PHF19 and the class 1 gene MACROD2 altered BMD in humans (S4 Table). The direction of BMD change matched between mouse and human for the five genes ACSF2, PKP4, TTC28, PHF19 and MACROD2, and as described below the latter was found to be expressed in osteoblasts.

Sexual dimorphism of the BMD genes

The observation that the direction of BMD change depended on sex in 13 genes prompted further investigation of a potential sexual bias in BMD. Of the 13 genes, 11 genes caused low BMD in females and high BMD in males, and 2 genes produced high BMD in females and low BMD in males (Fig 3A). We then examined genes causing a unidirectional change in BMD. Twenty-six, 29, and 68 genes caused low BMD in female, male and both sexes, respectively (Fig 3B). In comparison, 6, 29, and 29 genes caused high BMD in females, males, and both sexes, respectively (Fig 3C). A cross-comparison independent of the direction of the BMD change, demonstrated that a sub-set of 90 genes regulated BMD exclusively in either females or males. Of those 90 genes, only 11 genes have been previously described in the context of sexual dimorphism in organs other than bone (Fig 3D).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. BMD genes displayed sexual dimorphism.

(A-C) A subset of 90 genes affected BMD specifically in one sex. The 90 genes encompass in females (pale pink bars) 26 low BMD and 6 high BMD genes plus in males (magenta) 29 low BMD and 29 high BMD genes. (D) The identity of all sexually dimorphic genes grouped by class and direction of BMD change. Genes with a previously reported role in sexual dimorphism other than bone are marked red.

https://doi.org/10.1371/journal.pgen.1009190.g003

Pleiotropy of the BMD genes

Taking advantage of the broad phenotyping approach of the IMPC, we screened the 200 BMD genes for additional phenotypes. Alterations in homeostasis/metabolism occurred most frequently (99 genes), followed by changes in growth/size/body region (94 genes), behavior/neurological (92 genes), and adipose tissue (72 genes) (Fig 4A and S5 Table). Next, we determined the inter-relationships of the co-phenotypes. With the exception of expected relationships, for example between homeostasis/metabolism, adipose tissue and growth/size, we observed a mostly uniform distribution of phenotype connections (Fig 4A). In contrast, there was variability in the pleiotropy distribution across the BMD genes. About half of the 200 BMD genes presented with up to 3 co-phenotypes, and only a few genes associated with large numbers of co-phenotypes (Fig 4B). Importantly, no co-phenotypes were detected for the eight class 2 genes Ang6, Dnase2b, Frmd7, Mlec, Pfkfb3, Stard5, Tmem42, and Zfp704, suggesting they affected the skeleton specifically. To assess if loss of any of these eight genes or any other class 2 genes produced gross skeletal abnormalities in addition to their BMD phenotype, we read the radiographs routinely performed as part of the IMPC phenotyping. None of the eight class 2 genes had both BMD changes and skeletal abnormalities. However, 19 other class 2 genes caused skeletal abnormalities in addition to changes in BMD (Fig 4C). Vertebrae abnormalities were the primary skeletal phenotype associated among the BMD genes (Fig 4C) and were seen in the 9 genes Adnp2, Arhgef4, Dscc1, Fam160a1, Hbs1l, Pex3, Sh3gl2, Slc25a30, and Slc38a10. As a combination of altered BMD and skeletal abnormality has the potential to compromise the biomechanical stability of bone, we explored the biomechanical testing data collected as part of the OBCD project. For this, we facilitated the import of skeletal phenotyping data from the OBCD into IMPC. Then, we surveyed the bones of the 19 knockout mice for compromised mechanical stability, a functional hallmark of osteoporosis. Of the 19 genes with both a BMD and skeletal phenotypes, the three genes Hbs1l, Slc25a30 and Slc38a10 matched to the OBCD data, but only loss of Hbsl affected skeletal stability, specifically it resulted in a decrease in vertebrae stiffness and yield load (Fig 4D–4F). This demonstrated a functional role of Hsb1l in skeletal maintenance and was consistent with the vertebrae abnormalities detected by radiographic screening (Fig 4G). Interestingly, loss of Hsb1l also caused cranio-facial abnormalities (Fig 4G). With the limited overlap between IMPC and OBCD in mind, we pursued another strategy to isolate good candidate genes. A weight change limited to about 10% and low number of maximal 2 co-phenotypes were used to identify genes more likely to directly affect the skeleton. The five genes Arl4d, Bbx, Fam160a1, Sdsl, and Tnfaip1 fulfilled these two criteria (Fig 4H).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. BMD genes were pleiotropic.

(A) Chord diagram plotting the inter-connectivity of co-phenotypes associated with the 200 BMD genes. The outer segments represent the co-phenotypes grouped as top-level Mammalian Phenotype (MP) terms. The segment size correlates with the number of genes exhibiting the co-phenotype. Homeostasis/metabolism was the primary co-phenotype detected. The size of the arcs connecting the phenotypes was proportional to the number of genes connected. Diagram colors were assigned clockwise starting with the primary phenotype. (B) The majority of BMD genes presented with up to 3 co-phenotypes, while eight genes were limited to a skeletal phenotype. Class 1 (known function in bone biology) and 2 (unknown function in bone biology) genes are distinguished in yellow and blue color, respectively. (C) Analysis of the IMPC phenotyping data showed that 19 of the 141 class 2 genes were also associated with skeletal abnormalities in addition to changes in BMD. Vertebrae abnormalities were the most frequent skeletal abnormalities associated with the deletion of one of the 19 genes. (D) OBCD data showed a difference in vertebra stiffness for Hbs1l compared to wt. (E) OBCD data detected a difference in vertebra yield load for Hbs1l compared to wt. (F) OBCD data measured a difference in vertebral height for Hbs1l compared to wt. (G) Dorsal images were taken from the IMPC x-ray screen. Representative female Hbs1l-deficient mice and wild type (wt) control mice are shown. In the upper panels, the position of the pubic body was denoted by a black line. In Hbs1l^-/- mice, a distal shift of the pubic body intersection relative to the S1 sacral vertebra was seen (p<0.0001 for vertebrae abnormality in the Hbs1l^-/- cohort versus wt mice). The two lower panels show a skull close-up on a second set of female animals. Identical outlines of the cranium were superimposed. Hbs1l^-/- mice showed a compressed cranium including nasal bones (p<1x10^-8 for vertebrae abnormality in the Hbs1l^-/- cohort versus wt mice). (H) Among the 19 genes with both BMD change and skeletal abnormality, a limited weight difference compared to wt controls and a low number of co-phenotypes (gray shaded area) indicated genes more likely to directly affect the skeleton. All 19 genes were plotted, genes in the gray shaded areas were colored in yellow according to their class 2 designation.

https://doi.org/10.1371/journal.pgen.1009190.g004

Function of novel BMD Genes in osteoclasts and osteoblasts

To offer insight into a potential direct role of the 200 BMD genes in bone we focused on the two principle bone cell types. Using available gene expression data, the 200 genes were assigned to osteoclasts and osteoblasts (Fig 5A and 5B). We found that of the 200 genes 29 and 16 genes were expressed in osteoclasts and osteoblasts, respectively (p<0.05), with an overlap in expression seen for Cbk, Ptafr, Satb1, and Sdsl. With respect to the five class 2 genes Arl4d, Bbx, Fam160a1, Sdsl and Tnfaip1 identified above, Arl4d and Sdsl were expressed in bone cells. Using STRING, we then tested (p<2.2x10^-16) for protein-protein interactions (PPIs) within the sets of 29 and 16 genes expressed in osteoclasts and osteoblasts, respectively (Fig 5C and 5D). Upon thresholding, 26 osteoclast and 13 osteoblast genes remained for construction of PPI networks. For several reasons we found such networks informative. Firstly, they identified proteins with a large number of interactions and hence hubs that route information. In osteoclasts, the class 1 gene product Ctnnb1 (ß-catenin) formed a central PPI hub (Fig 5C). In addition, the two class 2 genes products Pptc7 and Gtf2a1 served as PPI hubs in these cells. In contrast, in osteoblasts we observed that the class 1 gene product Col1a2 and to a lesser degree the class 1 gene product Plcg2 formed the central PPI hubs (Fig 5D). Secondly, patterns of PPIs uncovered prominent pathways in cells. For example, in osteoclasts we discovered a PPI chain connecting Pptc7, Excoc2, Rab3ip, Arf4, and Ncald (Fig 5C). The latter four proteins, participate in intracellular vesicle trafficking, highlighting the role of this pathway in BMD control via osteoclast function. Thirdly, PPI networks detected previously unknown links in BMD pathways. For instance, in osteoblasts (Fig 5D) the class 2 protein Arl4d interacted with Col1a2 via Kdelr3. Together, these findings support further investigation of the class 2 genes Pptc7, Rab3ip, Ncald, and Arl4d, of which the latter three produced a low BMD phenotype upon depletion (Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Subsets of BMD proteins mapped specific pathways in osteoclasts and osteoblasts.

Using curated NCBI Gene Expression Omnibus data, expression of the 200 BMD genes was mapped in (A) osteoclasts, the bone resorbing cells, and (B) osteoblasts, the bone-forming cells. Protein-protein interactions (PPIs) were identified using the STRING database allowing the annotation of 178 of the 200 BMD genes. The central hubs of the network where identified based on the betweenness centrality index and extracted from the main network. All PPIs were presented as black lines. Blue and yellow nodes represent class 1 (known function in bone biology) and class 2 genes (unknown function in bone biology), respectively. The zygosity of the mutants is indicated by a triangular (homozygous) or diamond (heterozygous) shape of the node. The resulting central protein-protein interactions sub-networks in osteoclasts (C) and osteoblasts (D) are shown. Proteins connecting the BMD gene nodes were represented as small pink nodes.

https://doi.org/10.1371/journal.pgen.1009190.g005

The finding that class 2 genes causing low BMD, such as Arl4d, Ncald, and Rab3ip, participated in osteoclast or osteoblast pathways, raised the possibility that they are involved in the skeletal mechanisms crucial for maintaining bone mass, particularly bone turnover. In the absence of established methods that permit rapid testing of potential bone turnover genes, without the need for animal experimentation, we devised a novel theoretical approach. As bone resorption and formation underlying bone turnover occur in a spatial and temporal fashion in the extracellular space, we exploited the known key proteins responsible for osteoclastic resorption of mineralized bone matrix and subsequent osteoblastic matrix formation and mineralization. The PPIs between these proteins were recorded and created an in silico bone turnover model (Fig 6A). We then introduced the three low BMD class 2 genes products Rab3ip, Ncald, and Arl4d, and found all three genes affiliated with essential extracellular bone turnover proteins in each case via a single intermediate. The model showed that during extracellular bone resorption, Rab3ip strongly interacted via Racgap1 with Mmp14 but also with Mmp9 and Atp6v0d2, suggesting an effect on both bone matrix and mineral degradation (Fig 6B). In contrast, Ncald interacted via Fgf2 only with Mmp9 and Mmp14 (Fig 6B), indicating a sole effect on matrix degradation. With respect to bone formation, Arl4d interacted via Kdelr3 with Pcolce (Fig 6B), possibly impacting bone matrix assembly. This further supported a role of Arl4d in osteoblast-mediated BMD maintenance.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. An in silico bone turnover model validated class 2 candidate BMD genes.

(A) The in silico model of bone turnover captured the extracellular break down of bone and the subsequent formation of new bone. The model is simplified excluding for example many of the intracellular proteins driving bone turnover or non-collagenous matrix proteins. While resorption of both bone matrix and mineral may occur concurrent, bone matrix formation typically precedes mineral deposition. The genes mainly responsible for bone turnover were grouped and annotated as follows. Osteoclastic bone matrix resorption: tartrate-resistant acid phosphatase type 5 (Acp5), essential for osteoclast motility; cathepsin K (Ctsk), matrix break down; matrix metallopeptidase 9, (Mmp9), matrix break down; matrix metallopeptidase 14 (Mmp14), matrix break down. Osteoclastic bone mineral resorption: ATPase H+ transporting V0 subunit D2 (Atp6v0d2), acidic pH and mineral removal; ATPase H+ transporting V0 subunit A3, (Tcirg1), acidic pH and mineral removal. Osteoblastic bone matrix formation: ADAM metallopeptidase with thrombospondin type 1 motif 2. (Adamts2), N-terminal processing of procollagen; ADAM metallopeptidase with thrombospondin type 1 motif 14 (Adamts14), N-terminal processing of procollagen; bone morphogenic protein 1 (Bmp1), C-terminal processing of procollagen; collagen type I alpha 1 chain (Col1a1) and collagen type I alpha 2 chain (Col1a2), procollagen formation; procollagen C-endopeptidase enhancer 1 (Pcolce), enhancer of C-terminal processing of procollagen; tolloid like 1 (Tll1), C-terminal processing of procollagen. Osteoblastic bone mineral deposition: tissue-nonspecific alkaline phosphatase (Alpl), removal of mineralization inhibitor inorganic pyrophosphate (PP_i); ectonucleotide pyrophosphatase/phosphodiesterase 1 (Enpp1), production of PP_i. Together, these proteins (red elliptic nodes) formed the PPI networks for osteoclastic resorption and osteoblastic formation during bone turnover. Connecting proteins distinct from the bone turnover genes are represented as pink elliptic nodes, while white nodes marked genes with an undefined role in bone. Black lines represent PPIs as determined by STRING or REACTOME analysis. Arrows indicated REACTOME-curated data. For validation, we used experimental gene expression data from osteoclasts and osteoblasts in combination with statistical testing. The resulting gaussian probability of interactions, i.e. FDR-adjusted p-values (q-values) were superimposed on the STRING PPIs as follows: large line width, q<0.0005; medium line width, q<0.05; small line width, q>0.05 or no q data. (B) Based on their cellular expression, the three class 2 low BMD candidate genes Rab31p, Ncald, and Arl4d were probed in the osteoclastic and osteoblastic compartments of bone turnover. Rab3ip and Ncald interacted via single intermediates with proteins degrading bone matrix and also mineral, while Arl4d, also via an intermediate, interacted with a regulator of the procollagen cleavage essential for bone matrix formation.

https://doi.org/10.1371/journal.pgen.1009190.g006

Because BMD proteins and their networks are contingent on the transcription of the corresponding BMD genes, we examined whether BMD gene expression was regulated on the transcriptional level. Within the pool of the 200 BMD genes, we probed for enrichment of TFs in promoters of genes causing low compared to high BMD phenotypes upon deletion (Table 2). A total of 49 and 32 TFs were identified for low and high BMD, respectively. In addition, we examined the 200 BMD genes for TFs and found 17 and 13 TFs that caused low and high BMD phenotypes, respectively (S6 Table). Together, these data (1) demonstrated the transcriptional regulation of BMD phenotypes, (2) recognized a total of 127 TFs involved in the regulation of low and high BMD genes, and (3) discovered 74 TFs that were not reported previously to be skeletal TFs.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. TFs enriched in low or high BMD genes (i-CisTarget analysis).

https://doi.org/10.1371/journal.pgen.1009190.t002

Ethics statement

All animal work described in this study was carried under the auspice of approved animal protocols (Baylor College of Medicine, #AN-5896; German Mouse Clinic Helmholtz Zentrum München, #144–10, 15–168; Institut Clinique de la Souris Mouse Clinical Institute, #4789-2016040511578546v2; Medical Research Council Harwell, #30/3384; Nanjing University, #NRCMM9; Rikagaku Kenkyūjo Tsukuba Institute, #Exp11-011, 12–011, 13–011, 14–009, 14–017, 15–009, 16–008; The Centre for Phenogenomics, #0153, 0275, 0277, 0279; The Jackson Laboratory, #11005, UC Davis, #20863).

Generation of mutant mouse strains

The IMPC systematically phenotypes mice that are homozygous for a single-gene knockout or heterozygous when homozygotes are lethal or sub-viable. IMPC members and non-members are free to nominate genes for deletion, and mouse production was coordinated by iMits (Fig 1A). The gene-targeting strategies that are used can be accessed via http://www.mousephenotype.org/about-ikmc/targeting-strategies. For every IMPC deleted gene, specifics on the targeting strategy can be obtained through a gene search on the IMPC website (http://www.mousephenotype.org) and use of the “Order Mouse and ES Cells” tab. In the case of single copy genes, hemizygous knockout mice are studied. IMPC mouse models are available to the research community via the IMPC website.

Mouse phenotyping

Data was derived from postnatal mice that were phenotyped under the adult and embryonic phenotype pipeline (https://www.mousephenotype.org/impress). Briefly, this pipeline conducted 16 tests, each with a set of phenotyping procedures. Phenotyping recorded appearance, behavior, or organ function across a spectrum of organs and tissues. All IMPC phenotyping data is shared with the public through the IMPC website. Experimental procedures are detailed also under IMPRESS (https://www.mousephenotype.org/impress). Although we mined the entire phenotyping data when investigating co-phenotypes, our study primarily investigated the X-ray phenotyping data obtained in postnatal week 14. The skeletal X-ray exams encompass radiographs and DXA-based measures of BA and BMC. Radiography and DXA are mandatory and voluntary IMPC measures, respectively. Specifically, for the data reported in this study, they were carried out using the UltraFocus DXA (Faxitron Bioptics LLC), the Lunar Piximus II (GE Medical Systems), or the pDexa sabre (Norland Stratec) instruments. For DXA, a minimum of 7 male and 7 female mutants were phenotyped, giving a minimum of 14 mice per line, while for radiography phenotyping was performed on a minimum of 4 male and 4 female mutants. In all experiments, mutants were matched with wild-type (wt) animals, with matching genetic background. Mutant and wt mice were phenotyped in the same manner based on the IMPReSS protocol, and detailed experiment characteristics captured in the procedure metadata. BMD was calculated from the measures of BA and BMC. For data analysis, mutants were matched with wt animals from the same center that used the same metadata.

IMPC data analysis

The IMPC uses a bespoke statistical package, implemented in R, called PhenStat, which was developed for identification of abnormal phenotypes from high-throughput pipelines [54]. The bone mineral data was analyzed by PhenStat using linear, mixed models, which take factors such as sex and body weight into account. S1 Table details the number of wt control mice analyzed for each of the 200 BMD genes. PhenStat identifies phenotype hits by sex to enable identification of sex-specific phenotypes which have been previously shown to be relevant in a large proportion of IMPC lines across many procedures [17].

Phenotype hits were identified by searching the IMPC database (version 6.0, released November 10^th, 2017, data access: https://www.mousephenotype.org/data/previous-releases/6.0) for Mammalian Phenotype (MP) terms as listed in IMPReSS. For bone mineral density, that includes the terms Increased Bone Mineral Density (MP:0000062) and Decreased Bone Mineral Density (MP:0000063). For those lines found to have BMD phenotypes, radiography-related MP terms, as listed in the IMPReSS X-ray procedure, were also queried to identify any skeletal phenotypes (www.mousephenotype.org/impress/procedureontologies/91/7).

Gene classification and assessment of sexual dimorphism

To learn more about the BMD phenotype genes and to classify genes as known or unknown to be associated with bone PubMed was searched using the gene names along with the search terms “osteoblast”, “osteoclast”, “osteocyte”, “osteogenic”, “skeletal” and “bone” NOT “bone marrow”. These broad terms were selected empirically to cover a wide range of bone biology. We also repeated the same searches, but with the inclusion of “mouse” and “human”. Our search did not encompass other skeletal tissues such as mesenchymal stem cells, cartilage, dental tissues, or malignancies. The literature search was closed on June 16, 2018.

Genes controlling BMD in a sex-biased fashion were retrieved from the 200 BMD phenotype genes. MGI and GeneCards were used to identify the chromosomal localization of genes in mice and humans, respectively. To identify genes described previously in the context of sexual dimorphism we used PubMed in combination with the search terms “sexual dimorphism”, “sexual”, “sex”, “male AND female”.

Mapping genes to gene ontology terms

GO terms for the genes of interest were extracted from the MGI database (www.informatics.jax.org), which contains information on genes and their annotated phenotypes. In order to avoid circular findings, as the MGI database includes IMPC phenotype data, IMPC entries were removed using MGI’s internal filters. We analyzed terms of the “molecular function” and “biological processes” domains as they provide information on the biological activity of genes of interest. To identify genes involved in skeletal function based on gene ontology term annotation, we selected all gene ontologies that included the terms “bone”, “skeletal”, “ossification”, “osteoclast” and “osteoblast”. MGI’s Gene eXpression Database www.informatics.jax.org/expression.shtml) was also searched to determine the known expression profiles of the BMD phenotype genes.

Genome-wide association studies (GWAS)

Orthologous/paralogous genes of mouse metabolism genes mapping to the human genome were used for analysis. GWAS analyses on human BMD were based on data acquired in the UK Biobank (www.ukbiobank.ac.uk) [9, 10, 55] and GEFOS (www.gefos.org) [8] studies. We used the Ensembl BioMart online tool to adjust all genomic data to the latest Assembly (GRCh38/hg38). The UCSC genome browser was used to analyze genes in detail. We searched for SNPs in a ±2 kb region in the SNiPA database. For each SNP occurring in or around genes, we evaluated the extent of sharing based on the three bone density phenotypes recorded by the GEFOS consortium, and used metadata without any link to individual IDs or data only. We used cross phenotype meta-analysis (CPMA), which detects association of a SNP to multiple, but not necessarily all, phenotypes [56]. The CPMA analysis applies the likelihood ratio test that measures the likelihood of the null hypothesis (i.e., that the significant SNP is uniformly distributed across consortiums) over the alternative hypothesis.

OBCD data

The OBCD (www.boneandcartilage.com) is a project performing bone and cartilage-related phenotyping on knockout mouse models. The bone-based procedures performed include three-point bend and vertebra compression, for which they currently have data for a total of 410 lines. These procedures were carried out in female animals only. In male animals, rapid-throughput joint phenotyping was performed instead. Analysis of the female OBCD bone data set was based on a tight reference data obtained from a population of >300 strain-matched female wt mice, giving it the power to detect differences between controls and small cohorts (two female animals per line) of mutant mice [5, 9, 57]. To identify lines with data differing from the wt, those with measurements greater than two standard deviations from the wt mean were considered significant.

Gene expression analysis

The NCBI Gene Expression Omnibus was used to analyze gene expression in osteoblasts and osteoclasts. To assure robustness of the expression data, the analysis was restricted to the curated data available under https://www.ncbi.nlm.nih.gov/sites/GDSbrowser/, and only data on mouse-derived cells was considered. To identify genes expressed in osteoclasts or osteoblasts, we assessed expression data acquired during differentiation of these cells. Genes significantly (p<0.05) increased during differentiation were considered. For osteoclastogenesis, data sets GSE43811 and GSE57468 (series GDS5398 and GDS5422, respectively) were used. GSE43811 compared osteoclast differentiation at 48 hrs between osteoclast precursor cells (GSM 1071626, 1071627, 1071628) and control cells (GSM 1071629, 1071630, 1071631). In Additionally, we selected GSE57468, which compared differentiation of primary osteoclast precursors between days 0 (GSM1383260, 1383262) and 3 (GSM1383261, 1383263). For osteoblastogenesis, the data set GSE2332 (series GDS1631) was selected because it permitted for the study of highly enriched primary pre-osteoblasts. Differentiation days 7 (GSM43188, 43202, 43209) and 17 (GSM43195, 43216, 43225) were included. All data was analyzed using the dedicated NCBI Geo2R suite (https://www.ncbi.nlm.nih.gov/geo/geo2r/).

Examination of protein-protein interaction (PPI) networks and the in silico bone turnover model

All PPI networks were built in STRING DB (https://string-db.org/). The proteins coded by the 200 BMD genes served as input or “seeds”. Our analysis utilized minimum networks generated by NetworkAnalyst. Based on the “String04 all sources” mode, up to 86% of the seeds were annotated (172 BMD proteins) with at least one PPI. Using these networks we calculated different centrality indexes using the cytoHubba app in Cytoscape to identify hubs central for maintaining the communication flow on the different parts of the network. To assess the significance of the network, it was then tested for the likelihood of observing a network with a the same degree of connectivity by comparing the reciprocity of the network versus 100,000 simulated networks of the same degree using the statnet and sna R packages (https://cran.r-project.org/web/packages/statnet/index.html, https://cran.r-project.org/web/packages/sna/index.html). For combined analysis of BMD proteins and TFs, the yfiles layouts were employed in Cytoscape. The radial layout algorithm was chosen. It portraits interconnected ring and star topologies. This algorithm produces layouts that emphasize sub-groups and tree structures within a network. It creates node partitions by analyzing the connectivity structure of the network and arranges the partitions as separate circles or disks. To develop networks representing an in silico bone turnover model, we used clinical markers of bone turnover as a guide for the identification of the most established proteins participating in bone turnover [58, 59]. Our final selection of bone turnover genes was supported by published data [60–62]. To experimentally validate the in silico model, we used gene expression data together with the R package GeneNet 1.2.13 (http://www.strimmerlab.org/software/genenet/). Based on gene expression data this package enabled the statistical assessment of the reliability of a given network (nodes and edges) and to compute q-values and posterior probabilities (= 1-local FDR) for each potential edge connecting in a pairwise manner the nodes of the network (S9 Table). We captured the bone resorption-formation sequence characteristic for bone turnover with the following expression data sets (1) GSE43811 (GSM 1071626, 1071627, 1071628 and GSM 1071629, 1071630, 1071631), (2) GSE57468 (GSM1383261, 1383263 and GSM1383260, 1383262), (3) GSE2332 (GSM 43191, 43212, 43219 and GSM 43184, 43205, 431198), (4) GSE37676 (GSM 25432, 25433, 25434 and GSM 25429, 25430, 25431), (5) GSE2332 (GSM43195, 43216, 43225 and GSM43188, 43202, 43209). The normalized log2FC expression values (GSE158151) from these datasets were used as input files for GeneNet [63].

Promoter analysis

To identify TFs that may regulate genes implicated in bone phenotypes, we surveyed for the bone transcriptional targets, which contain common TF binding sites in their cis-regulatory control elements. Using i-CisTarget and two datasets as input; genes with increased bone mineral density and genes with decreased bone mineral density, we detected motifs within 10Kb around the transcription start site. Every gene was scanned for motifs using a library of approx. 10k PWMs and 1120 Chip-seq tracks, each motif is scored with an algorithm called Cluster-buster and results in a subset of genes that are predicted as direct targets.