Genomic variants within chromosome 14q32.32 regulate bone mass through MARK3 signaling in osteoblasts

Osteoporosis, characterized by decreased bone mineral density (BMD)and increased risk of fracture, continues to affect global health with an alarming financial and social burden. BMD is a highly heritable, complex trait and one of the single strongest predictors of fracture. The application of GWAS over a decade ago offered great promise as a means to better define the genetic determinants of BMD and susceptibility to osteoporosis. GWASs have determined that BMD is influenced by hundreds of polymorphic genomic elements, most of which affect the expression of skeletal regulatory genes. However, despite an enormous amount of effort and capital expended to date, there is still very little information connecting the responsible causal variants to bone-active genes. This knowledge gap is in part the result of inherent limitations of the GWAS approach. For example, a typical GWAS locus may contain tens to hundreds of variants associated with a phenotype due to the phenomenon of linkage disequilibrium. In addition, the vast majority of loci identified by GWAS are associated with noncoding, regulatory variants, which are more difficult to link mechanistically to target genes. Thus, to use GWAS data to inform bone biology and the development of new therapies, it is critical to define the mechanisms linking genetic variation to causal genes and to understand how causal genes are regulated and function in bone-active pathways.

To address these limitations, we developed a strategy to predict potentially causal genes using bone coexpression networks (1, 2). Using this approach, we identified microtubule affinity regulating kinase 3 (MARK3) as a candidate gene for an estimated bone mineral density (eBMD) locus located on chromosome (chr) 14 at 103.4 Mbp. MARK3 was further implicated as being causal based on strong expression quantitative trait locus (eQTL) regulating its gene-level expression that colocalized with the eBMD association. MARK3 is a member of the conserved PAR (partitioning-defective) gene family (PAR1–6) first identified in Caenorhabditis elegans as essential in asymmetric cell division and polarized cell growth (3, 4). Mammalian PAR-1 encodes a serine/threonine kinase initially designated microtubule affinity regulating kinase (MARK) based on its ability to phosphorylate microtubule associated proteins (MAPs), such as MAP2, MAP4, and Tau (5, 6). In multiple organisms and cell types, PAR1 activity depends on LKB1 (PAR4), which phosphorylates Par-1 and the related AMP activated kinase AMPK on conserved activation loop domains. This property could explain the ability of LKB1 to regulate both cell polarization and metabolic functions.

In this study, we demonstrate that SNPs associated with BMD alter the splicing and transcript levels of MARK3. Functional assessment in mice showed that global and osteoblast-specific loss of function of Mark3 results in increased bone mass. Mechanistic studies in mouse and human osteoblasts indicate that MARK3-generated signals impinge on the Notch pathway to regulate osteoblast activity. Together, these findings reveal a mechanism whereby genetically regulated alterations in Mark3 expression perturb cell signaling in osteoblasts to influence bone mass.

Exon-and intron-level expression of MARK3. We previously reported an eQTL for MARK3 transcript levels summed across all exons of the 26 predicted transcript isoforms, which colocalized with an eBMD association on chr 14. Nonreference alleles of the lead SNPs were associated with increased MARK3 transcript levels and decreased BMD. To determine whether the eQTL was regulating all or a subset of isoforms, we evaluated exon- and intron-level expression using Genotype–Tissue Expression Project (GTEx) (version 6) thyroid RNA-Seq data stratified by genotype at rs11623869, one of the SNPs strongly associated with BMD (P = 7.5 × 10^–110). We used data from GTEx for this project due to the paucity of publicly available eQTL data on bone or bone cells and the fact we have previously shown that the MARK3 colocalizing eQTL was present in a large number of GTEx tissues, suggesting that its genetic regulation is under similar regulation across many tissues, including bone and bone cells. There are a total of 26 MARK3 exons across isoforms. Of the 16 exons found in the canonical full-length isoform (MARK3-201; ENST00000216288.11), we observed that increased read counts were associated with the nonreference T allele of rs11623869 for all 16 (14 of 16 had a significant effect of genotype at an adjusted P value (P_adj) < 0.05) (Figure 1, A and B), confirming our prior observation using gene-level eQTL data. The levels of the other 8 exons did not differ (P > 0.3) by rs11623869 genotype. Surprisingly, when we looked at intronic read counts, the rs11623869 genotype had a strong (P = 1.6 × 10^–16) effect on intron 5 (Figure 1C). We observed a significant increase in the number of reads mapping to intron 5 in individuals carrying the nonreference T allele of rs11623869 (Figure 1D).

Figure 1

Genetically regulated alterations in MARK3 splicing colocalize with an eBMD GWAS locus on chr 14 at 103.4 Mbp. (A) The levels of all exons of the canonical MARK3 transcript are affected by the rs11623869 genotype. (B) An example of exon 7 shows increased normalized RNA-Seq read counts in individuals homozygous for the nonreference alleles (TT) of rs11623869. (C) Reads mapping to intron 5 are influenced by rs11623869 genotype. (D) Normalized intron 5 reads were increased in individuals homozygous for the nonreference alleles (TT) of rs11623869. (E) sQTL influencing the excision of intron 5 in GTEx (version 8) pancreas tissue colocalizes with an eBMD GWAS locus on chr 14 at 103.4 Mbp. (F) sQTL influencing the excision of intron 5 in GTEx (version 8) salivary gland tissue colocalizes with an eBMD GWAS locus on chr 14 at 103.4 Mbp. (G) sQTL influencing the excision of intron 5 in GTEx (version 8) testis tissue colocalizes with an eBMD GWAS locus on chr 14 at 103.4 Mbp. (H) Cell lines with at least 1 rs11623869 T allele have a larger number of reads mapping to intron 5. (I) Schematic of the exons formed in intron 5 of MARK3. PPH4, colocalization posterior probability of hypothesis 4, which states that a single causal variant drives both the sQTL and eBMD association.

Two distinct MARK3 splicing QTLs colocalize with an eBMD association on chr14. In the latest GTEx release (version 8) the consortium identified splicing QTLs (sQTLs) in addition to gene-level eQTLs. This analysis was performed using the leafcutter algorithm (7), which measures intron-excision ratios as a measure of alternative splicing. For MARK3, there were a total of 3 significant sQTL involving intron 5. Two of the sQTL were observed in pancreas (P = 4.4 × 10^–31) (Figure 1E) and salivary gland (P = 5.4 × 10^–9) (Figure 1F) tissues, and both represented the excision ratio of the entire intron 5. Individuals with the nonreference T allele at rs11623869 had decreased levels of intron 5 excision, suggesting that intron 5 was inefficiently spliced in those individuals. The other sQTL, found in testis (P = 2.5 × 10^–73) (Figure 1G), represented a novel intron (novel exon-exon junction) that began at the 5′ end of intron 5 and ended in the middle of intron 5 (hg38; chr14:103451983:103455319). All 3 sQTL strongly (PPH4 > 0.95) colocalized with the eBMD GWAS association.

