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

Summary

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.

Keywords:

Genetic loci – Genome-wide association studies – Mouse models – Alizarin staining – Bone fracture – osteoporosis – Osteoblast differentiation – Osteoblasts

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 (h²>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 (h²) of mineralization was 0.87 (P<2.2 x 10⁻¹⁶).

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⁻⁶ (Fig 2A). The three associations were independent (r² 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].

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⁻⁶) 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, r²>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).

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).

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 Wnt4^fl/fl Prrx1-Cre mice at maturity; however, there were reductions in trabecular number at the femoral metaphysis in both female and male Wnt4^fl/fl femurs (S1 Fig). Furthermore, female Wnt4^fl/fl mice presented with reductions in femoral areal bone mineral density (aBMD) and cortical area compared to Wnt4^wt/wt mice (S1 Fig).

To further investigate the function of Zbtb40, a mouse model was created using the CRISPR/Cas9 technology, generating a strain (Zbtb40^mut/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 Zbtb40^mut/mut mice differentiated in vitro show reduced ALP staining and activity as well as reduced Alizarin Red staining compared to osteoblasts from Zbtb40^wt/wt mice (Fig 5C–5E). As was observed in our siRNA knockdown studies, the Zbtb40^mut/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 Zbtb40^mut/mut compared to Zbtb40^wt/wt mice, we found no changes in femoral areal BMD (Fig 5H) but did identify a significant reduction in lumbar spine areal BMD in Zbtb40^mut/mut male mice (Fig 5I). As measured by micro CT, trabecular bone in the L5 vertebrae from male Zbtb40^mut/mut mice exhibited reductions in bone volume/total volume (BV/TV), trabecular number (Tb.N), and connectivity density [Conn. Dens] compared to Zbtb40^wt/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.

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-Wnt4^tm1.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⁺_Itpr3^tf/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/cm² 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/cm² 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 dH₂0 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 dH₂0 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).

Zbtb40^mut strain generation

Zbtb40^mut 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 Zbtb40^mut transcript. The translated ZBTB40^mut 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 5^th 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 [a]

S2 Fig [tif]

S1 Table [xlsx]

S2 Table [xlsx]

S1 Data [xlsx]