A systematic dissection of human primary osteoblasts in vivo at single-cell resolution

Osteoblasts identification

We applied an established protocol for in vivo human osteoblast isolation [19] to obtain the alkaline phosphatase (ALPL)^high/CD45/34/31^low cells from the femur head-derived bone tissue of one human subject (31-year-old Chinese male) with osteoarthritis and osteopenia through fluorescence-activated cell sorting (FACS). Several studies [20, 21] have demonstrated that this isolation protocol can successfully recover a high proportion of osteoblasts based on the elevated expression levels of osteoblastic markers via quantitative polymerase chain reaction (qPCR) or bulk RNA sequencing. In total, 9,801 single cells were encapsulated for cDNA synthesis and barcoded using the 10x Genomics Chromium system, followed by library construction and sequencing (Figure 1A and Supplementary Figure 1A [22]). After quality control (Supplementary Figure 1B [22]), we obtained 8,557 cells, with an average of 2,365 and median of 2,260 genes detected per cell (Supplementary Table 1 [22]). A recently developed dimension reduction technique for scRNA-seq analysis, uniform manifold approximation and projection (UMAP) [23], was applied to project the gene expression profiles on a two-dimensional panel for visualization of cellular heterogeneity (Supplementary Figure 1C [22]). Compared with t-SNE [24], another commonly used dimension reduction method in single cell analysis, UMAP provides the useful and intuitively pleasing feature that it preserves more of the global structure [23], which is important to identify the systematic differences of gene expression patterns between different clusters through the inter-cluster distances. UMAP has already been successfully and now commonly used for scRNA-seq data dimension reduction in recent scRNA-seq studies [25–27]. After clustering the cells into six distinct subsets (C1-C6) by the k-nearest neighbor algorithm [28], we used pairwise differential expression analysis for comparing each individual cluster against all the others to identify the DEGs of each subtype (Supplementary Figure 1D [22]).

Figure 1. scRNA-seq analysis of human osteoblasts. (A) Study overview. (B) Three osteoblast clusters. UMAP visualization of 5,329 osteoblasts, colored by clustering. (C) Proportion of three osteoblast clusters. Colored by clustering. (D) Cluster signature genes. Violin plots showing the log-transformed normalized expression levels of the two most significant marker genes in clusters O1, O2, and O3, respectively. (E) Log-normalized expression of adipocyte and chondrocyte biomarkers in osteoblast clusters. (F) Log-normalized expression of LEPR⁺ mesenchymal cell related markers in osteoblast clusters. (G) Log-normalized expression of OPG, RANKL, IL6, and IL7 in three clusters.

While ALPL was enriched in clusters C1, C2 and C5, the osteoblast specific marker, RUNX2 (encoded by RUNX2) [29], was only enriched in clusters C1 and C2, suggesting the presence of contamination by other cell types during the cell isolation process (Supplementary Figure 1D [22]). Based on the identified marker genes and gene ontology (GO) enrichment analysis of the differentially expressed genes in each cluster (Supplementary Figure 1D–1F [22]), the contaminant cells were classified as; 1) two nucleated erythrocyte groups (C3 and C4, expressing hemoglobin coding genes HBB and HBA1) [30]; 2) one smooth muscle cell group (C5, expressing transgelin coding gene, TAGLN [31], and CNN1, which is specific to differentiated mature smooth muscle cells) [32]; and 3) one neutrophil group (C6, expressing neutrophil related genes S100A8 and ELANE) [33, 34]. These findings suggest that the protocol for human osteoblast isolation based on ALPL^high/CD45/34/31^low may not be sufficient for specifically isolating osteoblasts alone, since ALPL is also highly enriched in other mesenchymal cell-derived cells such as vascular smooth muscle cells [35]. Notably, by comparing the expression profiles between osteoblasts and contaminant cells, we found that the average fold change (avg_FC) of alpha-1 type I collagen (COL1A1, encoded by COL1A1) was large (avg_FC =4.497, Supplementary Table 2 [22]), suggesting that positive selection based on the combination of ALPL and COL1A1 would be a more appropriate choice for human primary osteoblast isolation. While several studies [17, 36] in mice have used COL1A1 for osteoblast sorting, few (if any) studies have adopted this marker for human primary osteoblast isolation.

scRNA-seq identifies multiple cell subtypes in human osteoblasts

To investigate the cellular diversity within osteoblasts, we extracted the clusters C1 and C2 with high expression levels of osteoblast related markers, i.e., RUNX2 and COL1A1 (Supplementary Figure 1D [22]).