To confirm these observations, we used nuclear RNA-Seq on cell lines from the ENCODE project (8). IMR90, a fibroblast cell line, was heterozygous for the top eBMD-associated SNPs, including rs11623869. In IMR90, we observed a large number of reads mapping to intron 5 (Figure 1H). Many were split reads spanning exons 5 and 6, supporting inefficient splicing of intron 5. Additionally, there was a relative absence of reads in intron 5 in NHEK cells, which are homozygous reference (GG) at rs11623869, and a large number of reads in HEPG2 cells, which are homozygous nonreference (TT) at rs11623869 (Figure 1H). A total of 134 SNPs were in linkage disequilibrium (r² > 0.5) with rs11623869 in individuals of European ancestry (Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/JCI142580DS1). Six of these SNPs are located in intron 5. In IMR90 nuclear RNA-Seq data from ENCODE (GEO GSE90262), all 6 showed significant allele-specific expression with more than 80% of the reads overlapping the 6 SNPs coming from the nonreference haplotype (Table 1). The human fetal osteoblast (hFOB) cell line, which was also heterozygous for rs11623869, demonstrated the same allele-specific expression pattern for the SNPs in intron 5 (Table 1).

Table 1

Allele-specific read counts for SNPs located in intron 5 in hFOB and IMR90 cells

In both IMR90 and hFOB, we also confirmed the formation of an exon in intron 5. In RNA-Seq data from both cell lines, we observed 111 reads supporting the presence of an exon that began at chr 14:103,455,317 and was preceded by a 3′ splice site AG (Figure 1I). We identified different ending sites for the exon. One ended at chr 14:103,455,353 and was followed by rs7150317 (chr 14:103,455,354), one of the 6 BMD associated intron 5 SNPS, which created a canonical 5′ splice site AG|GT on the nonreference allele (reference allele was AG|CT) (Figure 1I). The other version of the exon ended at chr 14:103,455,382 and was followed by a noncanonical 5′ slice site (AG|TTAATT) (ref. 9 and Figure 1I). Both forms of the exon contained a premature stop codon (PTC) in the middle of the MARK3 kinase domain, suggesting transcripts with the exon would not produce a functional protein and be targeted for nonsense-mediated decay (NMD).

The data regarding altered intron 5 splicing on the nonreference haplotype would seem to be at odds with the increased protein-coding transcripts, given that both the inefficient splicing and formation of exons led to the formation of PTCs. However, as shown above, exons of the full-length MARK3 isoform showed higher levels in GTEx thyroid tissue in individuals possessing the nonreference haplotype. To confirm this, we amplified full-length MARK3 cDNA using a primer pair starting in exon 1 and ending in the last exon. A total of 400 clones were sequenced. Of the 400, 180 (45%) represented full-length MARK3 transcripts arising from the reference haplotype (determined using the coding SNP rs13987) and 220 (55%) from the nonreference haplotype (1-tailed hypergeometric P = 0.025). Together, these data further support MARK3 as the causal gene at this locus. Additionally, they indicate that BMD-associated variants influence both MARK3 splicing and overall transcript levels.

Loss of MARK3 function in mice is associated with increased bone mass. The association between MARK3 expression and BMD in human tissue suggested by our previous work (1) was investigated in detail in genetically modified mice. Measurement of Mark3 mRNA expression levels across several mouse tissues (Figure 2A) showed high levels in mouse brain and kidney, as previously described (10, 11). Mark3 mRNA was also enriched in samples of cortical bone, growth plate, and bone marrow (Figure 2A). In primary mouse osteoblasts, Mark3 was expressed at levels 25-fold greater than that of 2 other family members, Mark1 and Mark2 (Figure 2B).

Figure 2

Unrestricted loss of Mark3 in mice results in reduced fat mass and higher bone acquisition. (A) Mark3 mRNA expression across mouse tissues (n = 4). (B) mRNA expression levels of Mark1, Mark2, and Mark3 in primary osteoblasts (n = 5). (C) Growth curve of Mark3 global KO (Mark3^–/–) mice and WT (Mark3^+/+) littermates from 4 weeks to 12 weeks of age (n = 6–10). (D) Body length and gonadal fat mass normalized by body weight at 12 weeks of age (n = 7–10). (E and F) Representative computer renderings and quantitative analysis of femoral trabecular (E) and cortical bone (F) by micro-CT in Mark3^+/+ and Mark3^–/– mice at 12 weeks of age (n = 7–8). BV/TV (%), bone volume/tissue volume; Tb. Th (mm), trabecular thickness; Ct. Th (mm), cortical thickness. (G) Histomorphometric analysis of Mark3^+/+ and Mark3^–/– mice at 12 weeks of age (n = 7–10). N. Ob/B.Pm, osteoblast numbers per bone perimeter; N. Oc/B. Pm, osteoclast numbers per bone perimeter; MAR, mineral apposition rate. (H) Serum P1NP and CTX levels (n = 7–9). Data are represented as mean ± SEM. *P < 0.05, Student’s t test between genotypes.

In order to evaluate the effect of Mark3 deletion on the bone, we first examined the phenotype of mice with unrestricted loss of the Mark3 gene. Immunoblots performed on the global KO mice indicated successful deletion of MARK3 protein in various tissues, including femoral cortical bone, bone marrow, white adipose tissue (WAT), hippocampus, and kidney (Supplemental Figure 1A). Global Mark3-knockout (Mark3^–/–) mice presented with significantly reduced body weight at all indicated ages (Figure 2C) compared with WT littermates (Mark3^+/+), which was attributed to decreased body length and abdominal fat depots (Figure 2D), and in agreement with the original description of Mark3^–/– mice (12). Morphometric analysis of the femoral bone by micro-CT showed significantly increased trabecular (Figure 2E) and cortical bone thickness (Figure 2F) in Mark3^–/– mice compared with WT littermates. Histomorphometry performed at the femoral metaphysis showed a significant increase in both osteoblast and osteoclast number per bone perimeter and increased bone-formation rate/bone surface (BFR/BS) in Mark3^–/– mice when compared with WT (Figure 2G). Mineralization surface/BS (MS/BS), which is a measure of the proportion of BS upon which new mineralized bone was being deposited during the period of calcein labeling, was also increased in Mark3^–/– mice (Figure 2G). Serum levels of the bone formation marker total procollagen type 1 N-terminal propeptide (P1NP) were significantly increased in Mark3^–/– mice compared with controls (Figure 2H), whereas carboxy-terminal collagen crosslinks (CTX) levels were similar between groups (Figure 2H). These results suggest that the increased bone volume following loss of Mark3 function results primarily from increased osteoblast number and activity.

