Genetic analysis of osteoblast activity identifies Zbtb40 as a regulator of osteoblast activity and bone mass

Abstract

Osteoporosis is a genetic disease characterized by progressive reductions in bone mineral density (BMD) leading to an increased risk of fracture. Over the last decade, genome-wide association studies (GWASs) have identified over 1000 associations for BMD. However, as a phenotype BMD is challenging as bone is a multicellular tissue affected by both local and systemic physiology. Here, we focused on a single component of BMD, osteoblast-mediated bone formation in mice, and identified associations influencing osteoblast activity on mouse Chromosomes (Chrs) 1, 4, and 17. The locus on Chr. 4 was in an intergenic region between Wnt4 and Zbtb40, homologous to a locus for BMD in humans. We tested both Wnt4 and Zbtb40 for a role in osteoblast activity and BMD. Knockdown of Zbtb40, but not Wnt4, in osteoblasts drastically reduced mineralization. Additionally, loss-of-function mouse models for both genes exhibited reduced BMD. Our results highlight that investigating the genetic basis of in vitro osteoblast mineralization can be used to identify genes impacting bone formation and BMD.

About the Authors:

Madison L. Doolittle

Contributed equally to this work with: Madison L. Doolittle, Gina M. Calabrese

Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing

Affiliation: Center for Musculoskeletal Research, University of Rochester, Rochester, New York, United States of America

ORCID logo http://orcid.org/0000-0003-0912-0095

Gina M. Calabrese

Contributed equally to this work with: Madison L. Doolittle, Gina M. Calabrese

Roles Formal analysis, Investigation

Affiliation: Center for Public Health Genomics, University of Virginia, Charlottesville, Virginia, United States of America

Larry D. Mesner

Roles Formal analysis, Investigation, Methodology

Affiliations Center for Public Health Genomics, University of Virginia, Charlottesville, Virginia, United States of America, Department of Public Health Sciences, University of Virginia, Charlottesville, Virginia, United States of America

Dana A. Godfrey

Roles Conceptualization, Investigation, Methodology

Affiliation: Department of Orthopedics, University of Colorado, Aurora, Colorado, United States of America

ORCID logo http://orcid.org/0000-0002-5092-8099

Robert D. Maynard

Roles Formal analysis, Methodology

Affiliations Center for Musculoskeletal Research, University of Rochester, Rochester, New York, United States of America, Department of Orthopedics, University of Colorado, Aurora, Colorado, United States of America

ORCID logo http://orcid.org/0000-0002-7590-2442

Cheryl L. Ackert-Bicknell

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected] (CLAB); [email protected] (CRF)

¶‡ These authors also contributed equally to this work.

Affiliations Center for Musculoskeletal Research, University of Rochester, Rochester, New York, United States of America, Department of Orthopedics, University of Colorado, Aurora, Colorado, United States of America

ORCID logo http://orcid.org/0000-0003-4253-4775

Charles R. Farber

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Software, Supervision, Validation, Writing – original draft, Writing – review & editing

* E-mail: [email protected] (CLAB); [email protected] (CRF)

¶‡ These authors also contributed equally to this work.

Affiliations Center for Public Health Genomics, University of Virginia, Charlottesville, Virginia, United States of America, Department of Public Health Sciences, University of Virginia, Charlottesville, Virginia, United States of America, Department of Biochemistry and Molecular Genetics, University of Virginia, Charlottesville, Virginia, United States of America

ORCID logo http://orcid.org/0000-0002-6748-4711

Introduction

Osteoporosis is a metabolic disease characterized by progressive bone loss leading to skeletal fragility and fracture [1]. Approximately 200 million people worldwide have or are at risk of developing osteoporosis [2], leading to ~8.6 million fractures annually [3]. As the proportion of aged persons worldwide is increasing, osteoporosis is becoming an even greater public health burden [4]. Importantly, one in three women and one in five men will suffer an osteoporotic fracture [5–7]. Family history remains the strongest risk factor for development of osteoporosis and studies in animal models reinforce that this is a complex genetic disease [8–10]. Furthermore, fracture-related traits, such as bone mineral density (BMD), are among the most heritable disease associated quantitative traits (h2>0.50) [10–12]. Thus, increasing our understanding of the genes influencing osteoporosis is critical for the development of approaches for its treatment and prevention.

Bone is not a static tissue, but rather it is constantly changing to adapt to the current needs of the organism. Turnover in bone is accomplished by three main cell types: the osteoblasts which are responsible for forming new bone, the osteoclasts, which are the cells responsible for bone resorption, and the osteocytes which orchestrate bone turnover through paracrine regulation of both the aforementioned cell types [13–15]. During normal remodeling, bone formation equals resorption and net bone mass is constant. When remodeling is imbalanced, such as occurs in osteoporosis, bone resorption exceeds formation leading to a net loss of bone and a decrease in BMD. Primary imbalances in bone remodeling occur due to hypogonadism, such as in menopause in women [16], and with age [17]. Secondary bone remodeling imbalances can result from numerous conditions which alter systemic processes such as hormone regulation, or nutritional intake and absorption. These conditions include, but are not limited to, Type 2 Diabetes [18], inflammatory bowel disease [19], anorexia nervosa [20], and chronic kidney disease [21]. Furthermore, several environmental factors such as smoking [22], air pollution [23], alcohol consumption [24], high fat diets [25] and use of common medications such as proton pump inhibitors [26], barbiturates [27], glucocorticoids [28], and nonsteroidal anti-inflammatory drugs (NSAIDS) [29] are associated with reduced bone mass and/or fracture incidence. Also, bone mass is influenced by genetic signals that emanate from a plethora of organ systems [30]. While it is important to identify all factors that influence bone mass, the complexity of BMD means that any individual genetic association may exert its influence through many potential mechanisms, confounding our understanding of the associated biology.

Genome-wide association studies (GWAS) have been extremely successful in identifying loci affecting osteoporosis-related traits. These studies have primarily focused on BMD and surrogates for BMD, such as estimated BMD (eBMD) measured using ultrasound of the calcaneus [31]. To date, the largest GWAS for eBMD identified 1103 independent associations, collectively explaining ~20% of the variance for this phenotype [32], emphasizing the idea that much remains to be identified. Additionally, mouse mapping studies have identified large numbers of quantitative trait loci (QTL) [8] and genome-wide associations for bone traits [33].

Reverse genetic studies provide valuable information about the function of a gene in physiology and provide an insight into pathology. For this reason, reverse genetic approaches remain unreplaceable for the study of individual genes, but there is no consensus on which genes at which loci to prioritize for follow up and why one should focus on one gene over another gene which may also be near to the lead single nucleotide polymorphisms (SNP). While all genetic influences for osteoporosis will be eventually identified, there is an urgent need for studies that tell us how aspects of BMD such as the regulation of bone acquisition, maintenance and loss are affected. While there are effective drug therapies that inhibit bone resorption, these are not without side effects and there are limited anabolic therapies for osteoporosis and rare bone diseases [34].

The goal of this study was to identify genetic regulators of a specific component of bone remodeling, namely mineralization by the osteoblast. After confirming that osteoblast mineralization is indeed heritable, we used this phenotype to performed the first cell-specific genome wide association study of osteoblast function. Our study was designed to minimize environmental variation and eliminate extrinsic physiological influencers that are known to act on the osteoblast. We identified three loci for this phenotype on Chromosomes 1, 4 and 17 respectively and determined that these loci were all concordant with previously identified loci for estimated BMD [32]. We then used reverse genetic methods to confirm that the previously under-characterized gene, Zbtb40, played a role osteoblast mineralization ex vivo. This work shows that our methods is efficient for finding genes that impact osteoblast function.

Results

Genetic analysis of osteoblast activity

Measuring osteoblast activity in a panel of inbred mouse strains.