To ensure the high quality of the osteoblasts data for downstream analysis, we further removed osteoblasts with more than 5% of the transcripts attributed to mitochondrial genes, following the quality control criterion applied in several scRNA-seq studies in bone field [37–39]. 5,329 high quality sequenced osteoblasts were retained for downstream analysis. We noticed that the biomarkers of other mesenchymal-derived cells such as adipocytes and chondrocytes were not enriched in the remaining cells [38] (Figure 1E). Further, although the remaining cells highly expressed the leptin receptor positive (LEPR⁺) mesenchymal cell related markers EBF transcription factor 3 (encoded by EBF3), C-X-C motif chemokine ligand 12 (encoded by CXCL12), and platelet derived growth factor receptor alpha (encoded by PDGFRA) [40] (Figure 1F), the high enrichment of osteoblast related markers (e.g., COL1A1, ALPL, RUNX2, etc.) suggests that these cells are actually osteoblasts rather than LEPR⁺ mesenchymal cells or CXCL12-abundant reticular (CAR) cells [41]. We identified three cell subtypes of osteoblasts (Figure 1B–1D) based on their transcriptional profiles and by comparing classical osteoblastic markers [42] with three published scRNA-seq datasets of mouse osteoblasts [16, 17, 41]. The cell subtypes were annotated as; 1) O1 (65.77%), pre-osteoblasts, with relatively high expression levels of LEPR⁺ mesenchymal cell markers LEPR (encoded by LEPR) [43] and vascular cell adhesion molecule 1 (encoded by VCAM1) [44]; 2) O2 (25.39%), mature osteoblasts, which showed highest expression levels of COL1A1 (COL1A1) and osteogenesis-associated genes including BGLAP (also known as osteocalcin, encoded by BGLAP), SPP1 (also known as osteopontin, encoded by SPP1), IBSP (encoded by IBSP) and interferon induced transmembrane protein 5 (encoded by IFITM5) [45, 46]; 3) O3 (8.84%), undetermined osteoblasts, which not only expressed several LEPR+ mesenchymal cell-associated genes (e.g., LEPR and VCAM1) but are also distinguished from the other subtypes by distinctively expressing high levels of nuclear receptor subfamily 4 group A member 1 and 2 (encoded by NR4A1 and NR4A2). Although pre- and mature osteoblasts are vague concepts largely suggested by in vitro studies rather than discrete in vivo cell types, we use these terms to better describe the specific characteristics of clusters O1 and O2. Notably, few osteoblasts in this cell population expressed OPG and RANKL (Figure 1G). This is consistent with the result proposed by Tat et al. [47] that the expression levels of OPG and RANKL are significantly reduced in the osteoblasts from human osteoarthritic bone. In addition, we also noticed the high expression of the inflammation markers IL6 and IL7 in cluster O3 and O1, respectively (Figure 1G). Previous studies have proposed that IL6 is one of several pro-inflammatory cytokines present in individuals with confirmed clinical diagnosis of osteoarthritis [48, 49].

We further examined the transcriptional profiles of the three identified osteoblast subtypes and found that in addition to the cell markers, some bone development regulators were highly enriched in the pre-osteoblast and mature osteoblast clusters. For instance, pre-osteoblasts expressed insulin-like growth factor-binding protein 2 and 4 (IGFBP2 and IGFBP4, encoded by IGFBP2 and IGFBP4) (Figure 2A). Previous studies have considered IGFBP2 as a stimulator for osteoblast differentiation through positive regulation of the AMP-activated protein kinase (AMPK) [50], and demonstrated that IGFBP4 could stimulate adult skeletal growth in males [51]. Meanwhile, tenascin, a glycoprotein encoded by TNC that modulates osteoblast mineralization via matrix vesicles [52], was highly enriched in mature osteoblasts. On the other hand, the cluster O3 showed high expression level of osteomodulin (OMD) (Figure 2A). It has been proposed that OMD induces endochondral ossification through PI3K signaling, which is an essential process for long bone formation [53].

Figure 2. Osteoblast subtypes and cellular functions in bone formation. (A) Osteoblasts related genes expressed in clusters O1, O2, and O3, respectively. (B–D) GO enrichment for the three osteoblast subpopulations, O1, O2, and O3, respectively. The length of the bar indicates the gene ratio (number of DEGs enriched in the GO term / total number of DEGs). The color indicates the adjusted p values for enrichment analysis. (E) Bone formation related GO terms enriched in clusters O1, O2, and O3. The size of dot indicates the gene ratio, which is the ratio of functional related genes and the total number of the differential expression genes compared with other clusters. The color indicates the adjusted p-value for enrichment analysis. (F) Immunofluorescence of mouse femur. The osteoblast marker ALPL was stained by green, while the cluster O3 marker NR4A1 was stained by red. The undetermined osteoblasts were located on the bone surface, co-stained by green and red (yellow).