To assess the impact of Mark3 deletion on bone mass without the potential confounding effects of alterations in body composition on skeletal acquisition, we generated mice lacking Mark3 selectively in osteoblasts by crossing mice expressing osteocalcin-Cre (Oc-Cre) with mice carrying Mark3 floxed alleles. Immunoblots performed in Mark3^fl/fl;Oc-Cre mice demonstrated a marked reduction in Mark3 protein expression specifically in bone samples, but not in other tissues (Supplemental Figure 1B). In contrast to the growth retardation observed in the global Mark3^–/– mice, Mark3^fl/fl;Oc-Cre mice had no discernible changes in overall growth, body weight, and length (Figure 3, A–C) at 12 weeks and 18 months of age. In addition, 12-week-old Mark3^fl/fl;Oc-Cre mice did not show alterations in gonadal fat mass (Figure 3D) as opposed to Mark3^–/– mice (Figure 2, C and D). Moreover, 12-week-old Mark3^fl/fl;Oc-Cre mice under a normal chow diet displayed unaltered insulin sensitivity and glucose tolerance (Figure 3E), which is consistent with the glucose metabolism of Mark3^–/– mice, as described previously (12). Similarly to the femoral phenotype in global Mark3^–/– mice, bone thickness in Mark3^fl/fl;Oc-Cre mice was increased in both trabecular and cortical bone (Figures 3, F and G), with concordant increases in the mean polar moment of inertia (MMI-polar), indicating the bone’s capability of resisting against torsional loads was increased in Mark3^fl/fl;Oc-Cre mice when compared with the control Mark3^fl/fl mice. The increased cortical thickness in the Mark3^fl/fl;Oc-Cre femoral bone resulted in increased resistance to fracture when assessed by 3-point bending (Figure 3, H and I). Representative stress-strain curves for Mark3^fl/fl;Oc-Cre and control mice group are shown in Figure 3H. The ultimate moment was shown here because it considers the span across which the bone was tested and is therefore more translatable across studies. Compared with control, 12-week-old Mark3^fl/fl;Oc-Cre mice exhibited an increase in the ultimate moment with a decrease in Young’s modules (Figure 3I). These changes were primarily driven by an increase in preyield strain leading to an increase in preyield energy (Figure 3I).

Figure 3

Mice lacking Mark3 in osteoblasts exhibit normal fat mass, but higher bone mass. (A) Growth curves for Mark3^fl/fl;Oc-Cre and control mice (Mark3^fl/fl) over 12 weeks (n = 5–6). (B) Body weight, (C) body length, and (D) gonadal fat mass of Mark3^fl/fl;Oc-Cre and Mark3^fl/fl mice at 12 weeks and 18 months of age (n = 5–6 for 12 weeks; n = 9–10 for 18 months). (E) Glucose tolerance test (GTT) and insulin tolerance test (ITT) in Mark3^fl/fl and Mark3^fl/fl;Oc-Cre mice at 12 weeks (n = 6–7). (F and G) Representative computer renderings of bone structure and quantitative analysis of the femoral trabecular (F) and cortical bone (G) generated from Mark3^fl/fl and Mark3^fl/fl;Oc-Cre mice at 12 weeks and 18 months old (n = 5–6 for 12 weeks; n = 9–10 for 18 months). (H) Representative stress-strain curves and stress-strain curves of Mark3^fl/fl;Oc-Cre and Mark3^fl/fl mice. (I) Ultimate moment, Young’s modulus, ultimate strain, ultimate stress, preyield energy, preyield strain, preyield stress and postyield energy as measured by 3-point bending test of the 12-week-old femur (n = 10). (J) Histomorphometric analysis of Mark3^fl/fl and Mark3^fl/fl;Oc-Cre mice at 12 weeks of age (n = 7–10). (K) Serum P1NP and CTX levels at 12 weeks (n = 10). Data are represented as mean ± SEM. *P < 0.05, Student’s t test between genotypes.

Histomorphometry performed in 12-week-old male Mark3^fl/fl;Oc-Cre and control mice showed no difference in either osteoblast and osteoclast number per bone perimeter (Figure 3J). However, there was a trend toward higher bone BFR/BS and MS/BS in Mark3^fl/fl;Oc-Cre mice compared with Mark3^fl/fl (Figure 3J), which may explain the increased bone thickness in Mark3^fl/fl;Oc-Cre mice compared with control. In agreement with increased bone thickness and similarly to Mark3^–/– mice, serum levels of P1NP were significantly increased (Figure 3K), while CTX levels were unchanged in Mark3^fl/fl;Oc-Cre mice compared with the control Mark3^fl/fl mice (Figure 3K).

MARK3 signals in osteoblasts impinge on the Notch pathway. To explore the mechanism of action of MARK3 in osteoblasts, we determined the effect of MARK3 loss of function on the differentiation of primary calvarial osteoblasts isolated from Mark3^fl/fl neonates. Adenoviral CRE-mediated knockout of Mark3 ( referred to herein as Ad-CRE) (Figure 4A) accelerated the onset of osteoblast differentiation, as indicted by increased Oc mRNA expression (Figure 4B) and increased alizarin red–stained mineral deposition at 14 days (Figure 4C) as compared with adenoviral-GFP transfected control primary osteoblasts (referred to herein as Ad-GFP). Analogous results were obtained in human osteoblasts following siRNA-mediated knockdown of Mark3 (Figure 4, D–F). BrdU incorporation (Figure 4G) and immunoblot of cleaved caspase 3 (Figure 4H) to assess osteoblast proliferation and apoptosis, respectively, revealed no differences between Ad-GFP– and Ad-Cre–transfected primary calvarial osteoblasts.

Figure 4

