C57BL/6J and C57BL/6-Tg (CAG-EGFP)1Osb/J male mice were obtained from The Jackson Laboratory. Bgn-deficient (KO) mice were generated as previously described (Xu et al., 1998) all mice were bred and housed at the NIDCR/NIH/DHHS animal facility with standard conditions and genotyped in our laboratory. All procedures using mice were approved by the NIDCR/NIH/DHHS ACUC (protocol numbers 18-865, 18-871).

After euthanasia, femurs and tibias were dissected from 6 to 7 week-old WT or Bgn KO male mice. Muscle and connective tissue were removed under sterile conditions and the epiphyses were coated with 5% low melting point agarose (Invitrogen) to protect the cartilage from digestion. PDCs were isolated by 1 h digestion in 3 mg/mL collagenase type II (Gibco) and 4 mg/mL dispase (Gibco) in α-minimal essential medium (αMEM) with 3% glutaMAX (Gibco) at 37°C (van Gastel et al., 2012; van Gastel et al., 2014). The cell suspension was passed through a 70 μm cell strainer (Falcon) and plated in 25 cm² flasks with growth medium composed of αMEM (Gibco) supplemented with 20% FBS (Gemini Bio Products), 100 U/mL each of penicillin and streptomycin and 2 mM glutaMAX (Gibco), and the medium was replenished twice a week until cells reached confluence.

When PDC cultures reached 80%–90% confluency (10–14 days post isolation), the cells were trypsinized and seeded in RAFT™ 3D 96 well cell culture plates according to the manufacturer’s instructions (RAFT™ 3D Cell Culture bundle Kit, Lonza) at a concentration of 120,000 cells/well. Briefly, cells were seeded in a chilled mixture of collagen solution containing x10 MEM medium, 2 mg/mL rat tail collagen type I and neutralizing solution and incubated for 18 min at 37°C to form a hydrogel. Next, an absorber device was placed on top of the hydrogel for 15 min to condense and concentrate the mix of cells and collagen to create the RAFT™ disc. 0.24 mL/well growth medium was immediately added and thereafter replaced twice a week. To induce osteogenic differentiation, 10^–8 M dexamethasone, 100 μM L-ascorbic acid phosphate and 2 mM β-glycerophosphate were added to the growth medium containing only 10% FBS for 14 days.

For all surgeries, 6–7 week-old WT or Bgn KO male mice were anesthetized with 2%–3% isoflurane. 1 mg/kg of buprenorphine SR-LAB (ZooPharm) was administrated subcutaneously before the procedure and 72 h post-surgery. For open stabilized fractures, the right hindleg was shaved and scrubbed with 10% povidone-iodine and alcohol solutions. A 5 mm medial parapatellar incision was created and the patella was dislocated to expose the femoral condyle. A hole was burred into the femoral epiphysis using a 25 gauge (25 G) needle and enlarged using a 23 G needle. A 26 G cannula (Millpledge Veterinary) was inserted into the femoral shaft through the burred hole. The muscles and tendons were shifted to expose the bone shaft and a horizontal fracture was created at the center of the shaft using surgical scissors. The soft tissue was repositioned and the protruding end of the needle was cut off. The incision was closed using 5-0 absorbable sutures (Ethicon). X-ray was performed to visualize the fracture (Faxitron® Ultra Focus). For subcutaneous transplantation, a small incision was made in the skin on the back of the mouse. 3D scaffolds with un-differentiated PDCs (that served as control) or differentiated PDCs were implanted subcutaneously and the incision was closed using 5-0 absorbable sutures (Ethicon). The 3D scaffolds were analyzed at time points 0 (the day of implantation), and 4 and 8 weeks after implantation.

The 3D scaffolds recovered after implantation and mice femurs before and after fracture were fixed in Z-fix (170; Anatech, LTD) for 24 h at room temperature (RT) and stored in 70% ethanol. 3D transplants were scanned at 45 kV, 200 μA, 600 ms integration time (IT), 14.8 µm voxel size. Femurs were scanned at 70 kV, 85 μA, 300 ms IT, 10 µm resolution (µCT50; Scanco Medical AG, Brüttisellen, Switzerland). Mineralized tissues were reconstructed using the global approach and segmented by a global thresholding software (Scanco Medical AG, Brüttisellen, Switzerland). Standardized nomenclature was used for the bone parameters measured. For the femurs, trabecular parameters were measured at the secondary spongiosa of the distal metaphysis and cortical parameters were determined in a 1 mm ring at the mid-diaphyseal region according to previously published guidelines (Bouxsein et al., 2010). The callus of the fractured femurs was analyzed along the entire shaft of the bone. Cross section measurements of the callus were calculated based on the entire stack of 2D images.