In addition to the known osteoblast subtypes (pre-osteoblast and mature osteoblast), we also identified one potential novel osteoblast subtype (cluster O3), which is a major source for orphan nuclear receptors coding genes NR4A1 and NR4A2. Although previous studies have reported the presence of NR4A1 in mice osteoblasts induced by parathyroid hormone (PTH) in vitro [54], the enrichment of this gene in osteoblasts in vivo has not yet been proposed or identified earlier in mice or in humans. Several studies have reported that NR4A receptors have an essential impact on bone metabolism by regulating the expression of several osteoblastic marker genes such as SPP1, BGLAP, COL1A1, and ALPL [55–57]. Additionally, the overexpression of NR4A2 led to increased expression of ALPL, COL1A1, and BGLAP in osteoblasts [55–57]. The immunofluorescence on mouse femur illustrated the co-staining of the osteoblast marker ALPL with the biomarker NR4A1 on the bone surface, thereby supporting the expression of NR4A1 in osteoblasts in vivo (Figure 2F).

While all the osteoblast clusters highly expressed the osteoblastogenic gene RUNX2, they may affect bone development from different aspects. To further investigate the distinct biological functions related to osteogenesis among the three clusters of osteoblasts, we performed GO enrichment analysis based on the DEGs (avg_FC > 1.200 and p_val_adj < 0.05, Supplementary Table 3 [22]) in each cluster. All the clusters were enriched in GO terms related to bone development including “extracellular matrix organization” and “extracellular structure organization” with a relatively high gene ratio (Figures 2B–2E). The pre-osteoblasts and mature osteoblasts are involved in the extracellular matrix (ECM) formation through production of different types of collagen, such as types 3, 14, and 18 in pre-osteoblasts and types 1, 12, and 13 in mature osteoblasts (Supplementary Tables 4, 5 [22]). The mature osteoblasts were uniquely enriched for “bone mineralization”, while mature and undetermined osteoblasts were both enriched for “osteoblast differentiation” (Figure 2C–2E and Supplementary Tables 5, 6 [22]). Additionally, the undetermined osteoblasts were uniquely enriched for “regulation of osteoblast differentiation” and “positive regulation of osteoblast differentiation” (Figure 2D). Surprisingly, few GO terms related to osteoblast differentiation are enriched in the pre-osteoblasts, suggesting that the undetermined and mature osteoblasts may modulate the osteoblast differentiation processes to a larger extent than pre-osteoblasts.

Dynamic gene expression patterns reveal the differentiation relationship between different osteoblast subtypes

In order to reveal the differentiation dynamics of the osteoblast cell population, we reconstructed the developmental trajectory of the three identified clusters of osteoblasts. Since all the three clusters expressed the LEPR⁺ mesenchymal cell related markers EBF3, CXCL12, and PDGFRA (Figure 1F), we suggested that all osteoblasts are in the LEPR cells derived osteoblastic lineage. All cells were contained within one cellular lineage without any bifurcations (Figure 3A), suggesting that only one terminal subtype exists in the osteoblastic population. While pre-osteoblasts and undetermined osteoblasts were highly enriched in the early stages of pseudotime, mature osteoblasts were distributed in the terminal stages of the osteoblastic lineage (Figure 3A, 3C). To strengthen the trajectory inference, we also assessed the transcriptional continuum of the cell lineage. Since the scRNA-seq data is difficult to model with a parametric curve (e.g., least squares regression, etc.) due to the high sparsity and variability, we used the loess curves to reveal the trends of gene expression along with the pseudotime. The results showed that the expression of the LEPR⁺ mesenchymal cell associated genes LEPR and VCAM1 [58] as well as cluster O3 specific markers, NR4A1 and NR4A2, decreased over pseudotime, while the expression levels of osteogenic markers such as RUNX2, BGLAP, SPARC [59], and COL1A1 were the highest at the end stages of pseudotime (Figure 3B). This result is consistent with the findings from other studies that have examined the transcriptional continuum in osteoblastic differentiation [42]. The analysis of the gene expression profiles of mouse osteoblasts cultured in vitro showed that the expression levels of NR4A1 rapidly increased from day 5 to day 7 whereas the expression of ALPL remained relatively stable throughout the osteoblast development (Figure 3D), further supporting the conclusion that the undetermined osteoblasts are involved in the early stages of the osteoblastic lineage.