Mark3 deletion accelerates osteogenic differentiation in vitro via DVL/JNK/JAG1 pathway. (A) Mark3 and (B) Oc mRNA levels during osteogenic differentiation at days 0, 7, and 14 in Mark3^fl/fl primary osteoblasts transfected with Ad-GFP and Ad-CRE (n = 4). (C) Representative images and quantification of Alizarin red (ARS) staining after 14 days of osteogenic differentiation in Ad-GFP and Ad-Cre transfected primary osteoblasts (n = 3). (D) Mark3 mRNA levels after 48 hours of siRNA treatment in hBMSCs. (E) Representative images and quantification of ARS staining after 14 days of differentiation in control and siRNA transfected primary hBMSCs (n = 3). (F) qRT-PCR analysis of Oc expression after 14 days of osteoblast differentiation (n = 4) in primary hBMSCs. (G) BrdU incorporation for primary osteoblasts isolated from Mark3^fl/fl and Mark3^fl/fl;Oc-Cre mice (n = 5-6). (H) Western blot to detect cleaved caspase 3 in Ad-GFP– and Ad-Cre–transfected osteoblasts (n = 5). (I) Bulk RNA-Seq of Ad-GFP– and Ad-Cre–transfected Mark3^fl/fl primary osteoblasts identifies Notch signaling and cytoskeletal organization as potential pathways regulated by MARK3. (J) Relative mRNA expression of Jag1 and Hes1 at days 0, 7, and 14 of osteogenic differentiation Ad-GFP– and Ad-Cre–transfected Mark3^fl/fl osteoblasts (n = 4). (K) Western blot and quantification of JAG1 48 hours after adenovirus treatment of primary osteoblasts (n = 8). (L) Hes1 luciferase reporter assay on Ad-GFP– and Ad-Cre–transfected Mark3^fl/fl primary osteoblasts (n = 4–5). (M) Egr1 mRNA levels after 48 hours of Ad-GFP and Ad-Cre transfection (n = 5). (N) Egr1 and Jag1 mRNA after 48 hours of transfection of si-Egr1 in primary osteoblasts (n = 5). (O) Phos-tag SDS-PAGE of DVL to indicate phosphorylated forms of DVL after 48 hours in Ad-GFP– and Ad-Cre–transfected Mark3^fl/fl primary osteoblasts (n = 4). (P) Phos-tag SDS-PAGE of DVL to indicate phosphorylated forms of DVL after overexpression of the human MARK3 kinase domain in primary hBMSCs (n = 3). (Q) Immunoblot of p-JNK and t-JNK after 48 hours of Ad-GFP and Ad-Cre transfection (n = 6). Data are represented as mean ± SEM. *P < 0.05, Student’s t test between genotypes.

To identify possible downstream targets of MARK3 in osteoblasts, we profiled global RNA transcripts by RNA-Seq in Ad-GFP versus Ad-CRE Mark3^fl/fl primary osteoblasts after 24 hours of infection. Supplemental Table 3 shows all differentially expressed genes (P < 0.05). Pathway analyses of RNA-Seq conducted on control and Mark3-deficient osteoblasts revealed significant downregulation of genes involved in cytoskeleton organization, gap junctions, insulin signaling, apoptosis, and Notch signaling (Figure 4I). In particular, the expression of Jagged1 (Jag1) mRNA was downregulated, as were several transcription factors known to regulate Jag1 expression, including Neurl1a, Stat3, Tfcp2l1, Egr1, and Foxl (Figure 4I). All gene symbols, fold changes, and adjusted P values of the analyzed Notch pathway from RNA-Seq are included in Supplemental Table 2. In accordance with these predictions, Cre-mediated knockdown of Mark3 in primary osteoblasts reduced Jag1 and Hes1 mRNA expression (Figure 4J), decreased JAG1 protein levels (Figure 4K), and reduced the activity of a Hes1 reporter construct expressing luciferase (Figure 4L). As predicted from our pathway analysis (Figure 4I), the expression of Erg1, a known transcriptional regulator of Jag1 (13), was reduced in Mark3-efficient osteoblasts (Figure 4M), whereas siRNA-mediated knockdown of Egr1 reduced Jag1 mRNA (Figure 4N).

We next investigated the impact of Mark3 deficiency on downstream targets of Jag1. Mark3-deficient osteoblasts had decreased levels of phospho disheveled (DVL) (Figure 4O), a known downstream mediator of the MARK3 signaling cascade (14, 15). Conversely, overexpression of the human MARK3 kinase domain increased phosphorylation of DVL in primary human bone marrow stromal cells (hBMSCs) (Figure 4P). In addition, levels of phosphorylated JNK, an established target of DVL (16, 17), were significantly reduced in Mark3-deficient osteoblasts (Figure 4Q) whereas phosphorylated β-catenin (p–β-catenin) and total β-catenin (t–β-catenin) were not changed compared with that in controls (Supplemental Figure 2). These observations indicated a regulatory association between MARK3 and the Notch signaling pathway. The change of Egr1 is in accordance with the inhibited phosphorylation of DVL-JNK, as Egr1 is transcriptionally regulated by the JNK/c-Jun pathway (18, 19).

Osteoblast-targeted overexpression of Jag1 normalizes bone mass in Mark3^fl/fl;Oc-Cre mice. To further test the concept of MARK3 as a key regulator of Notch signaling in mature osteoblasts, we overexpressed Jag1 in osteoblasts derived from Mark3^+/+ and Mark3^–/– calvaria (Figure 5A) to rescue the effects of Mark3 deletion. Indeed, the enhanced alkaline phosphatase and mineral accumulation upon Mark3 deletion were reversed by overexpression of Jag1 plasmid in 14 day–differentiated osteoblasts (Figure 5B), which further suggests that JAG1 is a downstream mediator of MARK3 signaling during osteogenesis.

Figure 5

Targeted overexpression of Jag1 in osteoblasts reverses high bone mass in Mark3^fl/fl;Oc-Cre mice. (A) Efficiency of Jag1 overexpression after transfection measured by qRT-PCR (n = 3). (B) Left: ALP and ARS staining after 14 days of differentiation in Mark3^+/+ and Mark3^–/– primary osteoblasts transfected with a control or Jag1 overexpression plasmid. Right: quantification (n = 6). (C) Deletion of Mark3 and overexpression efficiency of Jag1 in femur of Mark3^fl/flJag1LSL^/+;Oc-Cre mice (n = 6–7). (D) Growth curves for Mark3^fl/fl, Jag1LSL, Jag1LSL^/+;Oc-Cre, Mark3^fl/flJag1LSL, Mark3^fl/fl;Oc-Cre, and Mark3^fl/flJag1LSL^/+;Oc-Cre over 12 weeks (n = 6–10). (E) Representative computer renderings of bone structure and quantitative analysis in the femoral trabecular and cortical bone in Mark3^fl/fl, Jag1LSL, Jag1LSL^/+;Oc-Cre, Mark3^fl/flJag1LSL, Mark3^fl/fl;Oc-Cre, and Mark3^fl/flJag1LSL^/+;Oc-Cre mice at 12w of age (n = 9-10). Data are represented as mean ± SEM. *P < 0.05 by Student’s t test between 2 groups or 1-way ANOVA followed by Student-Newman-Keuls test between multiple groups.