One of the primary goals of this study was to “simplify” the genetic analysis of complex skeletal traits by focusing on a cellular phenotype. Primary osteoblasts isolated from neonatal calvariae, when cultured in osteogenic media, form “nodules” of mineralization that are similar to the mineralized matrix found in bone and in vitro mineralization directly reflects the primary activity of bone forming osteoblasts. We began by isolating primary calvarial osteoblasts from 22 inbred strains of mice, differentiated them into mineralizing osteoblasts, and measured the amount of mineral after 10 days in culture (S1 Table). We observed high levels of inter-strain variation in mineralization (Fig 1A), with a 3.9-fold increase between the strain with the lowest (BALB/cJ) and highest (C3H/HeJ) mineralization values. The narrow-sense heritability (h2) of mineralization was 0.87 (P<2.2 x 10−16).

[Figure omitted. See PDF.]

Fig 1. In vitro osteoblast mineralized nodule formation is a heritable trait and correlates with in vivo BMD.

(A) Quantification of Alizarin Red-stained cultures of calvarial osteoblasts isolated from 22 inbred mouse strains after 10 days of osteogenic differentiation. Data are presented as mean ± S.E.M. (B) Alizarin Red values correlate significantly with areal BMD in mouse strains (r = 0.54 p = 0.03; BMD data from Ackert1 [35], https://phenome.jax.org/).

https://doi.org/10.1371/journal.pgen.1008805.g001

We previously generated total body BMD data for 16 of the 22 strains (http://phenome.jax.org, Project: Ackert1, Dataset:25011) [35]. As would be expected, mineralization and BMD were positively correlated (r = 0.54; P = 0.03, Fig 1B). These data indicate that in vitro osteoblast mineralization is highly heritable and confirm that variation in osteoblast activity is an important contributor to BMD.

Genome-wide association study for osteoblast activity.

We next sought to identify genome-wide associations influencing osteoblast activity. To accomplish this, associations between mineralization and ~196K SNPs across the 22 strains were calculated using the Efficient Mixed Model Algorithm (EMMA, [36]). Three loci on Chromosomes (Chrs) 1, 4, and 17 exceeded the permutation-derived significance threshold of P = 2.5 x 10−6 (Fig 2A). The three associations were independent (r2 among all top SNPs <0.25). Osteoblasts from strains homozygous for the reference (C57BL/6J) allele on the Chrs. 1 and 17 associations were less active, whereas osteoblasts from strains homozygous for the reference allele on Chr. 4 were more active (Fig 2B, 2C and 2D). Mouse QTL for BMD and other bone traits have been found overlapping these associations, providing further support that they are true associations [8].

[Figure omitted. See PDF.]

Fig 2. Genome-wide association analysis identifies three loci for osteoblast activity.

[A) Association results (-log10 of the P-value) for 196,330 SNPs and Alizarin Red quantification values across 22 inbred mouse strains. The dotted black line represents the genome-wide significant threshold. (B) Effect plots for the lead SNPs at each of the three significant loci. Ref = reference (B6) allele; alt = non-reference allele.

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

Characterization of GWAS loci.

For the purpose of characterizing each locus, we conservatively defined the location of each association as the region spanning the most significantly associated haplotype flanked by 500 Kbp upstream and downstream. Based on this definition, the associations spanned 1.0 (Chr. 1), 1.2 (Chr. 4), and 2.1 Mbp (Chr. 17). We first evaluated whether any of the most significantly associated SNPs or SNPs in high LD (r2>0.8) were potentially functional (missense, frameshift, etc.). There were no genes, however, across the three regions that harbored such a variant.

We then used two expression datasets to prioritize genes: RNA-seq time-course profiles [days 0–18, every 2 days, zero aligned reads at all time point was considered not-expressed] from flow assisted cell sorting (FACs) sorted calvarial osteoblasts isolated from C57BL/6J mice (GSE54461, previously described in detail in [37] and local eQTL in demarrowed cortical bone GSE27483, previously described in detail in [33, 38] in the Hybrid Mouse Diversity Panel (HMDP), a set of 100 inbred strains, including 13 of the strains used in our pilot [39]). Of the genes determined to be expressed in osteoblasts across differentiation (data set GSE54461), Jarid1b, 9530009M10Rik, and Ptprv on Chr 1 and Zbtb40, A430061O12Rik, 2810405F17Rik, Rap1gap, and Ece1 on Chr 4 were regulated by local eQTL (P<1.0 x 10−6) in cortical bone (S2 Table). However, of these, only the lead SNP for the Ptprv eQTL on Chr 1 was in strong linkage disequilibrium (LD, r2>0.8) with the lead SNP for the mineralization association. Specifically, the lead SNP, rs3086265, was the same for both associations and the allele of rs30862651 associated with lower mineralization, was associated with lower levels of Ptprv in bone and this is consistent with the known role of Ptprv in osteoblast activity (other synonyms for this gene name include Esp, mOST-PTO and OST, [40, 41]). Together these data suggested that the Ptprv local eQTL was responsible for the Chr 1 association.

Mineralization associations overlap with human BMD loci.

We next queried the human regions syntenic with the mouse associations for evidence of association with BMD. Overlap would suggest convergence on the same target genes impacting osteoblast activity in mice and BMD in humans. To accomplish this, we used data from the largest GWAS for heel estimated BMD (eBMD) performed to date (N~426K individuals) using data from the UK BioBank [32]. The three human regions syntenic with the mouse associations harbored a total of 12 genome-wide significant associations (Fig 3). The human region syntenic with the mouse Chr. 1 mineralization locus harbored an association influencing eBMD overlapping LGR6 and PTPRV (Fig 3A). The human region syntenic with the mouse Chr. 4 mineralization locus contained eight independent associations, most (six) of which were centered in the intergenic region between ZBTB40 and WNT4 in both species (Fig 3B). We also identified three human eBMD associations in the region syntenic with the mouse Chr. 17 mineralization locus, one centered over PKDCC (also known as Adtk1 and Vlk) and the other overlapping the ABCG5/8 cluster (Fig 3C).

[Figure omitted. See PDF.]

Fig 3. Human regions syntenic with mouse associations harbor genome-wide significant associations for BMD.

Plots of GWAS results for mineralization associations on mouse Chrs. 1 (A), 4 (B), and 17 (C). eBMD associations [32] overlapping the human genomic regions syntenic with the three osteoblast mineralization loci on human Chrs. 1 (D), 1 (E), and 2 (F). Note the orientation is reversed in humans for the mouse associations on Chrs. 1 and 4. The horizontal solid black line in each plot is the genome-wide significant threshold. The dashed black vertical lines in panels D, E, and F denote the locations of independent eBMD loci in humans. The vertical grey solid lines denote the approximate peak location of the lead mouse mineralization SNP in humans. SNPs in all plots are color-coded based on LD with the most significant SNP in the genomic interval.

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

Validation of candidate genes flanking chr. 4 locus

Given that the immediate region in humans around ZBTB40—WNT4 harbored six eBMD GWAS associations and the locus has also been associated with fracture [31, 32, 42–46], the Chr. 4 association in the mouse syntenic region provided grounds for further investigation. In mice, both Zbtb40 and Wnt4 were located within or flanking the boundaries of the mineralization locus (same LD block) and in humans both genes were located in the LD block of at least one eBMD association (Fig 3). Therefore, we predicted that either Wnt4 or Zbtb40 was the causal gene underlying the association for this Chr. 4 locus and sought to experimentally validate the putative roles of these two genes in osteoblast mediated mineralization of the bone matrix.

Knockdown of Zbtb40, but not Wnt4, disrupts in vitro osteoblast mineralized matrix formation.

To investigate the effect of Wnt4 knockdown on mineralization, we isolated primary calvarial osteoblasts from mice with Wnt4 fl/fl, wt/fl, and wt/wt genotypes. To delete Wnt4, the cells were then transfected in vitro with cre-expressing or empty control plasmids. Using this approach, approximately 38% (P<0.001) and 20%(P<0.001) of the Wnt4 floxed alleles were excised in fl/fl and wt/fl cells, respectively (Fig 4A). We then quantified mineralized nodule formation at 10 days post-differentiation by Alizarin Red staining and found no significant difference in mineralization as a function of Wnt4 genotype (P = 0.2077; one-way ANOVA, Fig 4B and 4C).

[Figure omitted. See PDF.]

Fig 4. Knockdown of Zbtb40, but not Wnt4, inhibits osteoblast differentiation and mineralization.

(A) Relative expression of Wnt4 mRNA levels in calvarial osteoblasts from differing genotypes isolated and transfected with a cre-expressing plasmid in vitro. (B) Quantification of Alizarin Red-stained cultures of Wnt4 knockdown (wt/wt n = 12, wt/fl n = 45, fl/fl n = 21). (C) Immunoblot image and (D) quantification of ZBTB40 normalized protein levels in MC3T3 preosteoblast cells indicating knockdown. (E) Relative expression of Zbtb40 mRNA levels measured by qRT-PCR. (F) Quantification of Alizarin Red-stained Control (Ctrl) or Zbtb40 siRNA cultures. (G) Cell monolayer images of siRNA-treated cultures stained for Alkaline Phosphatase (ALP) or Alizarin Red measured at day 4 and day 21, respectively. (H) Quantified ALP per well measured by colorimetric activity assay. (I) Relative Expression of mRNA levels of Msx2 at day 0. Relative mRNA levels of early (J) and late (K) osteogenic factors measured at day 4 and day 21, respectively. n = 3 independently performed experiments. +++ = p<0.0001 –One-way ANOVA with Dunnett’s, * = p<0.05 ** = p<0.01 *** = p<0.001 –Unpaired T test.

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

As Wnt4 loss-of-function appeared to have no intrinsic effect on osteoblast mineralization, we next looked to Zbtb40. For functional studies, siRNA knockdown of Zbtb40 was performed in MC3T3-E1 Subclone 4 pre-osteoblast cells [47]. After delivery of either scrambled control [Ctrl] or Zbtb40 siRNA (Fig 4D–4F) the cells were treated with osteogenic media and studied for markers of osteoblast differentiation. As these cells matured, the Zbtb40 siRNA-treated cells exhibited a complete absence of mineralized nodule formation, as measured by Alizarin Red staining (Fig 4G and 4H). In the Zbtb40 siRNA group, there was a reduction in staining intensity and activity (Fig 4H and 4I) of alkaline phosphatase [ALP], a marker of osteoblast differentiation. siRNA knockdown of Zbtb40 led to reduced mRNA expression of the early transcription factor Msx2 (Fig 4J), a factor that stimulates mesenchymal progenitor cell fate commitment towards the osteoblast lineage [48]. Furthermore, there were reductions in transcript expression of the early osteoblast marker, Col1a1, as well as osteogenic transcription factors Runx2 and Sp7 (Fig 4K). Upon maturation, there was also a reduction in expression of Bglap, the gene that encodes for the mature osteoblast marker osteocalcin (Fig 4I).

Genetic manipulation of Wnt4 and Zbtb40 in vivo lead to reductions in bone mass in mice.

We generated a loss-of-function mouse model to further test the effects of Wnt4 deletion in the intact bone environment. Wnt4 global knockout mice have been documented to die within 24 hours of birth due to disrupted kidney function [49]. Thus, a conditional knockout mouse model was created using the Cre-Loxp system [50], with the Cre recombination driven by the Prrx1 promoter to achieve targeted Wnt4 deletion in early osteochondro-progenitors [51]. Most bone parameters were unaffected in the Wnt4fl/fl Prrx1-Cre mice at maturity; however, there were reductions in trabecular number at the femoral metaphysis in both female and male Wnt4fl/fl femurs (S1 Fig). Furthermore, female Wnt4fl/fl mice presented with reductions in femoral areal bone mineral density (aBMD) and cortical area compared to Wnt4wt/wt mice (S1 Fig).

To further investigate the function of Zbtb40, a mouse model was created using the CRISPR/Cas9 technology, generating a strain (Zbtb40mut/mut) homozygous for a mutant or truncated form of the ZBTB40 protein. Specifically, this allele resulted in a protein product which lacked the “BTB” protein-protein interaction domain, emulating a partial loss-of-function (Fig 5A and 5B). Calvarial osteoblasts from Zbtb40mut/mut mice differentiated in vitro show reduced ALP staining and activity as well as reduced Alizarin Red staining compared to osteoblasts from Zbtb40wt/wt mice (Fig 5C–5E). As was observed in our siRNA knockdown studies, the Zbtb40mut/mut osteoblasts exhibited reduced transcript expression of Col1a1, Runx2, and Sp7 early in differentiation (Fig 5F), as well reduced Bglap upon reaching maturity (Fig 5G). In analyzing the site-specific skeletal phenotypes of Zbtb40mut/mut compared to Zbtb40wt/wt mice, we found no changes in femoral areal BMD (Fig 5H) but did identify a significant reduction in lumbar spine areal BMD in Zbtb40mut/mut male mice (Fig 5I). As measured by micro CT, trabecular bone in the L5 vertebrae from male Zbtb40mut/mut mice exhibited reductions in bone volume/total volume (BV/TV), trabecular number (Tb.N), and connectivity density [Conn. Dens] compared to Zbtb40wt/wt mice (Fig 5J–5M). Investigating the femur in more detail, we found no differences between genotypes of either sex in cortical bone morphology at the femoral mid-diaphysis, or any differences in cortical area or thickness (S2 Fig). Femur cortical bone strength was measured by performing three-point biomechanical bending to failure. There were no differences in maximum load or calculated stiffness (S2 Fig) suggesting that there was no impact of this mutant allele on bone quality.

[Figure omitted. See PDF.]

Fig 5.

Mutation of ZBTB40 protein leads to disrupted osteogenesis and reduced vertebral bone mass in male mice: (A) Diagram of ZBTB40 protein from WT and Zbtb40mut genotypes, illustrating the truncation of the BTB functional domain. (B) Immunoblot for ZBTB40 indicating the lower molecular weight of the mutant protein (~8kDa loss). (C) Calvarial osteoblasts isolated from WT and ZBTB40mut neonatal mice stained for ALP or with Alizarin Red. (D) Quantified ALP and (E) Alizarin Red values. (F) Relative mRNA levels of early and (G) late osteoblast markers taken at day 2 and day 14, respectively. (H) Bone Mineral Density (BMD) of femurs or (I) L4-L6 vertebrae from WT or ZBTB40mut mice at 16 weeks of age (male wt/wt n = 14, fl/fl n = 9, female wt/wt n = 9, fl/fl n = 9). (J) 3D reconstructions and calculated (K) Bone Volume/Total Volume (BV/TV), (L) Connectivity Density (Conn. Dens), and (M) Trabecular Number (Tb.N) of trabecular bone in the vertebral body of the L5 vertebrae as measured by MicroCT (wt/wt n = 9, fl/fl n = 9) * = p<0.05 ** = p<0.01 *** = p<0.001—Unpaired T test.

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

Discussion

The traditional phenotypes used to determine the genetic etiology of osteoporosis and low bone mass are BMD, which is based on X-ray attenuation by the bone tissue, or a surrogate estimate of BMD obtained via ultrasound (eBMD). Both techniques yield information about the result of a biological process, but cannot inform on what parts of that process are impacted. Thus, GWAS conducted on these phenotypes can yield loci associated with many different aspects of physiology, making experimental validation of candidate genes challenging. In this study, we chose a different route to identify loci associated with bone physiology, namely the phenotype of mineralization by the osteoblast. In this study, we established that in vitro osteoblast mineralized matrix formation is a heritable phenotype. To attempt to find the genetic determinants affecting osteoblast function, we identified three loci associated with the in vitro phenotype of osteoblast mineralization, on mouse Chrs. 1, 4 and 17. All three human regions syntenic with the mouse associations harbored associations with eBMD. As bone formation by the osteoblast is contributory to BMD, it is likely that that candidate genes for the osteoblast function loci are also contributory to the concordant BMD loci.

The Chr. 1 locus, which is located at 136.6 Mb (mm9), has been previously identified in a large number of genetic mapping studies in mice [8]. Further follow up studies using congenic mice have suggested that there are actually multiple loci on the distal end of this chromosome [52]. Our locus overlaps with one of these locus, which is referred to in the mouse genome informatic database (MGI) as the quantitative trait locus (QTL), bmd5 [53]. This locus was found using a series of nested congenic mice for which C3H/HeJ (C3H) alleles were introgressed on a C57BL/6J background. Specifically, the c3h alleles conferred higher volumetric BMD of the femur but had no impact on femur length [53]. Based on the complex phenotypes of the various congenic lines, Beamer and colleagues concluded that this region of mouse Chr 1 contains multiple small effect size loci for BMD [52]. No candidate gene for BMD was found for bmd5 using the congenic approach. More recently a GWAS loci in the proximity of PTPRVP was identified for the phenotype of eBMD [32]. This gene is located in the region directly homologous with our region of interest [Fig 3A]. Interestingly, this gene is expressed at low levels in the maturing osteoblast (Gene expression omnibus series: GSE54461, [37]), but no expression was noted in the osteoclast (Gene expression omnibus accession number: GSM1873361). Further, the mouse Ptprv gene is regulated by a local eQTL, providing evidence of causality for this gene, for this locus.

The Chr. 17 locus interval encompasses a number of genes. Of specific interest is Pkdcc, as numerous GWAS in humans have identified this gene as a candidate for loci associated with whole body, femoral neck and heel ultrasound BMD [32]. Mice lacking this gene, which is also known as Adtk1, present with a large number of skeletal defects including craniofacial deformities, hypo-mineralization of long bones, and limb shortening [54]. Further, this gene is moderately expressed in the osteoblast (Gene expression omnibus series: GSE54461, [37]). Collectively, these data support Pkdcc as a likely candidate for this Chr 17 locus.

As the candidate gene underlying the Chr 4 association remained unclear, an experimental approach was used to study this locus. The absence of an effect of Wnt4 knockdown on in vitro mineralized matrix formation was an unexpected result. The WNT signaling pathway is well described in the literature to have a major role in osteogenesis and numerous WNT family members (Wnt16, Wnt3a, Wnt5a) [55–57] and canonical signaling proteins (β-catenin, LRP5/6, Dkk1/2, SOST) [58–60] have been shown specifically to be involved in this process. Moreover, others have shown that overexpression of Wnt4 in osteoblasts protects against several mechanisms of bone loss [61]; however, loss-of-function studies have been lacking. Although we did not see any effect of Wnt4 knockdown on in vitro osteoblast cultures, we recognize that this cell culture system removes the osteoblast from the niche in which it normally resides. When observing the effects of Wnt4 deletion in osteoblasts in vivo, there were noticeable reductions in trabecular bone mass in the femur, with substantial reductions in cortical bone and aBMD in femurs of the female mice. This data conflicts with our in vitro data, as we do not see any changes in function of isolated osteoblasts, but observed an effect on whole bone tissue. The reason for these paradoxical changes in vivo may be due to the previously characterized non-canonical, paracrine effects of osteoblast-secreted Wnt4 in altering osteoclast bone resorption through inhibiting NF-κB activation [61]. While loss of Wnt4 had no effect on intrinsic osteoblast function in isolation, loss of Zbtb40 substantially disrupted innate osteoblast gene expression, differentiation, and mineralized matrix formation. While we cannot exclude either gene [or other flanking genes] as candidates for the BMD locus, our data show that Zbtb40, unlike Wnt4, is able to modulate the phenotype of osteoblast mineralization ex vivo and therefor Zbtb40 is strongly supported as a candidate for this locus.

The Zbtb40 gene is mostly uncharacterized, with little known about its function in any tissue. Zbtb40 belongs to the BTB-ZF protein family, which consist of transcription factors regulating cell commitment, differentiation, and stem cell self-renewal [62, 63]. A recent study revealed that human ZBTB40 regulates osteoblast gene transcription in cell lines, as ZBTB40 knockdown reduced mRNA expression levels of COL1A1, RUNX2, SP7, ALP and even WNT4 [64], which is confirmatory of our results in mouse cell lines. Further, putative loss-of-function of ZBTB40 protein in primary cells from our mutant mouse model and in MC-3T3 cells led to abundant reductions of mRNA expression of early osteoblast factors (Runx2, Sp7, Msx2). These factors have been well documented to induce commitment to the osteoblast lineage and drive osteogenic differentiation [65–67]. Our data shows that ZBTB40 loss alters expression of these critical drivers, suggesting that, like other members of the BTB-ZF family, this gene may have an effect on cell differentiation. Msx2, for example, regulates early mesenchymal stem cell (MSC) fate by stimulating osteogenesis [48]. The fact that Zbtb40 silencing reduces Msx2 expression demonstrates a role for this gene in MSCs. However, more investigation is needed to more fully establish the function of and signaling pathway(s) acted on by Zbtb40.

In summary, we have conducted a GWAS for osteoblast function in vitro using cells isolated from inbred strains of mice. We have shown that these loci have high concordance with human GWAS loci for eBMD. Because our phenotype was restricted to single physiological process by a single cell type, our mouse-human concordant loci become informative for the human BMD loci as we are able to provide information regarding a mechanism by which given loci may be able to impact BMD. Our data suggest that genetic studies focused on osteoblast activity have the ability to provide significant insight into the genetic and molecular basis of bone formation and osteoporosis.

Methods

Ethics statement

All animal procedures were performed according to protocols reviewed and approved by The University Committee for Animal Resources (UCAR) Institutional Animal Care and Use Committee (IACUC) at the University of Rochester (Protocol number: 101717). Where possible, data from public repositories was used to reduce the number of animals used in this study. The forms of euthanasia that were used were done so in accordance with the guidance and recommendations provided by the Panel on Euthanasia of the American Veterinary Medical Association.

Animal models

All mice used in this study, with the exception of the Zbtb40 loss of function allele mice (described in detail below), were purchased from The Jackson Laboratory. The mice were fed Laboratory Autoclavable Rodent Diet 5010 (LabDiet, Cat no. 0001326) and had ad lib. access to food and water. The mice were maintained on a 12hr:12hr light:dark cycle. The stock numbers for the strains used were as follows: B6;129S-Wnt4tm1.1Bhr/BhrEiJ (Wnt4 floxed mice, #007032, [68]), B6.Cg-Tg(Prrx1-cre)1Cjt/J (Prrx1-cre driver strain, # 005584), BALB/cJ (#000651), DBA/2J (#000671), FVB/NJ (#001800), RIIIS/J (#000683), BTBR_T+_Itpr3tf/J (#002282), SJL/J (#000686), NZW/LacJ (#001058), C57BL/6J (#00064), BPN/3J (#003004), TALLYHO/JngJ (#005314), NZO/HlLtJ (#002105), NON/ShiLtJ (#002423), NOR/LtJ (#002050), A/J (#000646), LG/Jm(#000675), NOD/ShiLtJ (#001976), 129S1/SvImJ (#002448), CBA/J (#000656), DDY/JclSidSeyFrkJ (#002243), AKR/J (#000648) and C3H/HeJ (#000659). The ages and sexes of the animals used are described per experiment. Neonatal mice were euthanized by decapitation and adult mice via CO2 asphyxiation followed by cervical dislocation.

Measuring mineralization in inbred strains

Primary calvarial osteoblasts were isolated from 3–9 day old neonates (males and females pooled; N = 10–15 per pool; N = 1–6 pools per strain) from a panel of 22 inbred strains using sequential Collagenase P digestions [69]. Cells were plated into 6 well plates at 300,000 cells in 2 ml sterile plating media (DMEM, 10% heat-inactivated FBS, 100 U/ml penicillin, 100 μg/ml streptomycin) per well. After 24 hours, confluent cells were washed 1x with DPBS [GIBCO] and placed in sterile differentiation media (MEM alpha, 10% heat inactivated FBS, 100 U/ml penicillin, 100 μg/ml streptomycin, 50 μg/ml ascorbic acid, 4 mM B-glycerophosphate). Every 48 hours thereafter cells were washed one time with DPBS (GIBCO) and differentiation media was replaced until cells were collected for analysis at day 10. Mineralized nodule formation was measured by staining cultures at 10 days post-differentiation with Alizarin Red (40 mM, pH 5.6). The stained cells were imaged and nodule number was measured using ImageJ (NIH) [70]. Alizarin Red was quantified by destaining cultures with 5% Perchloric acid and determining the optical density (405 ηM) of the resulting solution against a standard curve. Alizarin Red staining was performed on three wells (technical replicates) per neonate pool [biological replicates] from each strain. The narrow-sense heritability of mineralization was calculated as described in [71] based on the results of a one-way analysis of variance (ANOVA) by dividing the between-strain sum of squares (SS) by the sum of the between-strain and within-strain SS.

Association analysis

To identify loci influencing mineralization, we used the Efficient Mixed Model Association (EMMA) algorithm [36]. For the analysis, alizarin red values for each pool per strain were rank Z transformed. SNPs were obtained from strains genotyped on the Mouse Diversity Array (http://churchill-lab.jax.org/website/MDA) [72]. SNPs with a minor allele frequency < 0.05 were removed, leaving 196,330 SNPs. These SNPs were used to generate a kinship using the ‘emma.kinship’ R script available in the EMMA R package (available at http://mouse.cs.ucla.edu/emma/) [36]. The emma.REML.t function of EMMA was used to perform all mapping analyses. Individual pool measures of Alizarin Red staining per strain were used for the EMMA analysis. The significance of the maximum association peak was assessed by performing 1,000 permutations of the data. In each permutation, the minimum p-value was recorded to produce an empirical distribution of minimum permutation p-values. The quantiles of this distribution were used to assign adjusted p-values. P-values exceeding a genome-wide significant of P<0.05 were used as thresholds to identify associated loci. GWAS resulted were visualized using the “qqman” R package [73].

Comparison between mouse mineralization loci and BMD associations in human syntenic regions

For each mouse locus (Chrs 1, 4, and 17), we conservatively defined the location of each association as the region spanning the most significantly associated haplotype flanked by 500 Kbp upstream and downstream. We then used the UCSC Genome Browser Lift-Over tool (https://genome.ucsc.edu/cgi-bin/hgLiftOver) to convert each mouse region (mm9) to its human syntenic region (hg19). Each human syntenic region was then queried for associations with BMD using data from the largest BMD (eBMD) GWAS performed to date [32].

Effect of Wnt4 knockdown on mineralization

The effect upon mineralization of the in vitro knockdown of the Wnt4 gene in differentiating primary calvarial osteoblasts was accomplished as outlined in Calabrese et al [74] with the following modifications. Primary calvarial osteoblast cells were isolated from the F2 generation of floxed Wnt4 mice, plated and differentiated exactly as outlined. DNA from the tails of the 3–9 day old mice, from which the calvaria were removed, were used for genotyping via PCR using a forward primer located 5’ of the loxP site in intron 1 (GCA GAG AGG CCC AGC CTG CCC CTC A) along with a reverse primer located 3’ of the loxP site in Intron 2 (CAT GTG CCT GGC CCT AGA AAT ATC AT; expected PCR product size: 642bp (wild type allele), 819bp (floxed allele), 239bp (recombined allele). Wnt4 knockdown was accomplished by transfecting cells with a Cre Recombinase expression plasmid (or empty vector control, [75]) 24h post-plating followed by the introduction of differentiation media 48h latter. Cells were collected for RNA and mineralization analysis 10 days after the initiation of differentiation as outlined except primers used for Wnt4 expression analysis were (GAA CTG TTC CAC ACT GGA CTC (Exon2); GTC ACA GCC ACA CTT CTC CAG (Exon3). DNA was also collected at day 10 and used to qualitatively analyze the extent of Cre-recombination activity by using the genotyping primers listed above as well as quantitative analysis via qPCR with the exon2 and 3 primers. Mineral formation was measured exactly as described and the effect of Wnt4 knockdown on mineralization was expressed as the ratio of moles of Alizarin Red bound by the Cre-transfected cells of wt/wt, wt/fl, and fl/fl genotypes.

Cell culture

Primary calvarial osteoblasts were cultured in the same manner as described above. The MC3T3-E1 mouse preosteoblast cell line (ATCC Subclone 4 CRL-2593) were cultured in ascorbic acid-free Alpha Minimal Essential Medium (α-MEM) supplemented with 10% Fetal Bovine Serum (FBS) and 100 U/ml penicillin/streptomycin. Cells were dissociated from the culture plate using 0.05% trypsin-EDTA and re-plated in experimental wells at 10,000 cells/cm2 for downstream siRNA transfection. Cells were further grown for 72 hours before reaching confluence, upon which media was changed to osteogenic medium consisting of basal culture media supplemented with 4mM β-Glycerophosphate and 50mg/ml Ascorbic Acid and differentiated until the indicated time point.

siRNA treatment

MC3T3-E1 cells were plated at 10,000 cells/cm2 in 12-well plates and allowed to adhere overnight. The following day the cells were transfected with 5nM scrambled control or Zbtb40 Silencer Select siRNA (ThermoFisher) delivered with Lipofectamine RNAiMAX (ThermoFisher). Knockdown efficiency was determined by measuring RNA (24 hours) or protein (72 hours) expression levels post-transfection. Remaining cells were then induced with osteogenic medium at 100% confluency (72 hours post-transfection) and differentiated until the indicated time point for phenotype analysis.

RNA isolation and qRT-PCR gene expression analysis

Cells were washed once with PBS and total mRNA was isolated using TRIzol (ThermoFisher) and purified using the GeneJET RNA Purification Kit (ThermoFisher). RNA was treated with DNAse I Amplification Grade (ThermoFisher) and reverse-transcribed into cDNA using the iScript cDNA Synthesis Kit (Bio-Rad). Real Time qPCR was performed using PerfeCTa SYBR Green FastMix (QuantaBio) and the Rotor-Gene Q Real-Time PCR Cycler (Qiagen). All reactions were performed in triplicate and normalized to β-2-Microglobulin using the primers F: 5’-TGACCGGCCTGTATGCTATC-3’, R: 5’- AGGCGGGTGGAACTGTGTTA-3’. Relative expression levels were calculated using the 2-ΔΔCT method [76]. Primer sets: Zbtb40 F: 5’-AGAGCCACAGCATGGAACTC-3’, R: 5’-CCGACGGAAATGGTGCAATC-3’.

Msx2 F: 5’-GATACAGGAGCCCGGCAGAT-3’, R: 5’-CTTGCGCTCCAAGGCTAGAA-3’.

Runx2 F: 5’-TGATGACACTGCCACCTCTGACTT-3’, R: 5’-ATGAAATGCTTGGGAACTGCCTGG-3’. Sp7 F: 5’-ATGGCGTCCTCTCTGCTTGA-3’, R: 5’-CTTTGTGCCTCCTTTCCCCA-3’. Col1a1 F: 5’-CGACCTCAAGATGTGCCACT-3’, R: 5’-GCAGTAGACCTTGATGGCGT-3’. Bglap F: 5’-TTCTGCTCACTCTGCTGACC-3’, R: 5’-TATTGCCCTCCTGCTTGGAC-3’.

Western blot

Cells were washed once with PBS then lysed and collected in Pierce RIPA buffer (ThermoFisher) supplemented with Protease/Phosphatase Inhibitor Cocktail (Cell Signaling). Lysates were vortexed three times and cleared by centrifugation (15 minutes at 15,000 rpm). Total protein was measured using the Pierce BSA Total Protein Kit (ThermoFisher) and samples were analyzed by SDS-PAGE in a 4–12% acrylamide gel (ThermoFisher). The gel was transferred to a nitrocellulose membrane using the iBlot 2 Dry Blotting System (ThermoFisher) and blocked in 5% milk in TBST for 1 hour at room temperature. The membrane was incubated overnight at 4° with Anti-ZBTB40 primary antibody (Abcam #Ab190185) or Anti-β-Actin Antibody (Sigma #A2228) and then incubated with either Goat Anti-Rabbit or Goat Anti-Mouse IgG-HRP Conjugate (Bio-rad) Secondary Antibody in 5% milk for 1 hour at room temperature. Membranes were developed using SuperSignal West Pico Sensitivity substrate (ThermoFisher) and imaged using a Bio-Rad Universal Hood II Chemiluminescent Imager. Bands were quantified using ImageJ and ZBTB40 protein levels were normalized to β-Actin.

Cell staining and phenotype measurements

At the indicated time point, cells were washed with PBS then fixed with 10% Neutral Buffered Formalin [NBF] for 10 minutes. For Alkaline Phosphatase (ALP) analysis, the fixed cells were stained with 1-Step NBT/BCIP Substrate Solution (ThermoFisher) in the dark for 30 minutes. Cells were then washed in dH20 and let dry. For quantitative analysis of ALP activity, separate cells were harvested and analyzed using the SensoLyte pNPP Alkaline Phosphatase Assay Kit (AnaSpec) according to manufacturer’s instructions. To detect mineralization, cells were fixed in 10% NBF and stained with Alizarin Red for 30 minutes then washed with dH20 and let dry. To obtain quantitative values for staining, plates were scanned and then de-stained with 5% Perchloric Acid. The supernatants were then collected and absorbance was measured at 405nm in a Synergy Mx Monochromator-Based Microplate Reader (Biotek).

Zbtb40mut strain generation

Zbtb40mut mice were created on a C57BL/6J background using CRISPR/Cas9 technology, as previously described [77]. Briefly, an insertion/deletion mutation was created on the DNA in the second coding exon of Zbtb40 to induce a premature stop codon, limiting full translation of the generated Zbtb40mut transcript. The translated ZBTB40mut protein is truncated at the N-terminus by 76 amino acids, lacking a majority of the BTB domain, with the remainder of the protein intact. This has been verified by both DNA and RNA sequencing as well as western blot analysis of protein products (Fig 5B). This strain has been continually backcrossed to the parental C57BL/6J strain.

Mouse phenotyping

Dual X-ray absorptiometry (DXA) was performed on all groups using the Lunar Piximus II (GE Heathcare) as described previously [52]. Areal BMD (aBMD), bone mineral content (BMC) and % body fat was measured at 16 weeks of age. BMD and BMC for the lumbar spine and femur was examined, as well as lean body mass. The left femur and 5th lumbar vertebrae were scanned using a high-resolution micro-computed tomography system (vivaCT 40, Scanco Medical AG) to assess bone architecture. Variables computed for trabecular bone regions include: bone volume, BV/TV, and trabecular number and thickness. For cortical bone, total cross-sectional area, cortical bone area, cortical thickness, and area moments of inertia about the principal axes were computed at the femoral midshaft. Femurs were tested in 3-point flexure to failure at 0.1mm/second, collecting force and displacement data at 10Hz (Bose EnduraTec 3200, Eden Prairie, MC). Specimen alignment and orientation were co-registered accurately between imaging and mechanical test configuration. Material properties were calculated using the cortical geometric indices obtained from μCT analysis.

Statistical analysis of Zbtb40 and Wnt4 experimental validation

Data are presented as the mean ± SD. Statistical significance of in vitro Wnt4 experiments was determined by one-way ANOVA followed by a Dunnett’s multiple comparison’s test using wt/wt as the control group. Statistical significance of all other experiments were determined by Student’s unpaired t test. p values <0.05 were considered significant.

Supporting information

S1 Fig.

Wnt4fl/fl Prrx1-Cre mice show reductions in femoral trabecular number in both sexes and cortical area in females: (A) Confirmation of Wnt4 deletion by floxed allele recombination in the femur of a Cre positive wt/fl mouse by RT-PCR (B) Whole body and (C) Femoral BMD of female (wt/wt n = 19, fl/fl n = 19) and male (wt/wt n = 5, fl/fl n = 17) mice measured by DXA at 16 weeks of age. (D,E,F) Trabecular number (Tb.N) values and 3D reconstructions of the femoral metaphysis from female (wt/wt n = 9, fl/fl n = 10) and male (wt/wt n = 4, fl/fl n = 9) mice. (G,H) Cortical area (Ct.Ar) calculations and 3D reconstructions of the femoral mid-diaphysis in female mice. * = p<0.05 ** = p<0.01 *** = p<0.001—Unpaired T test.

https://doi.org/10.1371/journal.pgen.1008805.s001

(TIF)

S2 Fig.

Zbtb40mut mice show no changes in femoral cortical bone mass or strength: (A, B, C) 3D reconstructions of femur mid-diaphyseal cortical bone and calculations for cortical area and cortical thickness (wt/wt n = 9, mut/mut n = 9 per sex, 16 weeks of age). (D, E) Max load and stiffness calculations from biomechanical 3-point bending test (female—wt/wt n = 9, mut/mut n = 6, male—wt/wt n = 8, mut/mut n = 9, 16 weeks of age).

https://doi.org/10.1371/journal.pgen.1008805.s002

(TIF)

S1 Table.

Quantification of mineralization by primary osteoblasts. Column A: Inbred strain name, Column B: Quantification of alizarin Red staining on a per well basis.

https://doi.org/10.1371/journal.pgen.1008805.s003

(XLSX)

S2 Table.

Local eQTL in the Hydbrid Mouse Diversity Panel: A: eQTL index. B. Probe name. C. Gene name. D. Mouse Chromosome. E. Base pair location. F. Cis Pvalue G. Cis Chromsome. H. Cis base pair location I. Chromosome Cis P value. J Chromosome Cis. K. Chromosome cis Base Pair. L. Lead Snp for the mineralization QTL. M. Lead SNP for the eQTL. N. Correlation between eQTL and mineralization QTL.

https://doi.org/10.1371/journal.pgen.1008805.s004

(XLSX)

S1 Data.

Raw data for Figs 1 through 5.

https://doi.org/10.1371/journal.pgen.1008805.s005

(XLSX)

Citation: Doolittle ML, Calabrese GM, Mesner LD, Godfrey DA, Maynard RD, Ackert-Bicknell CL, et al. (2020) Genetic analysis of osteoblast activity identifies Zbtb40 as a regulator of osteoblast activity and bone mass. PLoS Genet 16(6): e1008805. https://doi.org/10.1371/journal.pgen.1008805

References

1. Black DM, Rosen CJ. Postmenopausal Osteoporosis. N Engl J Med. 2016;374(21):2096–7.

2. Sozen T, Ozisik L, Basaran NC. An overview and management of osteoporosis. Eur J Rheumatol. 2017;4(1):46–56. pmid:28293453

3. Johnell O, Kanis JA. An estimate of the worldwide prevalence and disability associated with osteoporotic fractures. Osteoporos Int. 2006;17(12):1726–33. pmid:16983459

4. Cummings SR, Melton LJ. Epidemiology and outcomes of osteoporotic fractures. Lancet. 2002;359(9319):1761–7. pmid:12049882

5. Kanis JA, Johnell O, Oden A, Sembo I, Redlund-Johnell I, Dawson A, et al. Long-term risk of osteoporotic fracture in Malmo. Osteoporos Int. 2000;11(8):669–74. pmid:11095169

6. Melton LJ 3rd, Atkinson EJ, O'Connor MK, O'Fallon WM, Riggs BL. Bone density and fracture risk in men. J Bone Miner Res. 1998;13(12):1915–23. pmid:9844110

7. Melton LJ 3rd, Chrischilles EA, Cooper C, Lane AW, Riggs BL. Perspective. How many women have osteoporosis? J Bone Miner Res. 1992;7(9):1005–10. pmid:1414493

8. Ackert-Bicknell CL, Karasik D, Li Q, Smith RV, Hsu YH, Churchill GA, et al. Mouse BMD quantitative trait loci show improved concordance with human genome-wide association loci when recalculated on a new, common mouse genetic map. J Bone Miner Res. 2010;25(8):1808–20. pmid:20200990

9. Farber CR, Kelly SA, Baruch E, Yu D, Hua K, Nehrenberg DL, et al. Identification of quantitative trait loci influencing skeletal architecture in mice: emergence of Cdh11 as a primary candidate gene regulating femoral morphology. J Bone Miner Res. 2011;26(9):2174–83. pmid:21638317

10. Ralston SH, Uitterlinden AG. Genetics of osteoporosis. Endocr Rev. 2010;31(5):629–62. pmid:20431112

11. Boudin E, Fijalkowski I, Hendrickx G, Van Hul W. Genetic control of bone mass. Mol Cell Endocrinol. 2016;432:3–13. pmid:26747728

12. Richards JB, Zheng HF, Spector TD. Genetics of osteoporosis from genome-wide association studies: advances and challenges. Nat Rev Genet. 2012;13(8):576–88. pmid:22805710

13. Bonewald LF. The amazing osteocyte. J Bone Miner Res. 2011;26(2):229–38. pmid:21254230

14. Raggatt LJ, Partridge NC. Cellular and molecular mechanisms of bone remodeling. J Biol Chem. 2010;285(33):25103–8. pmid:20501658

15. Sims NA, Gooi JH. Bone remodeling: Multiple cellular interactions required for coupling of bone formation and resorption. Semin Cell Dev Biol. 2008;19(5):444–51. pmid:18718546

16. Eastell R, O'Neill TW, Hofbauer LC, Langdahl B, Reid IR, Gold DT, et al. Postmenopausal osteoporosis. Nat Rev Dis Primers. 2016;2:16069. pmid:27681935

17. Khosla S. Pathogenesis of age-related bone loss in humans. J Gerontol A Biol Sci Med Sci. 2013;68(10):1226–35. pmid:22923429

18. Patsch JM, Burghardt AJ, Yap SP, Baum T, Schwartz AV, Joseph GB, et al. Increased cortical porosity in type 2 diabetic postmenopausal women with fragility fractures. J Bone Miner Res. 2013;28(2):313–24. pmid:22991256

19. Szafors P, Che H, Barnetche T, Morel J, Gaujoux-Viala C, Combe B, et al. Risk of fracture and low bone mineral density in adults with inflammatory bowel diseases. A systematic literature review with meta-analysis. Osteoporos Int. 2018;29(11):2389–97. pmid:29909470

20. Mehler PS, Cleary BS, Gaudiani JL. Osteoporosis in anorexia nervosa. Eat Disord. 2011;19(2):194–202. pmid:21360368

21. Kim SM, Long J, Montez-Rath M, Leonard M, Chertow GM. Hip Fracture in Patients With Non-Dialysis-Requiring Chronic Kidney Disease. J Bone Miner Res. 2016;31(10):1803–9. pmid:27145189

22. Yan C, Avadhani NG, Iqbal J. The effects of smoke carcinogens on bone. Curr Osteoporos Rep. 2011;9(4):202–9. pmid:21874290

23. Prada D, Zhong J, Colicino E, Zanobetti A, Schwartz J, Dagincourt N, et al. Association of air particulate pollution with bone loss over time and bone fracture risk: analysis of data from two independent studies. Lancet Planet Health. 2017;1(8):e337–e47. pmid:29527596

24. Mikosch P. Alcohol and bone. Wien Med Wochenschr. 2014;164(1–2):15–24. pmid:24477631

25. Ackert-Bicknell CL, Demissie S, Marin de Evsikova C, Hsu YH, DeMambro VE, Karasik D, et al. PPARG by dietary fat interaction influences bone mass in mice and humans. J Bone Miner Res. 2008;23(9):1398–408. pmid:18707223

26. Vakil N. Prescribing proton pump inhibitors: is it time to pause and rethink? Drugs. 2012;72(4):437–45. pmid:22356286

27. Cosman F, de Beur SJ, LeBoff MS, Lewiecki EM, Tanner B, Randall S, et al. Clinician's Guide to Prevention and Treatment of Osteoporosis. Osteoporos Int. 2014;25(10):2359–81. pmid:25182228

28. Buckley L, Guyatt G, Fink HA, Cannon M, Grossman J, Hansen KE, et al. 2017 American College of Rheumatology Guideline for the Prevention and Treatment of Glucocorticoid-Induced Osteoporosis. Arthritis Rheumatol. 2017;69(8):1521–37. pmid:28585373

29. Vestergaard P, Hermann P, Jensen JE, Eiken P, Mosekilde L. Effects of paracetamol, non-steroidal anti-inflammatory drugs, acetylsalicylic acid, and opioids on bone mineral density and risk of fracture: results of the Danish Osteoporosis Prevention Study (DOPS). Osteoporos Int. 2012;23(4):1255–65. pmid:21710339

30. Karsenty G, Ferron M. The contribution of bone to whole-organism physiology. Nature. 2012;481(7381):314–20. pmid:22258610

31. Kemp JP, Morris JA, Medina-Gomez C, Forgetta V, Warrington NM, Youlten SE, et al. Identification of 153 new loci associated with heel bone mineral density and functional involvement of GPC6 in osteoporosis. Nat Genet. 2017;49(10):1468–75. pmid:28869591

32. Morris JA, Kemp JP, Youlten SE, Laurent L, Logan JG, Chai RC, et al. An atlas of genetic influences on osteoporosis in humans and mice. Nat Genet. 2019;51(2):258–66. pmid:30598549

33. Farber CR, Bennett BJ, Orozco L, Zou W, Lira A, Kostem E, et al. Mouse genome-wide association and systems genetics identify Asxl2 as a regulator of bone mineral density and osteoclastogenesis. PLoS Genet. 2011;7(4):e1002038. pmid:21490954

34. Russow G, Jahn D, Appelt J, Mardian S, Tsitsilonis S, Keller J. Anabolic Therapies in Osteoporosis and Bone Regeneration. Int J Mol Sci. 2018;20(1).

35. Mesner LD, Calabrese GM, Al-Barghouthi B, Gatti DM, Sundberg JP, Churchill GA, et al. Mouse genome-wide association and systems genetics identifies Lhfp as a regulator of bone mass. PLoS Genet. 2019;15(5):e1008123. pmid:31042701

36. Kang HM, Zaitlen NA, Wade CM, Kirby A, Heckerman D, Daly MJ, et al. Efficient control of population structure in model organism association mapping. Genetics. 2008;178(3):1709–23. pmid:18385116

37. Zheng HF, Forgetta V, Hsu YH, Estrada K, Rosello-Diez A, Leo PJ, et al. Whole-genome sequencing identifies EN1 as a determinant of bone density and fracture. Nature. 2015;526(7571):112–7. pmid:26367794

38. Calabrese G, Bennett BJ, Orozco L, Kang HM, Eskin E, Dombret C, et al. Systems genetic analysis of osteoblast-lineage cells. PLoS Genet. 2012;8(12):e1003150. pmid:23300464

39. Ghazalpour A, Rau CD, Farber CR, Bennett BJ, Orozco LD, van Nas A, et al. Hybrid mouse diversity panel: a panel of inbred mouse strains suitable for analysis of complex genetic traits. Mamm Genome. 2012;23(9–10):680–92. pmid:22892838

40. Chengalvala MV, Bapat AR, Hurlburt WW, Kostek B, Gonder DS, Mastroeni RA, et al. Biochemical characterization of osteo-testicular protein tyrosine phosphatase and its functional significance in rat primary osteoblasts. Biochemistry. 2001;40(3):814–21. pmid:11170399

41. Wheeler MA, Townsend MK, Yunker LA, Mauro LJ. Transcriptional activation of the tyrosine phosphatase gene, OST-PTP, during osteoblast differentiation. J Cell Biochem. 2002;87(4):363–76. pmid:12397596

42. Duncan EL, Danoy P, Kemp JP, Leo PJ, McCloskey E, Nicholson GC, et al. Genome-wide association study using extreme truncate selection identifies novel genes affecting bone mineral density and fracture risk. PLoS Genet. 2011;7(4):e1001372. pmid:21533022

43. Estrada K, Styrkarsdottir U, Evangelou E, Hsu YH, Duncan EL, Ntzani EE, et al. Genome-wide meta-analysis identifies 56 bone mineral density loci and reveals 14 loci associated with risk of fracture. Nat Genet. 2012;44(5):491–501. pmid:22504420

44. Rivadeneira F, Styrkarsdottir U, Estrada K, Halldorsson BV, Hsu YH, Richards JB, et al. Twenty bone-mineral-density loci identified by large-scale meta-analysis of genome-wide association studies. Nat Genet. 2009;41(11):1199–206. pmid:19801982

45. Styrkarsdottir U, Halldorsson BV, Gretarsdottir S, Gudbjartsson DF, Walters GB, Ingvarsson T, et al. Multiple genetic loci for bone mineral density and fractures. N Engl J Med. 2008;358(22):2355–65. pmid:18445777

46. Zhang L, Choi HJ, Estrada K, Leo PJ, Li J, Pei YF, et al. Multistage genome-wide association meta-analyses identified two new loci for bone mineral density. Hum Mol Genet. 2014;23(7):1923–33. pmid:24249740

47. Quarles LD, Yohay DA, Lever LW, Caton R, Wenstrup RJ. Distinct proliferative and differentiated stages of murine MC3T3-E1 cells in culture: an in vitro model of osteoblast development. J Bone Miner Res. 1992;7(6):683–92. pmid:1414487

48. Ichida F, Nishimura R, Hata K, Matsubara T, Ikeda F, Hisada K, et al. Reciprocal roles of MSX2 in regulation of osteoblast and adipocyte differentiation. J Biol Chem. 2004;279(32):34015–22. pmid:15175325

49. Stark K, Vainio S, Vassileva G, McMahon AP. Epithelial transformation of metanephric mesenchyme in the developing kidney regulated by Wnt-4. Nature. 1994;372(6507):679–83. pmid:7990960

50. Kuhn R, Schwenk F, Aguet M, Rajewsky K. Inducible gene targeting in mice. Science. 1995;269(5229):1427–9. pmid:7660125

51. Logan M, Martin JF, Nagy A, Lobe C, Olson EN, Tabin CJ. Expression of Cre Recombinase in the developing mouse limb bud driven by a Prxl enhancer. Genesis. 2002;33(2):77–80. pmid:12112875

52. Beamer WG, Shultz KL, Coombs HF 3rd, DeMambro VE, Reinholdt LG, Ackert-Bicknell CL, et al. BMD regulation on mouse distal chromosome 1, candidate genes, and response to ovariectomy or dietary fat. J Bone Miner Res. 2011;26(1):88–99. pmid:20687154

53. Shultz KL, Donahue LR, Bouxsein ML, Baylink DJ, Rosen CJ, Beamer WG. Congenic strains of mice for verification and genetic decomposition of quantitative trait loci for femoral bone mineral density. J Bone Miner Res. 2003;18(2):175–85. pmid:12568393

54. Goncalves L, Filipe M, Marques S, Salgueiro AM, Becker JD, Belo JA. Identification and functional analysis of novel genes expressed in the Anterior Visceral Endoderm. Int J Dev Biol. 2011;55(3):281–95. pmid:21553379

55. Esen E, Chen J, Karner CM, Okunade AL, Patterson BW, Long F. WNT-LRP5 signaling induces Warburg effect through mTORC2 activation during osteoblast differentiation. Cell Metab. 2013;17(5):745–55. pmid:23623748

56. Moverare-Skrtic S, Henning P, Liu X, Nagano K, Saito H, Borjesson AE, et al. Osteoblast-derived WNT16 represses osteoclastogenesis and prevents cortical bone fragility fractures. Nat Med. 2014;20(11):1279–88. pmid:25306233

57. Okamoto M, Udagawa N, Uehara S, Maeda K, Yamashita T, Nakamichi Y, et al. Noncanonical Wnt5a enhances Wnt/beta-catenin signaling during osteoblastogenesis. Sci Rep. 2014;4:4493. pmid:24670389

58. Karner CM, Long F. Wnt signaling and cellular metabolism in osteoblasts. Cell Mol Life Sci. 2017;74(9):1649–57. pmid:27888287

59. Westendorf JJ, Kahler RA, Schroeder TM. Wnt signaling in osteoblasts and bone diseases. Gene. 2004;341:19–39. pmid:15474285

60. Yavropoulou MP, Yovos JG. The role of the Wnt signaling pathway in osteoblast commitment and differentiation. Hormones (Athens). 2007;6(4):279–94.

61. Yu B, Chang J, Liu Y, Li J, Kevork K, Al-Hezaimi K, et al. Wnt4 signaling prevents skeletal aging and inflammation by inhibiting nuclear factor-kappaB. Nat Med. 2014;20(9):1009–17. pmid:25108526

62. Beaulieu AM, Sant'Angelo DB. The BTB-ZF family of transcription factors: key regulators of lineage commitment and effector function development in the immune system. J Immunol. 2011;187(6):2841–7. pmid:21900183

63. Siggs OM, Beutler B. The BTB-ZF transcription factors. Cell Cycle. 2012;11(18):3358–69. pmid:22894929

64. Mei B, Wang Y, Ye W, Huang H, Zhou Q, Chen Y, et al. LncRNA ZBTB40-IT1 modulated by osteoporosis GWAS risk SNPs suppresses osteogenesis. Hum Genet. 2019;138(2):151–66. pmid:30661131

65. Komori T, Yagi H, Nomura S, Yamaguchi A, Sasaki K, Deguchi K, et al. Targeted disruption of Cbfa1 results in a complete lack of bone formation owing to maturational arrest of osteoblasts. Cell. 1997;89(5):755–64. pmid:9182763

66. Matsubara T, Kida K, Yamaguchi A, Hata K, Ichida F, Meguro H, et al. BMP2 regulates Osterix through Msx2 and Runx2 during osteoblast differentiation. J Biol Chem. 2008;283(43):29119–25. pmid:18703512

67. Nakashima K, Zhou X, Kunkel G, Zhang Z, Deng JM, Behringer RR, et al. The novel zinc finger-containing transcription factor osterix is required for osteoblast differentiation and bone formation. Cell. 2002;108(1):17–29. pmid:11792318

68. Kobayashi A, Stewart CA, Wang Y, Fujioka K, Thomas NC, Jamin SP, et al. beta-Catenin is essential for Mullerian duct regression during male sexual differentiation. Development. 2011;138(10):1967–75. pmid:21490063

69. Dacic S, Kalajzic I, Visnjic D, Lichtler AC, Rowe DW. Col1a1-driven transgenic markers of osteoblast lineage progression. J Bone Miner Res. 2001;16(7):1228–36. pmid:11450698

70. Schneider CA, Rasband WS, Eliceiri KW. NIH Image to ImageJ: 25 years of image analysis. Nat Methods. 2012;9(7):671–5. pmid:22930834

71. Lariviere WR, Mogil JS. The genetics of pain and analgesia in laboratory animals. Methods Mol Biol. 2010;617:261–78. pmid:20336428

72. Yang H, Ding Y, Hutchins LN, Szatkiewicz J, Bell TA, Paigen BJ, et al. A customized and versatile high-density genotyping array for the mouse. Nat Methods. 2009;6(9):663–6. pmid:19668205

73. Turner SD. qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. bioRxiv Cold Spring Harbor Laboratory. 2014;005165.

74. Calabrese GM, Mesner LD, Stains JP, Tommasini SM, Horowitz MC, Rosen CJ, et al. Integrating GWAS and Co-expression Network Data Identifies Bone Mineral Density Genes SPTBN1 and MARK3 and an Osteoblast Functional Module. Cell Syst. 2017;4(1):46–59 e4. pmid:27866947

75. Mesner LD, Hamlin JL. Specific signals at the 3' end of the DHFR gene define one boundary of the downstream origin of replication. Genes Dev. 2005;19(9):1053–66. pmid:15879555

76. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25(4):402–8. pmid:11846609

77. Canalis E, Yu J, Schilling L, Yee SP, Zanotti S. The lateral meningocele syndrome mutation causes marked osteopenia in mice. J Biol Chem. 2018;293(36):14165–77. pmid:30042232