Figure 3. Trajectory inference of human osteoblasts. (A) Reconstructed cell differentiation trajectory of human osteoblasts. The upper-right trajectory plot in the square indicates the direction of pseudotime. (B) Expression levels (log-normalized) of indicated genes in the three osteoblast subtypes with respect to their pseudotime coordinates. The x-axis indicates the pseudotime, while the y-axis represents the log-normalized gene expression levels. The color corresponding to the three different osteoblast subsets. The loess regression was applied to fit the relationship between pseudotime and expression level. (C) Cell distribution based on the pseudotime coordinates. The x-axis is the pseudotime and the y-axis represents the osteoblast subtypes. (D) Expression levels of ALPL and NR4A1 in mouse osteoblasts in vitro at day 0, 5 and 7, respectively. N.S., not significant, *p-adjusted ≤ 0.05, **p-adjusted≤ 0.01, *** p-adjusted ≤ 0.005.

Pre-osteoblasts regulate angiogenesis through multiple signaling pathways

Although evidence has shown that osteoblasts can regulate angiogenesis [11], few studies have explored this relationship at the single-cell level. Based on functional enrichment of the highly expressed genes, we found that the GO term “regulation of angiogenesis” was enriched in pre-osteoblasts with a gene ratio of 0.10 but not in other clusters (Figure 4A). We further investigated the enriched genes in pre-osteoblasts and found that 22 genes are involved in angiogenesis (Supplementary Table 4 [22]). In particular, the pre-osteoblasts showed significantly higher expression levels of SFRP1, MDK, and THBS1 compared with the other clusters (Figure 4B). It has been demonstrated that Secreted Frizzled-related protein-1 (sFRP1, encoded by SFRP1), a modulator of the Wnt/Fz pathway, can modify mesenchymal cell capacities enabling mesenchymal cells to increase vessel maturation [60]. Midkine (MDK, encoded by MDK) is an enhancer of angiogenesis, and MDK expression has been shown to be positively correlated with vascular density in bladder tumors [61]. Thrombospondin 1 (THBS1, encoded by THBS1) is considered to be a potent endogenous inhibitor of angiogenesis through antagonization of vascular endothelial growth factor (VEGF) activity [62].

Figure 4. Potential functions of pre- and undetermined osteoblasts in angiogenesis and hemopoiesis. (A) Angiogenesis and hematopoiesis modulation related GO terms enriched in the three clusters. The x-axis represents the clusters and the y-axis is the GO terms related to the angiogenesis and hematopoiesis regulation. The size of the dot indicates the gene ratio and the color indicates the adjusted p-values. (B) Angiogenesis associated genes enriched in cluster O1. X-axis represents the three clusters and y-axis reflects log-normalized gene expression levels. Stars indicate the significance levels of the gene expression difference between two clusters by Wilcoxon signed-rank test. N.S., not significant, *p-adjusted ≤ 0.05, **p-adjusted ≤ 0.01, *** p-adjusted ≤ 0.005. (C) Gene expression pattern in enriched pathways. Squares show enriched DEGs in the corresponding terms (rows). Color indicates the expression value of the DEGs (average logFC). (D) Hemopoiesis factors enriched in clusters O1 and O3.

Next, we investigated the signaling pathways enriched in pre-osteoblasts that are related to angiogenesis. Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis revealed three highly enriched pathways including PI3K-Akt, phospholipase D, and relaxin signaling pathways, which are known to be implicated in angiogenesis modulation [63–65] (Figure 4C). Several genes related to the “regulation of angiogenesis” GO term were enriched in these pathways (Supplementary Table 4 [22], Figure 4C) including COL4A2, TGFBR2, and IL7. These findings suggest that pre-osteoblasts may potentially regulate angiogenesis to a larger extent than other osteoblast subtypes.

Osteoblast populations modulate the development of HSC niche

Osteoblasts are part of the stromal cell support system that provides critical regulators of hematopoiesis [10]. It has been reported that osteoblasts lining the bone surface can generate an extracellular environment supporting non-skeletal hematopoietic cells [66]. To identify the specific roles of each osteoblast subpopulation in hematopoiesis regulation, we first examined the expression patterns of hematopoiesis related factors in all three osteoblast subtypes. We found that while both pre-osteoblasts and undetermined osteoblasts highly expressed insulin-like growth factor 1 (IGF-1, encoded by IGF1), which is a growth factor that has been shown to be responsible for burst-like growth of early erythroid progenitor cells in vitro [67], the pre-osteoblasts may also play an important role in modulating macrophage development due to the high expression of M-CSF (encoded by CSF1, Figure 4D) [68]. Further, the high enrichment of the stem cell factor (SCF, encoded by KITLG, Figure 4D) in pre-osteoblasts suggests that these cells may be crucial for preventing HSC apoptosis [69]. Since Aguila et al. [70] proposed that the overexpression of human interleukin 7 (IL-7, encoded by IL7) increases the development of B cells in female mice, it is plausible that the pre-osteoblasts, with high expression of IL7 (Figure 4D), may affect the development of B cells in humans.