For scanning electron microscopy (SEM), 3D scaffolds were fixed at 4°C overnight with Karnovsky fixative (Electron Microscopy Sciences) followed by three washes with phosphate buffered saline (PBS) and dehydrated in a series of ethanol washes (0, 30, 50, 70, 85, 95, and 100%), with each wash conducted 3 times for a total of 15 min. The samples were then critical point dried (Autosamdri-814, Tousimis), mounted and sputter-coated with gold (75 mA for 60 s, Denton Vacuum Desk II) prior to imaging (2 kV, 10 Ua, S-4700-II FE-SEM, Hitachi). For histology, immunofluorescent and second harmonic generation (SHG) imaging, 3D scaffolds or femurs were decalcified prior to paraffin-embedding.

Femur sections were deparaffinized and rehydrated prior to SHG imaging of type I collagen. SHG imaging was performed on a Nikon A1R MP + HD confocal system using a ×40 Apo LWD objective (N.A. 1.15) in resonant scanning mode (512 × 512, 4X line average scanning). The two-photon beam excitation was provided by a Chameleon Vision II laser tuned to 820 nm (Coherent). The non-descanned detector used a 400–480 band pass emission filter. Z-stacks were taken of the middle 10 or 15 microns of 20 or 30 micron-thick tissue slices, respectively. Image tiles were created in a ∼250–500 µm wide by 1800 µm long region centered around the femur center. Femur sections were oriented to match the two-photon beams polarity (parallel with the Y-axis). After acquisition, images were denoised using NIS Elements Denoise A.I. software (Nixon). Line scans depicted in figures show SHG signal intensity normalized to the brightest pixels for comparison on a similar scale.

To analyze SHG fibril differences between WT and KO femurs, ctFIRE and curve align standalone plugins for MatLab (Bredfeldt et al., 2014) were used with the following settings: Min fiber length 100 pixels, Max fiber width 15.

For Bgn immunochemistry, deparaffinized and rehydrated sections were incubated for 1 h at 37°C with ABCase (Seikagaku biobusiness corp.), followed by antigen retrieval (Unitrieve, Innovex), and quenching of endogenous peroxidase activity with dual endogenous enzyme block (Dako). The sections were blocked with 10% normal goat serum in PBS for 1 h at 37°C and incubated with Bgn rabbit antisera (from Dr. Larry W. Fisher, NIH, ref. LF-159) 1:600 diluted in the blocking buffer overnight at 4°C. The sections were washed and incubated with Super PicTure™ Polymer detection kit (Invitrogen) for 20 min at RT and detected with ImmPACT™ AEC (Vector laboratories). Slides were scanned using an Aperio ScanScope slide scanner.