To gain additional support for NOTCH/JAG1 as mediators of MARK3 signaling in vivo, we analyzed bone mass acquisition in double-transgenic mice generated by crossing Mark3^fl/fl;Oc-Cre mice with osteoblast-specific heterozygous Jag1 overexpression (Supplemental Figure 3) (Mark3^fl/fl,Jag1LSL^/+;Oc-Cre). Quantitative real-time PCR (qRT-PCR) analysis showed effective deletion of Mark3 and overexpression of Jag1 in the respective mouse models (Figure 5C). Similarly to Mark3^fl/fl;Oc-Cre mice, both Jag1LSL^/+;Oc-Cre and the double-transgenic Mark^fl/fl;Jag1LSL^/+;Oc-Cre showed no significant changes in body weight and trabecular BV/TV (Figure 5, D and E). However, the increased trabecular and cortical thickness observed in Mark3^fl/fl;Oc-Cre mice was reversed in the double-transgenic Mark3^fl/fl;Jag1LSL^/+;Oc-Cre mice (Figure 5E), suggesting that JAG1 acts downstream of MARK3 and that its overexpression rescues Mark3 deletion.

MARK3 is a downstream target of LKB1 in osteoblasts. The results described above support a mechanistic link between the genomic regulatory elements that account for Mark3 expression and its downstream signaling, but do not provide information on how this kinase is regulated by upstream signals in osteoblasts. In this regard, previous studies have shown that the master kinase liver kinase B1 (LKB1) phosphorylates MARK3 in a variety of cell types (20). To determine whether this also occurs in osteoblasts, we compared the effect of overexpression of either WT Lkb1 or a mutant Lkb1-K78M construct deficient in kinase activity (21). Transfection efficiency was measured by RT-PCR (Figure 6, A and C). In line with our hypothesis, overexpression of Lkb1 increased the phosphorylation of MARK3 (Figure 6B), while the inactive kinase K78M-Lkb1 inhibited the phosphorylation of MARK3 (Figure 6D).

Figure 6

Phosphorylation of MARK3 Thr²³⁴ at by LKB1. (A–D) WT Lkb1 or the kinase dead (K78M) mutant Lkb1 plasmid were transfected into primary osteoblasts. Cells were collected 48 hours after transfection. (A and C) Real-time PCR with specially designed primers for WT and K78M Lkb1 to show transfection efficiency. (B and D) Immunoblot analysis of p-Thr²³⁴ MARK3 and total MARK3 after plasmid transfection. (E) Model of MARK3-regulated osteoblast differentiation.

A primary goal of our work is to apply systems-based approaches to linking BMD-associated variants to the expression of bone active genes. Using a network strategy, we previously identified MARK3 as a strong candidate gene for a BMD GWAS locus on chr 14q32.32 (1). BMD-associated variants in this region were found to influence the expression of MARK3 in several human tissues. In this study, we show that the previously identified gene-level MARK3 eQTL is due to multiple (inefficient splicing and the formation of 2 novel exons) splicing alterations, which are driven by BMD-associated variants on human chr 14 at 103.4 Mbp. We also confirmed that the BMD-associated reference haplotype in the region is associated with decreased levels of MARK3 protein-coding transcripts and increased BMD. Mechanistic studies in genetically altered mice were consistent with the inverse relationship between MARK3 and BMD and revealed that reduced Mark3 expression increases bone volume by impinging on the Notch pathway to stimulate osteoblast activity (Figure 6E).

Based on prior work, we initially predicted that the gene-level MARK3 eQTL would be due to differential transcription of one or more of the 26 MARK3 predicted transcript isoforms. However, an analysis of the effects of the rs11623869 genotype on the levels of each known MARK3 exon and intron demonstrated that all exons included in the canonical MARK3 transcript were increased in individuals carrying rs11623869 nonreference alleles. Surprisingly, the strongest effect of the rs11623869 genotype was on the number of reads mapping to intron 5. Individuals homozygous for reference alleles had few reads, whereas those homozygous for nonreference alleles had a relatively large number of reads mapping to intron 5. Using sQTL from GTEx and RNA-Seq data from IMR90 and hFOB cell lines, we confirmed altered intron 5 splicing. Our analysis identified 2 different intron-splicing defects due to genetic variation on the nonreference haplotype. The first was the inefficient splicing of intron 5, and the second was the formation of 2 novel exons. It is likely that inefficient splicing of intron 5 leads to the use of alternative splice sites, one of which is created by the nonreference allele of rs7150317 (one of the 6 eBMD-associated variants in intron 5). Importantly, both effects were driven by eQTL that strongly colocalized with the eBMD locus, suggesting that the intron 5–splicing defects are causal for the BMD association.

We predict that both splicing defects lead to a reduction in protein-coding transcripts from the nonreference haplotype. As we demonstrate, however, there is a clear increase in protein-coding messages from the nonreference haplotype. As a result, 1 of 2 scenarios have to be at play: either the intron 5–splicing defect is causing upregulation of MARK3 in cis or BMD-associated variants in the region are also increasing MARK3 expression in cis. Of the 134 eBMD-associated SNPs, 20 are located in putative regulatory elements and it is possible this is the mechanism whereby MARK3 expression is increased. Thus, while we have referred to rs11623869 as the “lead” SNP in the region, it is likely not a causal variant. Current studies are underway to systematically test each variant using CRISPR/Cas9 editing in HFOB cells.

The finding of increased bone volume in mice lacking Mark3 is consistent with the prediction that the eQTL at human chr 14 at 103.4 represses MARK3 function. Importantly, a side-by-side comparison of mice lacking Mark3 globally or selectively in osteoblasts revealed a similar bone phenotype characterized by increased bone volume that was most pronounced in the cortical compartment. This finding demonstrates that bone is an important site of MARK3 action and the genomic elements regulating MARK3 expression discussed above are manifest in osteoblasts rather than through indirect mechanisms (e.g., body composition or general metabolic alterations).