We also found that the undetermined osteoblasts represent a major source of C-C motif chemokine ligand 2 (CCL-2, encoded by CCL2) (Figure 4D), which is a critical factor for monocyte recruitment in acute inflammatory response [71]. Additionally, the undetermined osteoblasts expressed interleukin 6 (IL-6, encoded by IL6), which displays a strong synergy with IL-7 to stimulate naïve CD8⁺ T cells [72]. Therefore, the undetermined osteoblasts may coordinate with pre-osteoblasts to support the hemopoietic niche. However, few hematopoietic factors were enriched in mature osteoblasts, suggesting that they may have limited impact on the hematopoietic system.

We investigated the predicted functions associated with hematopoiesis regulation enriched in the pre-osteoblasts and undetermined osteoblasts. Notably, while the GO term “positive regulation of hemopoiesis” was enriched in pre-osteoblasts, the highly expressed genes of the undetermined cluster were related to the negative regulation of hemopoiesis (Figure 4A). Apart from the growth factors mentioned above, pre-osteoblasts also had a significantly higher expression level of growth arrest-specific 6 (GAS6, encoded by GAS6) (Figure 4D), which can induce natural killer cell development via its positive regulatory effect on fibromyalgia syndrome-like tyrosine kinase 3 (FLT3) signaling in CD34+ hematopoietic progenitor cells (HPCs) [73]. Meanwhile, the undetermined subgroup had significantly higher expression levels of zinc-finger protein 36 (ZFP36, encoded by ZFP36) (Figure 4D), which can inhibit the erythroid differentiation [74]. These findings suggest that these two osteoblast subgroups may play different roles in hemopoiesis regulation, and that the relative proportions and functional levels of these subtypes within the osteoblast cell population may be crucial for maintaining the homeostasis of hematopoiesis.

Transcriptional divergence of human and mouse osteoblasts

Due to the limited in vivo transcriptomics studies in mouse osteoblasts at different developmental stages, we first compared the expression profiles of osteoblasts acquired from humans in vivo and the osteoblasts from mice cultured in vitro. The in vitro mouse osteoblast data demonstrate the gene expression levels at three different time points (Day 5, 14, and 21), which is considered as the reference for us to identify the divergence of transcriptional profiles between human and mouse osteoblasts in multiple development stages. Several fundamental differences were observed between humans and mice. For instance, all three human osteoblast clusters in vivo expressed secretory leukocyte protease inhibitor (encoded by SLPI) (Figure 2A), a serine protease inhibitor that promotes cell migration and proliferation while also suppressing the inflammatory response [75]. This conflicts with the previous in vitro findings in mice [76], which demonstrated the enrichment of SLPI in mature but not in pre-osteoblasts. Additionally, we found that the human pre-osteoblasts isolated in vivo showed a significantly higher expression level of SLPI compared with mature osteoblasts (Figure 2A). In further contrast, the mature mouse osteoblasts cultured in vitro were highly enriched for myristoylated alanine rich protein kinase C (PKC) substrate (MARCKS, encoded by MARCKS), while this gene was highly expressed by the pre-osteoblasts but not the mature osteoblasts in humans in vivo (Figure 2A). It has been shown that MARCKS is the PKC-δ effector which modulates cathepsin K secretion and bone resorption in osteoclasts [77]. Remarkably, we also found that CPE, encoding carboxypeptidase E, was only enriched in mature osteoblasts in vivo in humans (Figure 2A), which is inconsistent with the previous mouse study in vitro finding that the expression level of CPE decreases throughout osteoblast development [76]. While CPE plays an important role in RANKL-induced osteoclast differentiation [78], it is plausible that a reduction in the proportion of mature osteoblasts could attenuate osteoclast proliferation.

To better understand the systematic differences in osteoblastic transcriptional profiles between human and mouse, we then integrated our sequencing data with the publicly available data from a previous scRNA-seq study [17] on mouse osteoblasts in vivo (Figure 5A). After correcting for the potential batch effects between independent experiments using canonical correlation analysis (CCA) [28], a moderate positive correlation (R = 0.53, p-value < 2.23e-16) was observed between the human and mouse osteoblastic transcriptomic profiles (Figure 5B). This suggests that while some common characteristics exist between human and mouse osteoblasts, there are considerable differences in the gene expression profiles at the single-cell level.