For F4/80, Ym1, and iNOS immunohistochemistry, deparaffinized and rehydrated sections were incubated for 48 h at 55°C with citrate buffer, pH 6.0 (BIOZOL, Germany). The sections were blocked with 3% H₂O₂ in PBS for 10 min and with Protein Block Serum-free solution (Dako, X0909) for 20 min. This was followed by incubation with rat anti-mouse F4/80 antibody (1:250, Biorad, MCA497G), rabbit anti-Ym1 antibody (1:100, Abcam, ab93034) or rabbit anti-iNOS (1:500, Enzo Life Sciences, ADI-905-431–1) diluted in antibody diluent (Dako S3022) overnight at 4°C. The sections were washed and incubated with Histofine Simple Stain Mouse MAX PO (Rat) (Nichirei Biosciences INC, 414311F) for F4/80, or with Histofine Simple Stain Mouse MAX PO (Rabbit) (Nichirei Biosciences INC, 414341F) for Ym1 and iNOS, for 30 min at RT and detected with Vector DAB Substrate (Vector laboratories, SK4100). Mayer’s hematoxylin was used as counterstain. The protein expression levels were determined by the reciprocal intensity of the chromogen stain using the open source ImageJ Fiji software (http://fiji.sc/Fiji).

For aggrecan (ACAN) and osteopontin (OPN) immunofluorescence, deparaffinized and rehydrated sections were incubated for 1 h at 37°C with ABCase (Seikagaku biobusiness corp). Following antigen retrieval (Unitrieve, Innovex), and incubation with 0.1% sodium borohydride, the sections were blocked with 10% normal goat serum in PBS for 1 h at 37°C, and incubated 1:200 with either rabbit anti-ACAN antibody (Millipore, ref. AB1031) or rabbit anti-OPN antibody (IBL, ref. 18,621) overnight at 4°C. The sections were washed and incubated with secondary antibody, anti-Rabbit IgG conjugated to Rhodamine Red 1:500 (Jackson, ref. AB_2340614) and DAPI 1:500 (Thermo Fisher Scientific, ref. D21490) for 1 h at RT. Finally, the slides were washed, mounted and imaged by A1R-MP + HD multiphoton confocal system (Nikon). Tiled (2 × 2, 6 × 6 or 7 × 7) Z-stacks (every 1.0 µm) were taken using a ×40 Plan Fluor (N.A 1.3) oil immersion objective in resonant scanning mode (×512512, 4X line average scanning) using the 405 nm and 561 nm lasers. Images were denoised using NIS Elements Denoise A.I. software after acquisition (Nixon). For quantification, we applied a Li threshold and then captured the mean fluorescence intensity of included area of the image using ImageJ Fiji software.

For osteoclast analysis, deparaffinized sections of fractured bones were stained with TRAP (Sigma).

Staining with hematoxylin and eosin (H&E) was performed using standard protocols. For periosteal width analysis, H&E stained slides were scanned using an Aperio ScanScope slide scanner and the entire enlarged periosteum was outlined, measured and the expanded area analyzed by Image Pro 7.0 software (Media Cybernetics Inc., United States). This resulted in five data points for each genotype used.

Femoral shafts were isolated and immediately frozen in liquid nitrogen. The bones were centered in a tissue tube (Covaris) under liquid nitrogen and crushed using a CP02 cryoPREP Dry Pulverizer. Total mRNA was extracted and purified from the pulverized tissue using TriPure (Sigma, United States) followed by RNeasy mini kit (Qiagen) hybrid protocol.

Total mRNA from cultured cells was extracted and purified using RNeasy mini kit (Qiagen).

Total mRNA was converted to cDNA using iScript cDNA Synthesis Kit (Bio-Rad) and qRT-PCR analysis was performed using iQ SYBR Green Supermix (Bio-Rad). Target gene (Bgn: F: 5′-AGA​CAA​ACC​GAC​AGC​CTG​ACA​AC-3′, R: 5′-GCC​AGC​AGC​AAG​GTG​AGT​AGC-3′) was normalized to S29 (F: 5′-GGA​GTC​ACC​CAC​GGA​GTT​CG-3′, R: 5′-GGA​AGC​AGC​TGG​CGG​CAC​ATG-3′) and relative expression data was calculated using the ΔΔCt method.

For RNA seq procedures, RNA was transcribed by Superscript IV (Thermo Fisher Scientific) and full-length 2nd strand cDNA amplified by LongAmp DNA polymerase (New England BioLabs). Sequencing libraries were prepared using a Nextera XT kit (Illumina), individually barcoded, pooled to a 2 nM final concentration, and sequenced on a NextSeq500 or NextSeq2000 instrument (Illumina) using 37 × 37 paired-end (NextSeq 500) or 55 × 55 (NextSeq 2000) paired-end read configurations. After sequencing, the base-called demultiplexed (fastq) reads from multiple sequencing runs were merged when appropriate and read qualities were determined using FastQC (v0.11.2), aligned to the GENCODE M11 mouse genome (GRCm38.p4) and gene counts generated using STAR (v2.5.2a). Post-alignment qualities were generated with Picard Tools RnaSeqMetrics (v1.129). An expression matrix of raw gene counts was generated using R and filtered to remove low count genes (defined as those with less than five reads in at least one sample). The filtered expression matrix was used to generate a list of differentially expressed genes (DEGs) between the sample groups using three statistical methods: DESeq2, EdgeR, and Limma-voom.

For single-cell RNA sequencing (scRNA seq), WT and Bgn KO embryos (E16.5) were used. After euthanasia, the posterior limbs were dissected and connective tissues were removed under sterile conditions. The bones were dissected into small pieces and digested with 2 mg/mL collagenase IV (Gibco) in Advanced DMEM (supplement with 1% Glutamine) for 40 min at 37°C with manual shaking every 10 min. When the dissociation procedure was complete, the cell-suspension was centrifuged, supernatants were carefully removed and the cell pellet suspended in PBS containing 0.1% BSA. The cells were then filtered using Flowmi filter tips (Bel-Art Products, #H13680-0040) to remove possible undissociated cells and debris. The cells were counted and resuspended to a final concentration of 1 × 10⁶ cells/mL.

The single-cell suspension was loaded onto the droplet-based single-cell barcoding system (10x Chromium Controller, 10x Genomics) and a Chromium NextGEM Single Cell 3’ Reagent Kit v3.1 (10x Genomics) was used to prepare single-cell, barcoded 3’ cDNA libraries according to the manufacturer’’s instructions. The libraries were sequenced on a NextSeq2000 instrument (Illumina). The sequenced data was processed by CellRanger (v5.0.0, filtering, barcode and UMI counting) using default command line options and gene-barcode matrices were generated. The sequenced reads were aligned to the mouse genome assembly provided by 10X Genomics (10X Genomics reference mm10-2020-A), based on Ensembl annotation. Downstream data analysis was performed on the CellRanger cell-gene filtered matrix, using the Seurat R package V4 (Hao et al., 2021). Seurat objects were generated for each of the four samples (2 WT, two Bgn KO), which were then merged into a single Seurat object. Cells with unusual numbers of genes (<200), UMI count (>15,000) or percentage of mitochondrial genes (>20%) were filtered out. The data was then normalized by the “LogNormalization” method implemented in the Seurat package (scale factor = 10,000) and scaled using Seurat’s default settings. The 2,000 top variable genes were identified using the “vst” method implemented in Seurat. Linear dimensional reduction (PCA) was performed, followed by batch correction using “Harmony” (“RunHarmony” command, using default settings) (Korsunsky et al., 2019). Nearest Neighbor analysis and clustering were performed using the Harmony embeddings. The clusters were visualized using Uniform Manifold Approximation and Projection (UMAP) (Becht et al., 2018). The DEGs between the clusters were identified using the Wilcoxon Rank Sum test implemented in Seurat’s “FindAllMarkers” command. The cell identities were then assigned according to the markers identified (Supplementary Tables S1, S2). Cell types that were over clustered were merged (e.g., mast cells). For pseudotime trajectory analysis, a new Seurat object containing the osteoblasts clusters (1–4) and the skeletal progenitor cell was generated. The single cell trajectory and branching point were detected using the reversed graph embedding algorithm implemented in Monocle3 R package (Qiu et al., 2017) using default parameters. The root of the trajectory was defined as the list of the skeletal progenitor cells. The full code for the analysis and figures will be found in https://github.com/ishainer/R_shainer_et_al_2022/. The raw scRNA seq data is available through the NCBI GEO repository, accession number GSE192542, upon publication (token will be given upon reviewer’s request).

Mouse serum was analyzed from peripheral blood, obtained by retro-orbital bleeding, 24 h post fracture. Circulating cytokine levels were determined by flow cytometry using a mouse inflammation CBA kit (BD Bioscience, BD 552364) according to the manufacturer’s protocol. Briefly, a mix of six bead populations coated with antibodies specific for IL-6, IL-10, MCP-1, IFN-γ, TNF, and IL-12p70 were incubated for 2 h at RT with serum samples and with PE conjugated detection antibodies to form sandwich complexes. PE fluorescence intensity for each of the antibodies was measured by Flow Cytometer and FCAP Array™ software was used to generate the results.

Femurs were embedded without demineralization in SamplKwick fast cure acrylic compound (Burhler). 7 mm thick cross-sections were cut from the bone’s mid-shaft using a diamond saw. The samples were ground and polished with Micro-Mesh sanding sheets followed by 0.25 µm diamond paste. Indents were preformed using HMV-G21DT micro-hardness testing machine (Shimadzu) with a diamond Vickers microindenter tip using 0.05 N of applied force for 10 s before unloading. 3-4 indents were performed for each sample and the results were averaged.

Differences were examined by two-tailed Student’s t-test for comparing two groups and by either one-way or two-way analysis of variance (ANOVA) test for comparing multiple groups. When significant differences were indicated by ANOVA, group means were compared to establish the source of the differences. p < 0.05 was considered statistically significant.

Bgn is an abundant matrix component expressed in skeletal tissues throughout bone development and maturation. To study the role of Bgn in vivo, a Bgn-deficient mouse line was generated (Xu et al., 1998). Although Bgn KO mice appear normal at birth, they display a phenotype characterized by reduced bone mass and age-dependent osteopenia. Lechner et al. (2006) measured the expression of Bgn mRNA in mouse hindlimbs during development and found that Bgn was present at E14, exhibited a striking 5-fold upregulation at E16, before dropping 5-fold by E18. Therefore, to understand when bone defects begin to arise in Bgn KO mice, we isolated cells from fetal hindlimbs of WT and Bgn KO embryos at 16.5 days of age, where Bgn transcript level is at its highest, and performed single-cell RNA sequencing (scRNA seq) analysis using the 10X Genomic’s Chromium platform. Clustering analysis of the sequenced cells revealed 38 transcriptionally distinct populations. The top 50 differentially expressed genes within each cluster were used to identify the cell-type represented by each cluster (Figure 1A; Supplementary Table S1). 14,781 WT cells and 15,054 Bgn KO cells were analyzed, and both genotypes expressed the same clustering pattern (Figure 1B; Supplementary Figures S1A, S2). Bgn was found to be expressed mainly in skeletal progenitor cells, osteoblasts, chondrocytes, stromal cells, embryonic fibroblasts, tendon cells and neuronal progenitor cell clusters (Figure 1C). Bgn mRNA wasn’t detected in Bgn KO cells as shown in Figure 1C and Supplementary Figure S1B. All the clusters included cells originating from both WT and Bgn KO animals, but the relative distribution of Bgn KO and WT cells differed in many of the clusters (Figure 1D; Supplementary Figure S2). From the total analyzed cells, 8.6% were chondrocytes. However, 5.6% came from WT cells and 3% came from Bgn KO cells. In addition, from the total analyzed cells, 18.6% were osteoblasts. In this cell type, 12.1% came from WT cells while only 6.5% came from Bgn KO cells. More specifically, in osteoblast clusters 1 and 2, the number of WT cells was almost equal to the number of Bgn KO cells but in the progression to osteoblast clusters 3 and 4 the number of Bgn KO cells dramatically decreased. In particular, osteoblast cluster 4 was comprised primarily of WT cells (88%) versus Bgn KO cells, which represented a minor fraction (12%) (Figures 1D–E; Supplementary Figure S2). Since the total number of cells that was analyzed equally represent both genotypes, the differences in specific cluster distribution most likely reflect the effect of Bgn deficiency on cell differentiation. To identify the differences between the 4 osteoblast clusters, enrichment analysis of biological pathways was performed on the differentially expressed genes. While the cells in all osteoblast clusters expressed genes that are involved in the ECM organization and the skeletal system development pathways, the genes specifically associated with the ossification pathway were highly expressed in osteoblast clusters 3 and 4, and genes associated with the osteoblast differentiation pathway were only found in osteoblast cluster 4 (Figure 1F). To further investigate the relationship between the osteoblast clusters, we analyzed their differentiation progress along a pseudotime trajectory using monocle3 (Qiu et al., 2017). The pseudotime trajectory was tested for the four osteoblast clusters, together with the skeletal progenitor cells, that were defined as the trajectory root. The global gene expression changes, enabled the detection of a single branching point leading the cell fate towards osteoblasts 3 and osteoblasts 4 (Supplementary Figures S1C, D), suggesting these clusters differentiate later than osteoblast cluster 2.

Taken together these data suggest that the differentiation and maturation levels of the osteoblasts are more advanced as they progress from osteoblast cluster 1 to osteoblast cluster 4. Cells from embryonic Bgn KO bone have less mature osteoblasts, less chondrocytes and more immature embryonic fibroblast cells (Figure 1D; Supplementary Figure S2), suggesting that Bgn is needed for early stages of bone development starting from the origin of the cells that generate the bone.

To understand how the changes we found in the scRNA seq of embryonic bones could affect skeletal development, µCT was performed on mature 6 week-old mice.

Our data showed that there was a reduction in cortical thickness and trabecular number, and an increase in trabecular spacing and bone mineral content (BMC) in Bgn KO mice compared with WT mice (Figures 2A, B). Since 90% of the organic ECM of bone is collagen type I, second harmonic generation (SHG) microscopy was used to visualize and quantify the structural integrity of the collagen fibers in the Bgn-deficient bones. High collagen signature, as indicated by a consistently high even SHG intensity, with minor spacing between the collagen fibers was observed in the femurs of WT mice, whereas Bgn KO femurs showed significantly less of a collagen signature and more spacing between the fibers (Figures 2C, D). To understand how these molecular structural changes were affecting bone hardness, we assessed the strength of the bone by micro-indentation and found a reduction in the hardness of the Bgn KO bones compared with WT (Figure 2E). These results correlate with previous observations by Wallace et al. (2006) who used a 4-point bending test to demonstrate a reduction in the bone strength in Bgn KO bones compared to WT. Taken together, our findings show that in the absence of Bgn, compromised bone formation occurs that, ultimately, appears to affect collagen integrity, and subsequently biomechanical phenotype.