Our mechanistic studies were designed to link the genomic control of Mark3 expression to downstream targets regulating osteoblast number and/or activity. To identify Mark3 downstream targets we focused on the Notch signaling pathway for several reasons. First, Notch lig ands, including Jag1, were among the transcripts most downregulated in the RNA-Seq data set. Second, previous studies in kidney epithelial cells have shown that Mark3 affects Notch signaling to influence cell differentiation (11). Our results strongly suggest that JAG1 lies downstream of MARK3 in osteoblasts and that the increased bone mass in Mark3-deficient mice occur as the result of decreased Notch signaling. Thus, disruption of Mark3 suppressed Notch signaling, increasing the pace of osteoblast differentiation. Forced expression of Jag1 in Mark3-deficient osteoblasts reversed these effects. Analogous results showed reversal of the high cortical bone mass in compound mice overexpressing Jag1 in a Mark3-deficient background. However, we fully recognize the limitations of “overexpression” models and therefore cannot exclude the possibility that Jag/Notch signaling operating in a parallel (Mark3 independent) pathway might have contributed in part to reversal of both in vitro and in vivo phenotypes.

Our finding of increased bone volume following loss of osteoblast Mark3 and attenuated downstream Jag/Notch signaling is broadly compatible with previous findings in mice lacking Jag1 or Notch in osteolineage cells. For example, Youngstrom et al. (22) reported increased cortical thickness, but no change in trabecular bone volume in Prx1-Cre Jag1^fl/fl mice at 2 months of age. Lawal at al. (23), using the same Prx1-Cre to delete Jag1 in a different genetic background (FVB/N), reported higher trabecular bone mass, but no change in cortical bone thickness. Yorgan et al. reported that Notch2 inactivation in osteoblasts using Runx2-Cre increased trabecular bone in the proximal femur and distal tibia at 26 weeks, but there was no change in cortical bone thickness at both 12 weeks and 26 weeks of age (24). In another study (25), mice lacking both Notch1 and Notch2 in Prx1-Cre–expressing cells had increased trabecular bone volume at 8 weeks, but then demonstrated a marked loss to 10% of control with age (26 weeks). Close inspection of our micro-CT data in 18-month-old Mark3^fl/fl;Oc-Cre mice suggests that trabecular bone volume was indeed trending down (Figure 3G), even though cortical bone remained elevated. It appears therefore that attenuation of Jag1/Notch signaling can affect both trabecular and cortical bone and that this effect is strongly dependent on time of deletion in the osteogenic lineage, genetic background, and age of the mouse.

Although the current studies were primarily designed to establish the mechanisms responsible for genomic control of Mark3 activity, our results implicate this kinase as an important signaling molecule in osteoblasts. In this regard, we felt it was important to also establish linkage of Mark3 to upstream signaling pathways that might control osteoblast activity in response to local or systemic queues during development or in response to hormonal stimuli. Our finding that LKB1 regulates MARK3 activity in osteoblasts accords with previous findings in other cell types, which have established the proximate position of this kinase in controlling fundamental processes, including cell-cycle progression, cellular energy homeostasis, and cell polarity (14, 20, 26). Indeed, previous studies have shown that mice lacking LKB1 (aka Stk11) in osteoblast precursors also developed increased bone mass (27). In addition, a recent study by Nashimori (28) showed that the MARK3-related salt-inducible kinases SIK2 and SIK3 (SIKs) mediate skeletal actions downstream of PTH1R and are inhibited by PTH. Moreover, mice lacking Sik2/3 in osteoblasts have increased bone mass with morphometric features similar to those of our mice that lack Mark3. Current studies are being conducted to determine whether PTH signaling alters MARK3 activity and to define potential overlap between these 2 related signaling components in osteoblasts.