Figure 5. Integrated cross-species analysis between human and mouse osteoblasts in vivo. (A) UMAP visualization of human and mouse osteoblast integration. The upper-left plot represents the three osteoblast clusters in mice, the upper-right plot indicates the three osteoblast subtypes in humans. The bottom plot represents the five clusters after the integration, colored by clusters. The middle-right is the integrated data colored by different species. (B) Correlation of gene expression among different osteoblast datasets after CCA integration. Each dot represents an individual gene. Axes measure the average gene expression levels (logFC) in the indicated subject. Correlations were measured by Spearman correlation coefficients (R) (p < 0.01). (C) Log-normalized expression of NR4A1 in human and mouse integrated data. (D) Proportion of the human and mouse osteoblasts in each cluster after CCA integration, colored by different species. The x-axis represents five different clusters and y-axis indicates the proportion.

To further investigate the shared and distinct features among human and mouse osteoblasts, we applied unbiased clustering analysis after dimension reduction by UMAP [23]. The result showed that human osteoblast cluster O2 overlapped with mouse osteoblast subtype M2, which represents the mature osteoblasts in mice [17] (Figure 5A, 5C), indicating that the overall gene expression patterns of human and mouse mature osteoblasts are highly similar. Furthermore, we also found strong correlations in the expression patterns between small subgroups in human pre-osteoblasts (cluster O1) and the mouse osteoblast subtype M1 [17] (Figure 5A, 5C). In contrast, the human osteoblast cluster O3 did not overlap with any osteoblast subtypes in mice (Figure 5A). However, we did observe that NR4A1 was expressed in cluster O1 and M1 and M3 (Figure 5B). Since few human osteoblasts in cluster O1 expressed NR4A1 (Figure 1D), mouse osteoblast cluster M1 is the major source of the NR4A1 in the cluster O1 and M1. Therefore, although NR4A1 is expressed by both human and mouse osteoblasts, heterogeneity still exists between human and mouse NR4A1⁺ osteoblasts. These findings suggest that the transcriptional profiles of human and mouse osteoblasts may demonstrate systematic differences in early cellular developmental stages but share similar features at the terminal stage of osteoblastic development lineage. In addition, the comparison between our data and osteoblasts from mice cultured in vitro [76] revealed that several genes (e.g., SLPI, MARCKS, CPE etc.), which are not correlated with osteoarthritis pathogenesis, show fundamental differences of expression levels between human and mouse osteoblasts. This suggests that the systematic variations may still exist between healthy human and mouse osteoblasts. Therefore, studies based on osteoblasts acquired from mouse models may introduce notable bias when inferring the biological characteristics of human osteoblasts, especially for those that have not yet reached a mature state.

Study population

The clinical study was approved by the Medical Ethics Committee of Central South University and written informed consent was obtained from the study participant. The study population included one 31-year-old male with osteoarthritis and osteopenia (BMD T-score: 0.6 at lumbar vertebrae, -1.1 at total hip), who underwent hip replacement at the Xiangya Hospital of Central South University. The subject was screened with a detailed questionnaire, medical history, physical examination, and measured for bone mineral density (BMD) before surgery. Subject was excluded from the study if he had preexisting chronic conditions which may influence bone metabolism including diabetes mellitus, renal failure, liver failure, hematologic diseases, disorders of the thyroid/parathyroid, malabsorption syndrome, malignant tumors, and previous pathological fractures [82]. The femur head was collected from the patient during hip replacement surgery. The specimen was shortly stored in 4° C and transferred to the laboratory within 30 mins, where it was then processed immediately after delivery.

Mice

2 months old female C57BL/6J wild-type mice were purchased from Jackson Laboratory (Bar Harbor, ME, USA). All mice were housed in pathogen-free conditions and fed with autoclaved food. All experimental procedures were approved by the Ethics Committee of Xiangya Hospital of Central South University.

BMD measurement

BMD (g/cm²) was measured using DXA fan-beam bone densitometer (Hologic QDR 4500A, Hologic, Inc., Bedford, MA, USA) at the lumbar spine (L1 –L4) and the total hip (femoral neck and trochanter) as described previously by our group [83, 84]. According to the World Health Organization definition [85] and the BMD reference established for Chinese population [86], subject with a BMD of 2.5 SDs lower than the peak mean of the same gender (T-score ≤ −2.5) were determined to be osteoporotic, while subject with -2.5 < T-score < -1 are classified as having osteopenia and subject with T-score > -1.0 are considered healthy.

