Rapid-throughput skeletal phenotyping of 100 knockout mice identifies 9 new genes that determine bone strength.

Osteoporosis is a common polygenic disease and global healthcare priority but its genetic basis remains largely unknown. We report a high-throughput multi-parameter phenotype screen to identify functionally significant skeletal phenotypes in mice generated by the Wellcome Trust Sanger Institute Mouse Genetics Project and discover novel genes that may be involved in the pathogenesis of osteoporosis. The integrated use of primary phenotype data with quantitative x-ray microradiography, micro-computed tomography, statistical approaches and biomechanical testing in 100 unselected knockout mouse strains identified nine new genetic determinants of bone mass and strength. These nine new genes include five whose deletion results in low bone mass and four whose deletion results in high bone mass. None of the nine genes have been implicated previously in skeletal disorders and detailed analysis of the biomechanical consequences of their deletion revealed a novel functional classification of bone structure and strength. The organ-specific and disease-focused strategy described in this study can be applied to any biological system or tractable polygenic disease, thus providing a general basis to define gene function in a system-specific manner. Application of the approach to diseases affecting other physiological systems will help to realize the full potential of the International Mouse Phenotyping Consortium.

Introduction

Studies of extreme phenotypes in humans have been instrumental in identifying molecular mechanisms underlying rare single gene disorders as well as common and chronic diseases including diabetes and obesity. Such studies have resulted in novel treatments that revolutionize the lives of affected individuals [1]–[4]. Collection of suitable cohorts, however, is expensive and takes many years to achieve, and progress has been limited to conditions in which simple and quantitative phenotypes can be defined [4]–[7]. By analogy, we hypothesized that an organ-specific extreme phenotype screen in knockout mice would more rapidly identify new genetic determinants of disease and also provide in vivo models to elucidate their molecular basis. The International Knockout Mouse Consortium (IKMC) has now established an ideal resource of mutant ES cells to test this hypothesis [8], [9].

The skeleton represents a paradigm organ system and osteoporosis is an important global disease ideally suited to such an approach. Osteoporosis is the commonest skeletal disorder affecting hundreds of millions of people worldwide and costing tens of billions of pounds each year [10]. Between 50 and 85% of the variance in bone mineral density (BMD) is genetically determined [11], but only 3% is accounted for by known genetic variation [12] and the vast majority of genes involved remain to be identified [13]. Current treatments reduce fracture risk by only 25–50% [14], [15] and thus there is urgent need to define new pathways that regulate bone turnover and strength in order to identify novel therapeutic targets. Accordingly, application of an extreme phenotype approach to study skeletal disorders in humans has already led to discovery of SOST (ENSG00000167941) and LRP5 (ENSG00000162337) as critical regulators of Wnt signaling in bone [5], [7], [16] and resulted in development of new drugs to stimulate bone formation [17].

The Wellcome Trust Sanger Institute Mouse Genetics Project (MGP) is undertaking high-throughput production of knockout mice using targeted ES cells generated by the IKMC. Knockout mice are generated using a knockout-first conditional gene targeting strategy (Figure S1), in which expression from the targeted allele can be investigated by X-gal staining for LacZ gene expression [9] (Figure S2). Each mouse undergoes a broad-based primary screen to identify developmental, anatomical, physiological and behavioral phenotypes [18]–[20]. A critical challenge now is to enhance this initial screening by developing organ- or disease-specific approaches [21] that are essential to identify biologically significant and functionally relevant phenotypes rapidly and cost-effectively for the benefit of the scientific community [18], [19], [21].

We, therefore, developed high-throughput skeletal phenotyping methods and prospectively studied 100 consecutive unselected mutant strains from the MGP. Using this approach, we discovered nine new genetic determinants of bone mass and strength and identified a novel functional classification of bone structure. These conditional knockout mice [9] can now be used to investigate cell-specific gene function and identify new regulatory pathways in the skeleton. The strategy can be applied to other physiological systems and complex diseases, thus realizing the full potential of the International Mouse Phenotyping Consortium.

Results

Broad primary phenotype screening

Mice generated by the MGP pipeline undergo a broad primary phenotype screen followed by terminal collection of blood and tissue at 16 weeks of age [20]. The screen is conducted on viable homozygote mutants, or heterozygotes in cases of embryonic lethality, and reports 233 variables relating to 28 physiological systems that include embyrogenesis; reproduction; growth; neurological; behaviour; sensory; skeleton; muscle; gastrointestinal and hepatobiliary; cardiovascular; endocrine; adipose; metabolism; haematopoietic; immune; skin and pigmentation; respiratory; and renal. Parameters relevant to the skeleton include body length, x-ray skeletal survey, dual energy x-ray absorptiometry (DEXA) analysis of BMD and biochemical measures of mineral metabolism. Tissues in which the targeted gene is expressed are determined by staining for lacZ reporter gene expression in heterozygous mice (Figure S2). To extend this broad initial screen, we incorporated novel imaging, statistical and biomechanical approaches for the specific and sensitive detection of functionally important skeletal abnormalities (Figure S3).

Generation of normal reference data

In order to establish strain-specific reference ranges for these new approaches, limbs from 16 week-old female C57BL/6 (B6Brd;B6Dnk;B6N-Tyr) wild-type mice (n = 77) were obtained from 18 control cohorts. Normal ranges for six independent parameters of bone structure were obtained using Faxitron x-ray point projection microradiography and micro-computed tomography (micro-CT) (Fig. 1 and Figure S4). Bone mineral content (BMC), bone length and cortical bone thickness were determined by x-ray microradiography and measures of trabecular bone volume per tissue volume (BV/TV), trabecular number (Tb.N) and trabecular thickness (Tb.Th) by micro-CT (Fig. 1). Reference data were also obtained for six biomechanical parameters (Fig. 2 and Figure S5). The yield, maximum and fracture loads, stiffness and the proportions of energy dissipated prior to maximum load and fracture were determined from load displacement curves obtained in destructive 3-point bend tests (Fig. 2).