Mouse and treatment protocols. Global Mark3-KO mice were created by crossing Mark3^fl/fl mice (Mark3^{tm1a(KOMP)Mbp}, https://www.mousephenotype.org/data/alleles/MGI:1341865/tm1a(KOMP)Mbp) to mice expressing a ubiquitous Cre deletion strain (B6.C-Tg (CMV-Cre)1Cgn/J; JAX, 006054) to remove the Mark3 gene globally. Mark3-KO mice were backcrossed to WT C57BL/6 mice to remove Cre and bred to homozygosity. Mice lacking Mark3 in osteoblasts (Mark3^fl/fl;Oc-Cre) were generated by crossing Oc-Cre transgenic mice (29) with Mark3^fl/fl mice. R26-LSL-JAG1 mice, which express JAG1 after exposure to Cre recombinase, were purchased from The Jackson Laboratory (catalog 030173). For this paper, we used Jag1LSL mice as controls, which were homozygous for the transgene but lacked Cre-expression, while Jag1LSL^/+;Oc-Cre mice were heterozygous for the transgene and served as a model of osteoblast-specific overexpression of Jag1. Results are reported for male mice on the same genetic background (C57BL6/J). Genotyping primers for Mark3 global KO and Mark3^fl/fl mice are available in the Supplemental Methods. Genotyping strategies for R26-LSL-JAG1 mice can be found on the Jackson Laboratory website (catalog 030173). Serum bone turnover markers P1NP (IDS, catalog AC-33F1) and CTX (IDS, catalog AC-06F1) were measured by ELISA.

Evaluating exon– and intron-level expression of MARK3. Raw RNA-Seq data for thyroid tissue were downloaded from GTEx for 188 individuals. BAM files specific for the MARK3 genomic region were created using HISAT2 (30). Read counts for each MARK3 exon and intron were generated using BEDTOOLS (31). The effect of rs11623869 was tested using a linear model that included the effects of subject age and RNA integrity number.

Characterizing MARK3 sQTLs. sQTLs for MARK3 were downloaded from the GTEx portal, and eBMD GWAS data were downloaded from the GEnetic Factors for OSteoporosis Consortium website (http://www.gefos.org/?q=content/data-release-2018). Colocalization analyses were performed using the coloc R package (32). Investigation of MARK3 intron 5 allele–specific expression and the formation of exons utilized data from the following cell lines: IMR90, hFOB, HEPG2, and NHEK. Aligned RNA-Seq data were used for HEPG2 and NHEK and were downloaded from the Encyclopedia of DNA Elements at UCSC (https://genome.ucsc.edu/ENCODE/). Fastq files for IMR90 generated through the ENCODE Project (8) were downloaded from the Encyclopedia of DNA Elements (https://www.encodeproject.org/). RNA used for RNA-Seq was isolated from the human osteoblast immortalized cell line hFOB 1.19 (ATCC) with the mirVANA RNA Isolation Kit. For hFOB RNA, total RNA-Seq libraries were constructed using Illumina TruSeq Stranded Total RNA HT sample prep kits and sequenced at Psomagen. The quality of IMR90 and hFOB RNA-Seq fastq files were verified using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and were aligned to a SNP aware version of the human reference genome (GRCh38) using HISAT2, version 2.0.5 (30).

MARK3 cDNA cloning. All reagents listed below were from Invitrogen and used in accordance with the manufacturer’s recommendations unless otherwise noted. RNA used in construction of the libraries was isolated from the human osteoblast immortalized cell line hFOB 1.19 (ATCC) with the mirVANA RNA Isolation Iit. cDNA was synthesized with ProtoScript II reverse transcriptase (NEB) and amplified in a PCR with platinum Taq polymerase and a forward primer approximately 50 bp 5′ of the start codon (30 of the 31 hMARK3 reference sequence transcript variants; E1F3: GCTGTGCTGTTTTGTTTTGACC) and a reverse primer approximately 30 bp 5′ of the poly A addition site (all of the 31 hMARK3 reference sequence transcript variants; E16R2: GAATCTACTTAGAAAAACACACACAAGC). Products from 8 separate 25 μl reactions were combined, “cleaned up” with the Pure Link Quick PCR Purification Kit, and cloned into plasmid PCR 2.1 Topo; this was used to transform One Shot Top10 chemically competent E. coli cells. Colonies that arose after 16 to 20 hours of incubation at 37°C on LB/kanamycin agar plates were used to inoculate 5 ml of LB/kanamycin media and plasmid DNA from 16- to 20-hour cultures was isolated with a Pure Link Quick Plasmid Miniprep Kit. Plasmids that contained inserts (2–3 kb) were sequenced from plasmid primers located on either side of the cloning insert site (M13 reverse and T7) or primers E1F3 and E16R2 as well as one or more primers whose sequences were homologous to the MARK3 coding region. This was done in order to obtain complete sequences for all inserts reported in this study. The protocol outlined above was repeated 2 additional times (3 total) from separate cell cultures to obtain the results reported here with at least 100 clones (plasmid DNAs with inserts) coming from each preparation.

Primary cell isolation, culture, and treatment. Mouse osteoblasts were isolated from calvaria of newborn mice by serial digestion in 1.8 mg/ml of collagenase. For flow cytometry and BrdU analysis, primary osteoblasts were isolated from Mark3^fl/fl and Mark3^fl/fl;Oc-Cre mice alternatively. For other in vitro analysis, primary osteoblasts isolated from Mark3^fl/fl mice were infected with adenovirus encoding Cre recombinase or green fluorescent protein (Vector Biolabs), as previously described (33). Infection with 100 MOI was used in all experiments. Human primary bone marrow mesenchymal cells (BMSCs) were a gift from Aaron James (Johns Hopkins University). Primary BMSCs of deidentified arthroplasty specimens of the human femur and tibia were flushed with PBS and then cultured. siRNA was purchased from Ambion, QIAGEN, and Integrated DNA Technologies Co. and transfected with Lipofectamine RNAiMAX Transfection Reagent (Thermo Fisher Scientific) after 3 days of osteogenic differentiation of BMSCs. For the forced human Mark3 kinase domain expression, MARK3_HUMAN_D0 plasmid (catalog 79740 Addgene) (34) and the control plasmid pET His10 TEV LIC cloning vector were used. For Lkb1 WT and mutant kinase overexpression, GFP-Lkb1 (catalog 21147, Addgene) and GFP-Lkb1-K78M (catalog 21148, Addgene), which encoded by an invariant nucleotide binding site mutant to produce a mutated LKB1 protein with deficiency in kinase activity, were used (21). The plasmids were transfected using Lipofectamine 3000 Transfection Reagent (L3000008, Thermo Fisher Scientific) following the manufacturer’s protocol. Osteoblast differentiation was induced by supplementing aMEM media with 10% FBS, 10 mM β-glycerol phosphate, and 50 μg/ml ascorbic acid. Alkaline phosphatase and Alizarin Red S staining were carried out according to standard techniques.

RNA preparation and RNA-Seq analysis. Total RNA was extracted with the TRIzol method, reverse transcribed using the iScript cDNA Synthesis Kit (Bio-Rad), and amplified by real-time PCR using SYBR GREEN PCR Master Mix (Bio-Rad). Primer sequences were obtained from PrimerBank (http://pga.mgh.harvard.edu/primerbank/index.html). RNA-Seq libraries were constructed using TruSeq RNA Library Prep Kit, version 2, sample prep kits (Illumina). Samples were sequenced (2 × 75 bp paired-end reads) on an Illumina NextSeq500 sequencer. Quality of the sequencing reads was verified using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Fastq files were aligned to the mouse reference (GRCm38) using HISAT2, version 2.0.5 (30). Raw read counts were generated using Stringtie (35) and the prepDE.py python script (http://ccb.jhu.edu/software/stringtie/dl/prepDE.py). Differential expression analysis was performed using DESeq2 (36). Pathway analysis was performed using DAVID Bioinformatics Resources (37).

Immunoblot analysis. To extract protein from long bones, the attached muscles and tissues were removed, the epiphyseal ends were cut, and the marrow was flushed out; then the bone was ground in liquid nitrogen in the presence of proteinase inhibitor and phosphatase inhibitor. Immediately after grinding, the ground femur was transferred to a microfuge tube containing lysis buffer. For cell samples, total protein was extracted using RIPA lysis buffer. Immunoblot was performed according to standard techniques. Primary antibodies were as follows: anti-MARK3 (1:1000, catalog 9311, Cell Signaling Technology), anti-p-MARKs (1:1000, catalog 4836, Cell Signaling Technology), anti-p-JNK (1:1000, catalog 4668, Cell Signaling Technology), anti-t-JNK (1:1000, catalog 9252, Cell Signaling Technology), anti-JAG1 (1:200, catalog sc-390177, Santa Cruz Biotechnology Inc.), anti-DVL (1:200, catalog sc-8025, Santa Cruz Biotechnology Inc.), p–β-catenin (Ser33/37/Thr41) (1:1000, catalog 9561, Cell Signaling Technology), β-catenin (1:1000, catalog 9562, Cell Signaling Technology), and anti-GAPDH (1:1000, 10494-I-AP, Proteintech). Antibodies were incubated overnight at 4°C.

Phos-tag SDS-PAGE and Western blotting. Cells were collected with ice-cold lysis buffer containing 50 mM Tris/HCl, pH 7.4, 1% Triton X-100, 10% glycerol, 1 mM sodium orthovanadate, 50 mM NaF, 10 mM 2-glycerophosphate, 5 mM sodium pyrophosphate, 0.1 mg/ml microcystin-LR, and EDTA-free protease inhibitor cocktail (Roche). Lysates were centrifuged at 20,800 g and then mixed with 4× SDS/PAGE sample buffer (250 mM Tris/HCl, pH 6.8, 8% SDS, 40% glycerol, 0.02% bromophenol blue, and 4% 2-mercaptoethanol) and heated at 95°C for 5 minutes. For phos-tag SDS-PAGE, samples were mixed with 10 mM MnCl₂ before loading gels. Phos-tag SDS-PAGE was performed following the protocol described previously, with slight modifications (38). Gels for Phos-tag SDS/PAGE were composed of a stacking gel (4% acrylamide, 125 mM Tris/HCl, pH 6.8, 0.1% SDS, 0.08% ammonium persulfate [APS] and 0.2% N,N,N0,N0-tetramethylethylenediamine [TEMED]) and a separating gel (12% acrylamide, 375 mM Tris/HCl, pH 8.8, 0.1% SDS, 0.05% APS, 100 mM Phos-tag acrylamide, 200 mM MnCl₂, and 0.1% TEMED). 30 mg of protein sample were run at 70V for the stacking gel and at 130 V for the separating gel with the running buffer (25 mM Tris/HCl, 192 mM glycine, and 0.1% SDS). After SDS-PAGE, gels were washed 3 times with transfer buffer (48 mM Tris/HCl, 39 mM glycine, and 20% methanol) containing 10 mM EDTA and 0.05% SDS, followed by a 10-minute wash with transfer buffer containing 0.05% SDS. Blotting transferring to PVDF membranes was carried out at 90 V for 180 minutes on ice in the transfer buffer without SDS/EDTA. Transferred membranes were blocked with 5% nonfat dry milk dissolved in TBS-T (20 mM Tris/HCl, pH 7.5, 150 mM NaCl, and 0.1% Tween 20) at room temperature for 45 minutes, then incubated with primary antibodies overnight at 4°C. After washing in TBST, membranes were incubated at room temperature for 1 hour with secondary antibody. After washing in TBST, protein bands were detected using an ECL solution (Amersham ECL Western Blotting Detection Reagents; GE Healthcare).

Flow cytometry. Osteoblast proliferation was measured by flow cytometry with a BrdU kit (559619, BD Biosciences — Pharmingen). Briefly, control and Mark3^–/– osteoblasts were cultured in α-MEM containing 10% FBS. At 30% confluence, BrdU (10 μM) was added to the medium for 16 hours before collecting the cells. Cells were stained with anti-BrdU-FITC and 7-amino-actinomycin D (7AAD) and assessed by LSR II (BD Biosciences). Ten thousand events were collected for every sample, and the results were analyzed by FlowJo, version 7.6.

Luciferase assay. Primary osteoblast cells from Mark3^fl/fl and Mark3^fl/fl;Oc-Cre KO mice were transfected with pHes1-luciferase reporter gene (catalog 43806, Addgene) and pGL2- Empty control-Luc (39). pCMV-RL (Promega) was cotransfected using Lipofectamine 3000 Transfection Reagent (L3000008, Thermo Fisher Scientific). Twenty-four hours later, lysate luminescence was measured using the Dual-Luciferase Reporter Assay System (E1910, Promega).

Imaging and histomorphometry. Mice were sacrificed at the indicated age, and bone volume in the distal femoral metaphysis was assessed using a desktop microtomographic imaging system (Skyscan 1275). For the trabecular bone analysis, a region of interest (ROI) was selected and delineated into 2 mm segments along the longitudinal direction. The vertical distance between the first image of the ROI and the growth plate was 0.5 mm. For the cortical bone analysis, the ROI was selected from the midpoint of the femur, 0.25 mm toward the proximal end and 0.25 mm toward the distal end. Histological analyses, using a semiautomatic method (Osteoplan II, Kontron Instruments Ltd.) were carried out on 12-week-old mice injected with 1% calcein (w/v) at 8 and 3 days before sacrifice. Histomorphometric parameters followed the recommended nomenclature of the American Society of Bone and Mineral Research (40). Parameters were measured at a standardized site below the growth plate using a semiautomatic method (1 optical field, 180 microns at ×20 magnification). This ROI covered the secondary spongiosa at approximately 50 microns from each cortice and was 0.8 mm × 1 mm in dimension.

Mechanical testing. The strength of the femoral midshaft was assessed by 3-point bending using a low force mechanical testing system (Bose Electroforce 3100) according to established methods (41–43). From force-displacement data and micro-CT imaging, structural and material properties were calculated using a custom MATLAB (MathWorks) program, as previously reported (44). Left femur was dissected, cleaned of soft tissue, wrapped in gauze soaked in PBS, and stored at –80°C. Before mechanical testing, femurs were thawed at 4°C overnight. On the morning of testing, all samples were quick scanned using a Skyscan 1275 micro-CT system, in order to obtain the cross-sectional geometry at the failure site and calculate the material properties from tests of femur. All the samples were allowed to achieve room temperature prior to mechanical testing. The femurs were loaded to failure by 3-point bending and placed on the miniature materials testing machine as previously described (44) with a span length of 7 mm. Preload was applied to the bones before testing to ensure that bones were secure on the loading span. The crosshead speed during testing was 0.1 mm/s, and force-displacement data were collected every 0.05 seconds. Young’s modulus was calculated from data in the linear region of the stress-strain curve.

Data access. Raw data and processed data of RNA-Seq in the study have been submitted to the NCBI’s Gene Expression Omnibus (GEO GSE160415).

Statistics. All data are expressed as mean ± SEM. Significant differences were assessed using unpaired, 2-tailed Student’s t test or 1-way ANOVA followed by the Student-Newman-Keuls test. P < 0.05 was considered statistically significant.

Study approval. All mouse experiments and hBMSCs were approved by the IACUC of Johns Hopkins University in accordance with institutional guidelines. All human samples were obtained under IRB approval at Johns Hopkins University with a waiver of informed consent.