Isolation of osteoblasts

Osteoblasts were extracted from the human femur head based on the widely used dissociation protocols [19] with a few minor adjustments. Briefly, bone tissue samples were chopped into small fragments and washed twice with phosphate-buffered saline (PBS). These fragments, containing a mixture of cortical and trabecular bone, were then incubated with a highly purified, endotoxin-free type II collagenase (1 mg/ml, Cat: A004174-0001, Sangon Biotech, Shanghai, China) at 37.0° C for 30 mins and 60 mins in the first and second round of digestion respectively. After the incubation, the solution was filtered through a 70 μm filter and incubated with red blood cell lysis buffer (Cat: R1010, Solarbio Science and Technology CO., Beijing, China) for 5 mins. The collected cells were washed twice with PBS.

Fluorescence-activated cell sorting (FACS) enrichment for osteoblasts

Before FACS, collected cells were incubated with human CD31/34/45-PE (Cat:303106, Cat:343606, Cat:304008, BioLegend, San Diego, CA, USA), human ALPL-APC (Cat: FAB1448A, R&D Systems, Minneapolis, MN, USA), and 7-AAD (488 nm) antibody (R&D Systems, Minneapolis, MN, USA). FACS was performed on a BD FACS Aria II sorter (Becton-Dickinson, Franklin Lakes, NJ, USA), where dead cells and debris were excluded by forward scatter (FSC) versus side scatter (SSC), and live cells were further enriched by negative selection of 7-AAD. ALP+/CD45/34/31- cells were collected as osteoblasts (Supplementary Figure 1) [19]. The percentage of live cells was 80% after filtering the dead cells. Gating schemes were established with fluorescence-minus-one (FMO) controls.

Single-cell RNA-seq (scRNA-seq) library preparation and sequencing