Figure 1

Faxitron x-ray microradiography and micro-CT.

A, Faxitron femur images from WT and Sparc mice (arrows indicate location for cortical bone thickness measurement, bar = 1 mm). B, Bone mineral content (BMC) in WT and Sparc mice. Pseudo-colored images in which lower BMC is in green and yellow and higher BMC is red and purple. C, Cumulative frequency histograms of BMC in n = 77 female, 16 week-old WT (mean ±2.0SD reference range in grey) and Sparc mice (red line). The median grey level is indicated by the dotted line. Graphs showing mean (solid line), 1.0SD (dotted line) and 2.0SD (grey box) for D, median grey level BMC, E, bone length and F, cortical thickness in WT (n = 77) mice. Values for Sparc (n = 2) in red. Micro-CT tibia images from G, WT and H, Sparc mice (bar = 1 mm). Graphs showing mean, 1.0SD and 2.0SD for I, BV/TV, J, trabecular number (Tb.N) and K, trabecular thickness (Tb.Th) in WT mice. Values for Sparc in red.

Figure 2

Biomechanical analysis.

A, Load-displacement curve from a WT tibia showing yield load, maximum load, fracture load and gradient of the linear elastic phase (stiffness). B, Curves from WT and Sparc. C, Graphs showing mean (solid line), 1.0SD (dotted line) and 2.0SD (grey box) for yield load, maximum load, fracture load and stiffness of WT (n = 77) mice. Values for Sparc in red. D, Energy dissipated prior to maximum load (DEML, purple) and elastic stored energy at maximum load (ESEML, yellow). E, Graph showing mean ±1.0SD and 2.0SD for the proportion DEML/(DEML+ESEML) prior to maximum load for WT mice. Value for Sparc in red. F, Energy dissipated prior to fracture (DEF, purple) and elastic stored energy at fracture (ESEF, yellow). G, Graph showing mean ±1.0SD and 2.0SD for the proportion DEF/(DEF+ESEF) prior to fracture for WT mice. Value for Sparc in red. y-axis scale reflects angular transformation to normalize data distribution.

Validation of phenotyping methods

Limbs from 16 week-old female knockout mice in an identical C57BL/6 genetic background were obtained prospectively from the MGP pipeline (n = 100 unselected knockout strains, 2–6 mice per strain) and analyzed by x-ray microradiography, micro-CT and 3-point bend testing. Serendipitously, one of the 100 unselected strains was a homozygous knockout of Sparc (ENSMUSG00000018593), which encodes the extracellular matrix glycoprotein osteonectin. Deletion of Sparc is known to cause low bone turnover osteopenia resulting in weak and brittle bone with a higher mineral-to-matrix ratio due to reduced bone matrix content [22], [23]. Thus, Sparc knockout mice represented a well-characterized positive control for validation of our approach. Consistent with the reported phenotype [22], [23], we identified that Sparc knockout mice had reduced BMD and BMC with loss of trabecular bone but preservation of cortical bone thickness (Fig. 1 and Table S1), resulting in weak and brittle bone of reduced stiffness (Fig. 2). We also identified short stature in Sparc knockout mice (Fig. 1A), a parameter not investigated in previous studies. These findings validate the use of complementary and multi-parameter imaging together with biomechanical methods as a rapid and specific phenotyping approach to identify biologically significant and functionally relevant skeletal abnormalities using a minimal number of animals (n = 2).

Phenotyping of 100 knockout strains

To identify new genetic determinants of bone mass and strength, limbs from 100 knockout strains were analyzed for each structural and biomechanical variable. X-ray microradiography and micro-CT imaging identified 19 knockout strains in which at least one structural parameter was >2.0 standard deviations (SD) from the reference mean (Fig. 3, Figure S4 and Table S1). To ensure that significant abnormal phenotypes resulting from simultaneous but smaller variances in any of the six parameters were not overlooked, Mahalanobis distances were calculated as detailed in the methods and principal component analysis performed to identify multivariate outliers [24]–[27]. These studies identified 40 strains with outlier Mahalanobis distances (P<0.025), 24 of which had not been identified by analysis of individual x-ray microradiography or micro-CT values alone (Fig. 3 and Table S1). The MGP broad primary phenotype screen independently annotated 17 of these knockout strains with skeletal abnormalities (Table S1). Nine of the strains were also identified as outliers by x-ray microradiography, micro-CT or statistical methods whereas 8 did not display any abnormalities (Fig. 3).

Figure 3

Knockout strains with abnormal skeletal phenotypes.

Venn diagram showing strains with at least one outlier structural parameter >2.0SD from the C57BL/6 reference mean determined by Faxitron, micro-CT, Mahalanobis distance calculation or primary phenotype screening. Strains with at least one outlier biomechanical parameter in blue. 10 strains with major phenotypes are highlighted in boxes.

The biomechanical significance of the 43 outlier phenotypes identified by imaging (n = 19, Faxitron and micro-CT) and statistical (n = 24, Mahalanobis analysis but excluding Faxitron and micro-CT) approaches, together with the 8 additional strains identified only in the MGP primary screen, was investigated (Fig. 3). Twelve of the 51 strains had at least one biomechanical parameter >2.0 SD from the reference mean (Fig. 3, Figure S5 and Table S1). However, 2 of these 12 strains had only minor abnormalities of bone morphology in the primary screen and were normal when investigated by x-ray microradiography, micro-CT and principal component analysis. Destructive 3-point bend testing of bones from the remaining 49 strains identified a further 7 with a single outlier biomechanical parameter but no other abnormality (Fig. 3 and Table S1).