scRNA-seq libraries were prepared using Single Cell 3’ Library Gel Bead Kit V3 following the manufacturer’s guidelines (https://support.10xgenomics.com/single-cell-gene-expression/library-prep/doc/user-guide-chromium-single-cell-3-reagent-kits-user-guide-v3-chemistry). Single cell 3’ Libraries contain the P5 and P7 primers used in Illumina bridge amplification PCR. The 10x Barcode and Read 1 (primer site for sequencing read 1) were added to the molecules during the gel bead-in emulsions reverse transcription (GEM-RT) incubation. The P5 primer, Read 2 (primer site for sequencing read 2), Sample Index and P7 primer were added during library construction. The protocol is designed to support library construction from a wide range of cDNA amplification yields spanning from 2 ng to 2 μg without modification. All constructed single-cell RNA-seq libraries were sequenced on the Illumina Novaseq6000 platform with a sequencing depth of at least 100,000 reads per cell for a 150bp paired end run.

Pre-processing of single-cell RNA-seq data

We demultiplexed the cellular barcodes and aligned reads to the human transcriptome (GRCh38/hg38) using Cell Ranger 3.0 (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger). To create Cell Ranger-compatible reference genomes, the references were rebuilt according to instructions from 10x Genomics (https://www.10xgenomics.com), which performs alignment, filtering, barcode counting and unique molecular identifier (UMI) counting. Next, a digital gene expression matrix (gene counts versus cells) was generated and converted to a Seurat object by the Seurat R package [28]. For quality control, we removed the cells with <200 genes or >5,000 genes detected, or cells where >15% of the transcripts were attributed to mitochondrial genes based on the quality control metrics (Supplementary Figure 1B) and criteria used by previous scRNA-seq studies [87–89]. We normalized the filtered gene expression matrix by the NormalizeData function in Seurat R package, in which the number of UMIs of each gene was divided by the sum of the total UMIs per cell, multiplied by 10,000, and then transformed to log-scale.

Dimension reduction and osteoblastic clusters identification

For data visualization and classification, we projected the normalized gene expression matrix on a two-dimensional panel. The 2,000 genes with the highest dispersion (variance/mean) were selected for principal component analysis (PCA). The first 18 principal components (number of PCs was chosen based on standard deviations of the principal components, corresponding to the plateau region of ''elbow plot'') were used for uniform manifold approximation and projection (UMAP) [23] dimension reduction. After data visualization, we applied an unbiased graph-based clustering method [90] with the resolution of 0.1 for clustering analysis. To identify the differentially expressed genes (DEGs) in each cluster, we used the Wilcoxon Rank-Sum test to find the genes showing significantly higher levels of expression (false discovery rate (FDR) < 0.05) in a specific cluster compared to the other clusters. After osteoblast identification, we further removed osteoblasts where > 5% of the transcripts were attributed to mitochondrial genes. 2,000 genes with the highest dispersion and the first 18 principal components of osteoblast data were chosen for the UMAP dimension reduction. We then applied the resolution of 0.15 for clustering analysis to reveal three clusters of osteoblasts.

Pathway enrichment and trajectory inference analysis

To investigate the biological processes and signaling pathways associated with each subtype, we performed gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis for the genes identified as important cluster DEGs (avg_FC >1.2000 and p_val_adj < 0.05) by using the clusterProfiler R package [91]. We then used the Monocle 2 v2.8.0 [92] R package to reconstruct the single-cell developmental trajectories in pseudo-time order. The principle of this analysis is to determine the pattern of the dynamic process experienced by the cell population and to order the cells along their developmental trajectory based on differences in the expression profiles of highly variable genes. We followed the official tutorial of Monocle 2 v2.8.0 to perform the pseudotime analysis. Briefly, 1,000 with the highest dispersion were selected as the highly variable genes. Based on the highly variable genes, we then used “reduceDimension” function for the dimension reduction and “orderCells” function with default parameters for the cell ordering.

Cross-species scRNA-seq data integration

One previous scRNA-seq dataset of mouse osteoblasts was acquired from GEO database under the accession numbers of GSE108891 [17]. After acquiring the expression matrix of the osteoblasts in mice, we clustered them into three subtypes (M1, M2, and M3), following the analysis pipeline proposed by Tikhonova et al. [17]. Next, we integrated the scRNA-seq data of mouse and human osteoblasts by canonical correlation analysis (CCA) using Seurat R package [93]. The biological variance of transcriptional profiles across humans and mice was then evaluated based on the Spearman correlation of average gene expression between each of the datasets.

Bone sectioning, immunostaining, and confocal imaging

Freshly dissected femur from female C57BL/6 wild-type mouse was fixed in 4% paraformaldehyde overnight followed by decalcification in 10% EDTA for 1 week. Samples were cut in 5-μm-thick longitudinally oriented sections. After deparaffinize and antigen retrieval, sections were blocked in PBS with 5% bovine serum albumin (BSA) for 1 hour and then stained overnight with goat-anti-Alpl (R&D: AF2910-SP, 10 μg/mL). Rabbit-anti-Nur77 (Proteintech:12235-1-AP, 1:2000) was used as secondary antibodies (from Invitrogen, 1:400). Slides were mounted with anti-fade prolong gold (Invitrogen) and images were acquired with a Zeiss LSM780 confocal microscope.

Osteogenesis induction

Murine mesenchymal stem cells (MSCs) (1.0 × 10⁴ per well) from STEMCELL (Seattle, WA, United States) were plated in 48-well plates and cultured in MesenCult^TM basal expansion medium with 10% 10x Supplement (Stemcell) for 72 h. Next, the cells were rinsed with PBS and the medium was replaced with osteogenic differentiation medium (Stemcell). MSCs cultured in expansion medium were served as the negative control. Half of the medium was changed every 3 days, and cells were harvested at 0, 5, and 7 days after induction.

qRT-PCR analysis

Total RNA was extracted using an RNA extraction kit (Qiagen, Hilden, Germany). and cDNA was synthesized from 1 μg of total RNA by using the Revert Aid First Strand cDNA synthesis kit (Thermo). Then, the cDNA was amplified with iTaq^TM Universal SYBR^® Green Supermix (BioRad, Hercules, CA, USA). in an ABI PRISM® 7900HT System (Applied Biosystems, Foster City, CA, USA). The relative standard curve method (2–∆∆CT) was used to determine the relative gene expression and GAPDH was used as a housekeeping gene for internal normalization. The PCR primers used in this study were as follows:

GAPDH: forward, 5′-CACCATGGAGAAGGCCGGGG-3′,

reverse, 5′-GACG- GACACATTGGGGGTAG-3′;

ALPL: forward, 5′-CCAACTCTTTTGTGCCAGAGA-3′,

reverse, 5′-GGCTACATTGGTGTTGAGCTTTT-3′;

NR4A1:forward, 5′-AGGGTGTGTGTGCATATGGA-3′,

reverse, 5′-CCGCCATCTTTTCCTGTACG-3′;

Data availability

The scRNA-seq data for primary osteoblasts from one human sample can be accessed with accession number under GSE147390 in GEO database. One previous scRNA-seq data of mice osteoblasts used in this study can be accessed with accession number of GSE108891.