In summary, the broad primary phenotype screen together with x-ray microradiography, micro-CT and statistical analysis identified knockout strains with at least one abnormal bone-related parameter. The addition of functional biomechanical testing demonstrated that 10 of these strains had major phenotypes affecting both the structure and strength of bone. Three of these carried heterozygous mutations (Asxl1 (ENSMUSG00000042548), Setdb1 (ENSMUSG00000015697) and Trim45 (ENSMUSG00000033233)) while the rest were homozygotes (Bbx (ENSMUSG00000022641), Cadm1 (ENSMUSG00000032076), Fam73b (ENSMUSG00000026858), Prpsap2 (ENSMUSG00000020528), Slc38a10 (ENSMUSG00000061306), Sparc and Spns2 (ENSMUSG00000040447)).

The primary phenotype database (http://www.sanger.ac.uk/mouseportal/) was interrogated for each of the strains identified with major skeletal phenotypes (Table 1).

Table 1

Functions of genes identified as determinants of bone strength.

Gene	Expression	Protein	Biological Function	Skeletal Function	Primary Phenotyping Features	Ref
Weak but Flexible (Low BMC)
Bbx(NM_027444)	Brain, spinal cord, peripheral nervous system, testis and cartilage	Bobby sox homolog (HMG box transcription factor)	Unknown, but related protein family members include SRY-related Sox genes and TCF-LEF mediators of Wnt signalling	Sox genes and Wnt signalling regulate skeletal development and bone formation, mineralization and turnover. Bbx may have a related function	Homozygotes viable and fertile.Decreased lean body mass and cardiac weight, increased plasma IgA.Skeletal abnormalities include decreased nose to tail length, dental asymmetry and low BMD.	[28]–[35]
Cadm1(NM_001025600)	Ambiguous pattern in stomach, testis, nasal epithelia and bone, but absent from cartilage	Cell adhesion molecule 1 (immunoglobulin superfamily cell adhesion molecule)	Tumour suppressor gene with a role in retinoid-regulated synapse development	Unknown	Homozygous femalesfertile, but males infertile.No additional abnormalities.	[36]–[38]
Fam73b(NM_001242407)	Most tissues including cartilage but not bone	Family with sequence similarity 73 member B (membrane protein with conserved 500–600 amino acid domain)	Unknown	Unknown	Homozygotes viable and fertile.Decreased body weight and lean body mass and increased susceptibility to infection following citrobacter challenge, decreased serum total protein and albumin.Skeletal abnormalities include decreased body length, low BMD in females and abnormal tooth morphology in males.
Weak and Brittle (Low BMC)
Prpsap2(NM_001164244)	Most tissues including cartilage but not bone	Phosphoribosyl pyrophosphate ynthetase-associated protein 2 (non-catalytic subunit of phosphoribosylpyrophosphate synthetase)	Catalyzes formation of phosphoribosylpyrophosphate substrate for synthesis of purine and pyrimidine nucleotides	Candidate oncogene in osteosarcoma	Homozygotes viable and fertile.Increased rearing in males and decreased plasma IgG1.	[39], [40]
Slc38a10(NM_001164798)	Most tissues including bone and cartilage	Solute carrier family 38 member 10 (sodium-coupled neutral amino acid transporter)	Proposed as a cell volume regulator in mesenchyme	Growth defect suggests role in hypertrophic chondrocytes, which mediate growth by cell volume expansion	Homozygotes viable and fertile.Decreased fat and leanmass, increased oxygen consumption and energy expenditure, reduced serum amylase in females, increased creatinine and low albumin in males.Skeletal abnormalities include reduced body length and decreased BMD.	[41]–[43]
Sparc(NM_009242)	Ubiquitous	Secreted protein acidic and rich in cysteine (osteonectin - extracellular matrix glycoprotein)	Cell migration and tissue remodelling during development and in response to injury	Regulation of collagen assembly during bone formation and turnover, regulation of osteoclast maturation and function	Homozygotes viable and fertile.Cataracts.Skeletal abnormalities include abnormal teeth and decreased BMD.	[22], [23]
Strong but Brittle (High BMC)
Asxl1(NM_001039939)	Most tissues including bone and cartilage	Additional sex combs-like 1 (conserved polycomb chromatin-binding protein)	Cooperates with heterochromatin protein-1 to repress retinoic acid signaling. Mutations cause Bohring-Opitz syndrome	Regulates Hox genes during axial patterning, suggesting role in skeletal development	Homozygous lethal due to craniofacial defects.Heterozygotes viable and fertile.Minor defects in lumbar and sacral vertebrae in heterozygotes.	[44]–[48]
Setdb1(NM_001163641)	Restricted pattern including cartilage but not bone	SET (Su(var)3-9, Enhancer-of-zeste, Trithorax) domain bifurcated-1 (histone-lysine N-methyltransferase)	Regulates gene silencing	Expressed in osteoblasts with possible role in lineage commitment and differentiation	Homozygous embryonic lethal.Heterozygotes viable and fertile.Increased natural killer lymphocytes and CD4+ T cells in female heterozygotes.	[49]–[51]
Spns2(NM_153060)	Restricted pattern including cartilage and bone	Spinster homolog 2 (spingosine-1-phosphate transporter)	Required for secretion of sphingosine 1-phosphate (S1P), which binds to the G-protein receptors, S1PR1 and S1PR2	Modulates osteoclast and osteoblast precursor cell recruitment and migration. Regulation of S1P secretion may represent new mechanism to control coupling of bone formation to resorption	Homozygotes viable and fertile.Abnormal vision and hearing, reduced blood glucose and increased bilirubin in males, abnormal lymphocyte, granulocyte and monocyte fractions in females.	[52]–[56]
Trim45(NM_001165952)	Brain and testis	Tripartite motif containing protein 45 (putative ubiquitin or SUMO E3 ligase)	Interacts with AP-1 and inhibits MAP kinase activity	Role in the skeleton unknown, although AP-1 regulates osteoblast differentiation and bone formation.	Homozygous lethal due to exencephaly. Heterozygotes viable and fertile.	[57]–[61]

Table summarizing the known and proposed functions of the 10 genes which, when deleted, result in major skeletal phenotypes affecting bone structure and strength. Patterns of gene expression and abnormalities identified by primary phenotype screening are also included.

Accuracy of rapid throughput skeletal phenotyping

The sensitivities, specificities and predictive values of each phenotyping method at a statistical threshold of >2.0 SD and 95% confidence limit were calculated to determine their ability to identify the 10 strains with major phenotypes. X-ray microradiography was the most accurate, identifying 8 of the 10 strains. Additional use of micro-CT and Mahalanobis analysis was required to identify the remaining 2 strains (Table S2). The MGP primary phenotype screen identified 5 of the 10 abnormal strains. Thus, addition of organ-specific imaging, statistical and biomechanical analyses to the primary phenotype screen resulted in increased sensitivity and specificity (Table S2), demonstrating the advantage of a complementary multi-parameter approach.

Correlations between imaging and biomechanical parameters were also determined to investigate relationships between bone structure and strength (Figure S6). Bone strength correlated strongly with BMC, cortical thickness and bone length but not with BV/TV, Tb.Th or Tb.N. Furthermore, there was no significant relationship between bone structural parameters determined by x-ray microradiography and micro-CT, thus demonstrating independence of the two techniques. Consistent with sensitivity, specificity and predictive value data (Table S2) these findings demonstrate that mineralization and cortical bone parameters determined by x-ray microradiography are excellent predictors of bone strength determined by 3-point bend testing, whereas trabecular bone structure determined by micro-CT is not. Importantly, trabecular bone parameters determined by micro-CT may represent good predictors of bone strength at sites of predominantly cancellous bone such as the vertebra or in response to age-related bone loss, although these possibilities were not investigated.

Functional classification of bone structure

Biomechanical analysis of wild-type mice and the 10 knockout strains with major phenotypes enabled four distinct biomechanical categories to be identified (Fig. 4). Normal bone had a stiffness of 30.2±4.1 N/mm (mean ± SD) with the capability to resist loading up to a yield load of 8.9±0.9 N (mean ± SD) and maximum load of 10.7±0.9 N (mean ± SD), prior to fracture at a load of 6.4±1.7 N (mean ± SD). These material properties of normal bone represent an optimised compromise between strength and flexibility that allows dissipation of 86.2% of energy prior to failure and limits the structural damage at fracture (Table S1).

Figure 4

Functional classification of bone structure.

A, Load-displacement curve from WT tibia showing 2.0SD distribution of C57BL/6 reference range in grey. B, Curves from Bbx, Cadm1 and Fam73b mice with weak but flexible bones and low bone mineral content (BMC). C, Curves from Sparc, Prpsap2, and Slc38a10 mice with weak and brittle bones and low BMC. D, Curves from Asxl1, Trim45, Spns2 and Setdb1 mice with strong but brittle bones and high BMC. E, Proportion of energy dissipated prior to fracture (DEF/DEF+ESEF) versus maximum load. The y-axis scale reflects angular transformation to normalise data distribution. Strains with major phenotypes in red and individual WT mice in black. The plot separates four functional categories of bone structure that include normal bone which is strong and flexible with normal BMC and the three abnormal categories in B, C, and D.

Bones from strains with major phenotypes clustered into three abnormal biomechanical categories. Bones from Bbx, Cadm1 and Fam73b knockout mice were weak but flexible with reduced maximum load but the capability to bend and dissipate energy prior to fracture (Fig. 4B). Sparc, Prpsap2 and Slc38a10 bones were weak and brittle with reduced maximum load and lacking the capability to bend and dissipate energy prior to fracture (Fig. 4C). Asxl1, Trim45, Spns2 and Setdb1 bones were strong but brittle with an increased maximum load but were unable to bend and dissipate energy prior to fracture (Fig. 4D). A plot of the proportion of energy dissipated prior to fracture versus maximum load clearly separates these categories of bone strength (Fig. 4E). Further analysis demonstrated that Bbx, Cadm1, Fam73b Sparc, Prpsap2 and Slc38a10 bones have low BMC, whereas Asxl1, Trim45, Spns2 and Setdb1 bones have high BMC (Fig. 5). To determine whether these three categories of abnormal bone strength were related to a further morphological parameter known to have an important role for the biomechanical properties of long bones, we investigated their relationship with mid-diaphyseal diameter (Fig. 6). As expected [28], in WT mice mid-diaphyseal diameter correlated with fracture load and the proportion of energy dissipated prior to fracture. However, bones from the mutant strains identified with major skeletal phenotypes clustered into the same three abnormal categories, further demonstrating the validity of this functional classification.

Figure 5

Knockout strains with major phenotypes affecting bone structure and strength.

A, Digital radiographs of femurs from WT mice and each of the 10 knockout strains with major phenotypes (bar = 1 mm). B, Magnified images of mid-diaphysis, the region where cortical thickness was determined (bar = 1 mm). C, Grey-scale images pseudo-coloured using a 16-colour palette in which lower BMC is in green and yellow and higher BMC is red and purple (bar = 1 mm). D, Cumulative frequency histograms of whole femur BMC in WT mice and knockout strains: Bbx, Cadm1 and Fam73b mice with weak but flexible bones and low BMC (left); Prpsap2, Slc38a10 and Sparc mice with weak and brittle bones and low BMC (middle); and Asxl1, Setdb1, Spns2 and Trim45 mice with strong but brittle bones and high BMC (right). E, Transverse sections of tibias from WT and knockout mice imaged by micro-CT (bar = 1 mm). F, Mid-sagittal longitudinal sections of tibias from WT and knockout mice imaged by micro-CT (bar = 1 mm).

Figure 6

Relationship between mid-diaphyseal cortical bone diameter and strength.

A, Graphs showing mid-diaphyseal cortical bone diameter mean (solid line), 1.0SD (dotted line) and 2.0SD (grey box) in mutant strains with weak but flexible, weak and brittle, and strong but brittle bones. B, Relationship between fracture load and mid-diaphyseal cortical bone diameter. Strains with major phenotypes in red and individual WT mice in black. The 2.0SD reference range for each variable is represented by the grey box. The plot separates four functional categories of bone structure that include normal bone which is strong and flexible with normal BMC and the three abnormal categories weak but flexible (low BMC, green), weak and brittle (low BMC, purple) and strong but brittle (high BMC, orange). C, Relationship between energy dissipated prior to fracture (DEF/(DEF+ESEF)) and cortical bone diameter. The y-axis scale reflects angular transformation to normalize data distribution. The same functional categories of bone structure are separated by this plot.

Investigation of the biological activities of the 10 genes and their possible roles in bone indicates a broad diversity of function that was not clearly related to phenotype (Table 1), thus reinforcing the importance of an unbiased screening approach for identification of gene function in both homozygous and heterozygous mutants. In summary, a functional classification of four categories of bone structure was defined. Normal bone is strong and flexible with a normal mineral content, whereas abnormal bone is either (i) weak but flexible with low BMC, (ii) weak and brittle with low BMC or (iii) strong but brittle with high BMC.

Discussion

We have identified 10 genes with diverse and unrelated functions, the deletion of which resulted in major skeletal abnormalities. By adopting a multi-parameter phenotyping approach, we identified a new functional classification of bone structure based on its mineral content, strength and ductility that clarifies understanding of skeletal physiology and pathology, and which maps directly to human disease. As a result of evolutionary pressure, bone structure represents an optimal compromise between strength and flexibility that requires contributions from many diverse genes. Continuous bone remodeling enables the skeleton to adjust this compromise in response to changing physiological and environmental pressures [29], [30]. The current studies demonstrate that loss of function of individual genes can disrupt this optimal compromise resulting in skeletal phenotypes that cluster into three functionally distinct categories.

Postmenopausal osteoporosis is characterized by weak but flexible bone with low mineral content [31] and three of the identified knockout strains (Bbx, Cadm1, Fam73B) had phenotypes in this category. Bbx encodes a conserved transcription factor that contains a SOX-TCF HMG-box [32]–[36]. Family members include SRY-related Sox genes that are implicated in skeletal dysplasias [37], [38] and the TCF/LEF transcription factors that mediate Wnt/β-catenin signaling [39], a key pathway implicated in osteoporosis and osteoarthritis [40]–[43]. Cadm1 encodes a trans-membrane glycoprotein adhesion molecule of the immunoglobulin superfamily [44] for which a number of disparate functions have been reported including; tumor suppression [45], synapse development [46], behavioral regulation [47], T cell adhesion [48], mast cell interactions [49], and spermatogenesis [50]. However, no function in the skeleton has been reported. Fam73B encodes a conserved membrane protein of unknown function. These findings indicate that deletion of genes encoding proteins with diverse and unrelated functions can result in similar defects of bone strength and mineralization.

Disorders of bone matrix as typified by osteogenesis imperfecta [51] are characterized by bone that is weak and brittle with low BMC, and three of the strains (Prpsap2, Slc38a10, Sparc) displayed this phenotype. Prpsap2 encodes the non-catalytic inhibitory subunit of phosphoribosylpyrophosphate synthetase [52], and is required for synthesis of purine and pyrimidine nucleotides, the amino acids histidine and tryptophan, and the coenzyme nicotinamide adenine dinucleotide [53]. Its function in the skeleton is unknown, although a recent study proposed PRPSAP2 as a candidate oncogene in osteosarcoma tumorigenesis [54]. Slc38a10 encodes a proposed sodium-coupled neutral amino acid membrane transporter [55] that may act as a cell volume regulator in mesenchyme [56]. The severe growth defect in Slc38a10 knockout mice suggests a critical function in chondrocytes, which mediate linear growth by cell volume expansion during hypertrophic differentiation [57]. Furthermore, related transporters have already been implicated in human skeletal disease. SLC35D1 is critical for chondroitin sulphate synthesis and mutations cause Schneckenbecken skeletal dysplasia [58]. Mutations in SLC26A2, cause four distinct chondrodysplasia syndromes [59] and emphasize the key role of these transporters in endochondral ossification. Sparc encodes the well-described extracellular matrix glycoprotein osteonectin and its deletion resulted in the characteristic and expected phenotype [22], [23] of weak and brittle bone with low BMC. These findings highlight the importance of enzymes, transporters and structural proteins to the functional integrity of bone matrix.

Diseases of high bone mass are rare and include sclerosteosis due to deletion of SOST [60] and autosomal dominant high bone mass due to gain-of-function mutations in LRP5 [61]. They are characterized by bone that is strong but brittle with high BMC, and four of the knockout strains (Asxl1, Setdb1, Spns2, Trim45) displayed such a phenotype. Asxl1 encodes a polycomb protein that interacts with heterochromatin protein-1 [62] and is required for regulation of Hox genes during axial patterning [63], suggesting a role in skeletal development [64], [65]. Indeed, ASXL1 heterozygous nonsense mutations were recently described to cause Bohring-Opitz syndrome [66], a developmental disorder characterized by mental retardation, impaired intrauterine growth, trigonocephaly and wrist and metacarpophalangeal joint abnormalities. Although the disease mechanism is unknown, craniofacial defects identified in homozygous Asxl1 knockout mice suggest that mutations in Bohring-Opitz syndrome result in a mutant protein with dominant-negative activity. Setdb1 encodes a histone H3 methyltransferase that regulates gene silencing [67], [68]. Although found to be expressed in cartilage but not bone in the primary phenotype screen, other studies demonstrated Setdb1 expression in osteoblasts and suggested a role in lineage commitment and differentiation [69], [70]. Spns2 encodes a sphingosine 1-phosphate (S1P) transporter [71] that is essential for S1P secretion. S1P binds to the G-protein coupled receptors, S1PR1 and S1PR2, and regulates osteoclast [72], [73] and osteoblast [74] precursor cell recruitment and migration. Thus, control of S1P secretion by Spns2 represents a novel mechanism that couples bone resorption and formation [75]. Trim45 is a member of the tripartite protein family, many of which act as ubiquitin or SUMO E3 ligases [76]–[78]. Although restricted to brain and testis in the primary phenotype screen, human studies demonstrate that Trim45 is more widely expressed [79]. Little is known about its function, although one study indicates Trim45 interacts with AP-1 and inhibits activity of the MAP kinase pathway [79]. The physiological significance of these findings and the role of Trim45 in the skeleton are unknown, although AP-1 proteins are key regulators of osteoblast and osteoclast differentiation and function [80], [81]. These findings emphasize the importance of lineage commitment, control of cell differentiation and coupling of both osteoblasts and osteoclasts in high bone mass disorders.

In the context of osteoporosis, our identification of many new genes that determine bone strength, and which otherwise could not be predicted, is consistent with studies indicating that diverse genetic polymorphisms result in small effects on phenotype [11], [12], [82], [83]. Accordingly, and in line with current understanding that only 3% of the heritability of BMD is accounted for by known genetic variation [12], none of the genes identified in this study have been recognized in osteoporosis genome-wide association studies [84]. We hypothesize, therefore, that unbiased multi-parameter and functional phenotyping of knockout mice has the power to identify many of the major genes that determine bone strength. Ultimately, this approach is likely to identify several genes from a single signaling pathway with an important role in the control of bone mass and strength. This has the advantage of independently confirming critical pathways and the potential to identify several alternative therapeutic targets. Importantly, however, the approach has limitations. The study of knockout mice can only identify phenotypes that result from gene deletion but cannot identify genes that only cause abnormalities when they harbor gain-of-function or dominant-negative mutations. Furthermore, the strategy does not include challenges such as ageing that may reveal additional phenotypes. However, if such provocative challenges were to be incorporated into screening approaches they would inevitably increase costs and limit throughput.

Our findings resulted from development and refinement of a rapid-throughput phenotyping algorithm to identify knockout mice with major abnormalities of bone structure and strength (Figure S3). The methods require bones from only two knockout mice, which first undergo digital point projection x-ray microradiography and micro-CT determination of six parameters of bone structure. Mahalanobis distance calculations and principal component analysis is performed and strains with at least one structural parameter >2.0 SD from the reference mean plus those with outlier Mahalanobis distances (95% confidence limit) are selected for biomechanical studies. Bones from selected strains undergo destruction 3-point bend testing to determine six measures of bone strength. Application of this unbiased approach to 100 consecutive knockout strains from the MGP pipeline identified 10% with major phenotypes affecting bone strength. Subsequent consideration of the results of primary phenotype screening and biological plausibility (Figure S3) allowed selection of mice to be refined.

Inherent in this approach is the capability to alter the statistical stringency threshold of analyses such that the number of strains for subsequent functional studies can be adjusted according to phenotype severity. For example, if the threshold for structural parameters is increased from 2.0 to 3.0 SD, then 9 outlier strains (rather than 19) are identified. Furthermore, if the confidence limit for Mahalanobis distance is increased from 95 to 99.7% then 21 multivariate outliers (rather than 40) would be identified. Of note, Trim45, which was recognized as an outlier only by Mahalanobis analysis, would still be identified if the confidence limit were to be increased to 99.7%, thus emphasizing the importance of a robust statistical method to ensure that all functional outliers are captured. Biomechanical analysis following application of these more stringent thresholds would detect 8 outlier strains including Trim45 and Sparc, resulting in the identification of 7 novel determinants of bone mass and strength rather than 9. The intrinsic flexibility of such a bespoke approach facilitates its application to other biological systems or polygenic diseases.

Materials and Methods

Ethics statement

All mouse studies were undertaken by Wellcome Trust Sanger Institute Mouse Genetics Project as part of the International Knockout Mouse Consortium and licensed by the UK Home Office in accordance with the Animals (Scientific Procedures) Act 1986 and the recommendations of the Weatherall report.

Primary phenotype screen

All mice generated by the MGP undergo a broad primary phenotype screen (http://www.sanger.ac.uk/mouseportal/) that includes measurement of body length, x-ray skeletal survey, DEXA analysis of bone mineral density and biochemical measures of mineral metabolism performed between 14–16 weeks of age, and determination of the normal tissue expression pattern of the targeted gene in 6–12 week old mice. Following primary phenotyping, lower limbs were fixed in 70% ethanol.

LacZ reporter gene expression

The pattern of LacZ reporter gene expression was determined in whole mount tissue preparations from heterozygous knockout mice between 6 and 12 weeks of age. Under terminal anaesthesia, mice were perfused with fresh cold 4% paraformaldehyde (PFA). Tissues were fixed for a further 30 min in 4% PFA, rinsed in phosphate buffered saline and stained with 0.1% X-gal for 48 hours at 4°C. Samples were subsequently fixed overnight in 4% PFA at 4°C, cleared with 50% glycerol and transferred to 70% glycerol. Specific and non-specific staining was determined in 41 tissues (Figure S2). A panel of 27 standardized images were recorded if expression was widespread (www.sanger.ac.uk/mouseportal).

Faxitron point projection digital x-ray microradiography

Bones from 16 week-old mice were fixed in 70% ethanol. Soft tissue was removed from the fixed bones and digital X-ray images were recorded at 10 µm pixel resolution using a Faxitron MX20 variable kV point projection x-ray source and digital image system (Qados, Cross Technologies plc, Sandhurst, Berkshire, UK) operating at 26 kV and 5× magnification [85]. Magnifications were calibrated by imaging a digital micrometer. Bone mineral content, bone length and cortical bone thickness were determined with coefficients of variation (CV) of 1.7%, 2.0% and 5.1%, respectively. The relative mineral content of calcified tissues was determined by comparison with standards included in each image frame, which comprised: a 1 mm steel plate; a 1 mm diameter spectrographically pure aluminum wire; and a 1 mm diameter polyester fiber. 2368×2340 16 bit DICOM images were converted to 8 bit Tiff images using ImageJ and the histogram stretched from the polyester (grey level 0) to steel (grey level 255) standards. Increasing gradations of mineralization density were represented in 16 equal intervals by a pseudocolor scheme. Cortical bone thickness was determined in at least 10 locations at the mid-femoral diaphysis. Bone length was determined using ImageJ 1.41 software (http://rsb.info.nih.gov/ij/).

Micro-CT

Tibias were analyzed by micro-CT (Skyscan 1172a, Skyscan, Belgium) at 50 kV and 200 µA using a 0.5 mm aluminum filter and a detection pixel size of 4.3 µm2 [85]. Images were captured every 0.7°, with 2× averaging, through 180° rotation of each bone and reconstructed using Skyscan NRecon software. A volume of 1 mm3 of trabecular bone was selected as the region of interest, 0.2 mm from the growth plate. Trabecular bone volume as proportion of tissue volume (BV/TV, %, CV 18.4%), trabecular thickness (Tb.Th, mm, CV 11.1%) and trabecular number (Tb. N, mm−1, CV 17.4%) were analyzed [86] using Skyscan CT analysis software.

Destructive 3-point bend testing

Bones were stored and tested in 70% ethanol. Destructive 3-point bend tests were performed on an Instron 5543 materials testing load frame (Instron Limited, High Wycombe, Buckinghamshire, UK) using custom built mounts with rounded supports that minimize cutting and shear loads [85]. Bones were positioned horizontally and centered on custom supports with the anterior surface upward. Load was applied vertically to the mid-shaft with a constant rate of displacement of 0.03 mm/second until fracture. A span of 12 mm was used. Load-displacement curves were plotted and yield load, maximum load and fracture load determined. Stiffness, the slope of the linear (elastic) part of the load-displacement curve, was calculated by the “least squares” method. Work energy was calculated from the area under curve at both maximum load and fracture. Elastic stored energy at maximum load was determined by calculating the area of a right angled triangle with the vertex at the point of maximum load and hypotenuse with a slope equal to that of the linear phase of the load-displacement curve. Elastic stored energy at fracture was similarly calculated but with the vertex of the triangle at the point of fracture (Fig. 2). Energy dissipated at maximum load or fracture was calculated by subtracting the elastic stored energy from the work energy at maximum load or fracture. CVs for each parameter were as follows: yield load (9.8%), maximum load (8.5%), fracture load (26.6%), stiffness (13.6%), the ratio of energy dissipated at maximum load to elastic stored energy at maximum load (25.1%), and the ratio of energy dissipated prior to fracture to elastic stored energy at fracture (11.0%).

Calculation of Mahalanobis distances

C57BL/6 reference ranges were generated for all Faxitron and micro-CT measures. Outliers in multivariate data were identified using robust Mahalanobis distances [27], which measure how far each observation is from the center of a data cluster, taking into account the shape of the cluster [87]. Robust Mahalanobis distances

were calculated for the vector of multivariate observations as described [27]. Here is a robust (i.e. relatively unaffected by outliers) estimate of the mean vector and is a robust estimate of the covariance matrix of the data set . Under the assumption of multivariate normality, if mouse is from the same population as the rest of the data then has a chi-squared distribution with p degrees of freedom (where p is the number of variables). Robust estimates of the mean and covariance matrix are used so that potential outliers are not masked. The masking effect, by which outliers do not necessarily have a large Mahalanobis distance, can be caused by a small cluster of outliers that attract the mean and inflate the covariance in its direction. By replacing the sample mean and covariance with a robust estimate, the influence of these outliers is removed and the Mahalanobis distance is able to expose all outliers. Robust estimates of the mean and covariance matrix were calculated using the minimum volume ellipsoid method [87]. Given observations and variables, the minimum volume ellipsoid method seeks an ellipsoid containing

points of minimum volume. All analysis was conducted in the statistical computing package R (http://www.R-project.org).

Principal Component Analysis

Principal component analysis was used as a method to visualize multivariate data and reveal outliers by describing variation in a set of correlated variables in terms of a new set of uncorrelated variables. These new variables or principal components are linear combinations of the original variables derived in decreasing order of importance so that the first component accounts for the most variation of all possible linear combinations. The second component is then selected so that it accounts for as much of the remaining variance as possible (subject to it being uncorrelated with the first component), and so on [24]. Since the first few principal components often contain most of the variation in the data set they can be used as a lower-dimensional summary of the original variables.

Statistics

Normally distributed data were analyzed by Student's t test, or ANOVA followed by Tukey's multiple comparison post-hoc test. Relationships between bone structure and biomechanical measures, and between individual bone structure parameters, were determined by Pearson correlation. P values<0.05 were considered significant. Frequency distributions of bone mineral densities obtained by Faxitron were compared using the Kolmogorov-Smirnov test, in which P values for the D statistic in 1024 pixel data sets were D = >6.01 P<0.05, D = >7.20 P<0.01, and D = >8.62 P<0.001.

International knockout mouse consortium gene targeting strategy. A, Global gene inactivation is achieved using a gene-trap targeting cassette containing a LacZ gene expressing β-galactosidase and a neomycin marker gene flanked by Flp recombinase target sites (FRT), together with a critical exon flanked by loxP sites (exon 3 in this example). Gene expression from the targeted allele can be determined by X-gal staining for β-galactosidase activity. B, Crossing targeted mice with a Flp deleter strain removes the lacZ and neomycin genes, resulting in reversal of the gene-trap knockout and generation of a floxed allele. C, Mice harboring the floxed allele can be crossed with appropriate Cre recombinase expressing strains to generate tissue-specific knockout mice for further study.

(TIF)

Click here for additional data file.

Determination of the tissue distribution of gene expression. A, Widespread tissue distribution of Asxl1 expression as determined by X-gal staining of tissue preparations for lacZ reporter gene expression in heterozygous mice. Table shows the annotations as reported on the Wellcome Trust Sanger Institute Mouse Genetics Project portal (http://www.sanger.ac.uk/mouseportal/). B, Restricted tissue expression pattern of Bbx. C, Expression of Setdb1 in rib cartilage (arrow), as demonstrated by X-gal staining (blue) of tissue preparation from a 7 week old heterozygote. D, Expression of Slc38a10 in tracheal cartilage (arrow) in an 8 week old heterozygote. E, Expression of Spns2in distal femur and proximal tibia (arrows) in an 8 week old heterozygote. Corresponding images from wild type (WT) control mice show background X-gal staining.

(TIF)

Click here for additional data file.

Multi-parameter and functional skeletal phenotyping algorithm. Summary of an unbiased and high-throughput phenotype screen to identify knockout mice with skeletal abnormalities and identify new genetic determinants of bone mass and strength.

(TIF)

Click here for additional data file.

Bone structure determined by Faxitron x-ray microradiography and micro-CT in 100 unselected knockout mouse strains. The mean value for each parameter obtained from n = 77 female 16 week old WT mice is shown as a horizontal line with the 2.0 SD reference range limits in grey. Black dots represent values for individual knockout mice. Grey vertical boxes highlight the distributions of values from mice of a single strain in which at least one animal lies outside the reference range, but the mean value lies within the reference range. Violet vertical boxes and red squares highlight the distributions of individual values from outlier strains in which the mean value (horizontal line in violet box) lies outside the reference range. Gene symbols for outlier strains are indicated. A–C, Faxitron x-ray microradiography measures of bone mineral content, bone length and cortical thickness.

(TIF)

Click here for additional data file.

Bone strength determined by destruction 3-point bend testing in 100 unselected knockout mouse strains. The mean value for each parameter obtained from n = 77 female 16 week old WT mice is shown as a horizontal line with the 2.0SD reference range limits shaded in grey. Black dots represent values for individual knockout mice. Grey vertical boxes highlight the distributions of values from mice of a single strain in which at least one animal lies outside the reference range, but the mean value lies within the reference range. Violet vertical boxes and red squares highlight the distributions of individual values from outlier strains in which the mean value (horizontal line in violet box) lies outside the reference range. Gene symbols for outlier strains are indicated. A–F, Yield load, maximum load, fracture load, stiffness, the proportion of energy dissipated at maximum load, and the proportion of energy dissipated at fracture.

(TIF)

Click here for additional data file.

Relationship between bone structure and strength. In each graph individual wild-type mice are identified by black spots. The 10 strains with major phenotypes are identified by red squares and gene symbols. The 2.0SD reference range for each variable is represented by the grey box. Pearson correlation coefficients for significant relationships are indicated and blue lines demonstrate significant linear correlations. A–C, Biomechanical and Faxitron parameters. Scatter graphs showing each of the six biomechanical parameters (yield load, maximum load, fracture load, stiffness, proportion of energy dissipated at maximum load, proportion of energy dissipated at fracture) plotted versus each of the three structural parameters (median grey value, cortical thickness, bone length) obtained by Faxitron x-ray microradiography. D–F, Biomechanical and micro-CT parameters. Scatter graphs showing each of the six biomechanical parameters (yield load, maximum load, fracture load, stiffness, proportion of energy dissipated at maximum load, proportion of energy dissipated at fracture) plotted versus each of the three structural parameters (BV/TV, Tb. N, Tb. Th) obtained by micro-CT. G–I, Faxitron and micro-CT parameters. Scatter graphs showing each of the three structural parameters obtained by micro-CT versus each of the three structural parameters obtained by Faxitron x-ray microradiography.

(TIF)

Click here for additional data file.

Phenotype data for all 100 knockout strains. Mean values (±SD) for all phenotype parameters were determined from n = 77 WT mice and are highlighted in yellow in the table header. The mean values (and distance in SD from the C57BL/6 reference mean) are shown for each parameter determined in each of the 100 mouse strains. Bone volume/tissue volume (BV/TV), trabecular thickness (Tb.Th) and trabecular number (Tb. N) were determined by micro-CT and bone mineral content (BMC), mean length of tibia and ulna, and cortical thickness (C.Th) by Faxitron x-ray microradiography. Mahalanobis distances were calculated to account for variance in all six structural parameters simultaneously and strains with significant values are indicated. Yield load, maximum load, fracture load, stiffness, the proportion of energy dissipated at maximum load and proportion of energy dissipated at fracture were determined by destruction 3-point bend testing. Note all SD values were rounded to 1 decimal place and only values >2.0SD outside the reference mean for each parameter are highlighted in grey. Measures of skeletal morphology (x-ray), body length, bone mineral density (DEXA) and serum biochemical parameters (Biochem) of relevance to bone and mineral metabolism (Ca2+, Mg2+, PO4 3−, creatinine and alkaline phosphatase) were determined as part of the broad primary phenotype screen conducted by the MGP. Primary phenotype parameters for which ≥60% of knockout mice were outside the 95% reference range are highlighted in grey (http://www.sanger.ac.uk/mouseportal/). N = normal; ND = not determined; ↓ = reduced. Normal expression of deleted genes in bone and cartilage was determined by whole mount staining for lacZ reporter gene expression in heterozygous mice as part of the primary phenotype screen and the outcome of these studies is also shown. Note that quality control analysis of gene targeting events performed after completion of these studies identified two strains (Psmb2 and 5031439G07Rik) that were incorrectly targeted and did not harbour a knockout allele. Accordingly, neither strain had an abnormal skeletal phenotype.

(XLSX)

Click here for additional data file.

Accuracy of methods to identify major phenotypes affecting the structure and strength of bone. The sensitivities, specificities and positive and negative predictive values of x-ray microradiography, micro-CT, Mahalanobis calculations and broad primary phenotyping are shown.

(DOC)

Click here for additional data file.