Research Paper Volume 13, Issue 16 pp 20629—20650
A systematic dissection of human primary osteoblasts in vivo at single-cell resolution
- 1 Tulane Center for Biomedical Informatics and Genomics, Deming Department of Medicine, School of Medicine, Tulane University, New Orleans, LA 70112, USA
- 2 Department of Orthopedics, Xiangya Hospital, Central South University, Changsha 410008, China
- 3 Laboratory of Molecular and Statistical Genetics, College of Life Sciences, Hunan Normal University, Changsha 410081, China
- 4 Center of Reproductive Health, System Biology and Data Information, Institute of Reproductive and Stem Cell Engineering, School of Basic Medical Science, Central South University, Changsha 410081, China
- 5 School of Basic Medical Science, Central South University, Changsha 410008, China
- 6 Xiangya Nursing School, Central South University, Changsha 410013, China
- 7 Department of Orthopedics and National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, Changsha 410008, China
- 8 Tulane Center for Bioinformatics and Genomics, Department of Biomedical Engineering, Tulane University, New Orleans, LA 70112, USA
- 9 Department of Cell and Molecular Biology, School of Science and Engineering, Tulane University, New Orleans, LA 70112, USA
Received: March 2, 2021 Accepted: June 19, 2021 Published: August 24, 2021https://doi.org/10.18632/aging.203452
How to Cite
Copyright: © 2021 Gong et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Human osteoblasts are multifunctional bone cells, which play essential roles in bone formation, angiogenesis regulation, as well as maintenance of hematopoiesis. However, the categorization of primary osteoblast subtypes in vivo in humans has not yet been achieved. Here, we used single-cell RNA sequencing (scRNA-seq) to perform a systematic cellular taxonomy dissection of freshly isolated human osteoblasts from one 31-year-old male with osteoarthritis and osteopenia after hip replacement. Based on the gene expression patterns and cell lineage reconstruction, we identified three distinct cell clusters including preosteoblasts, mature osteoblasts, and an undetermined rare osteoblast subpopulation. This novel subtype was found to be the major source of the nuclear receptor subfamily 4 group A member 1 and 2 (NR4A1 and NR4A2) in primary osteoblasts, and the expression of NR4A1 was confirmed by immunofluorescence staining on mouse osteoblasts in vivo. Trajectory inference analysis suggested that the undetermined cluster, together with the preosteoblasts, are involved in the regulation of osteoblastogenesis and also give rise to mature osteoblasts. Investigation of the biological processes and signaling pathways enriched in each subpopulation revealed that in addition to bone formation, preosteoblasts and undetermined osteoblasts may also regulate both angiogenesis and hemopoiesis. Finally, we demonstrated that there are systematic differences between the transcriptional profiles of human and mouse osteoblasts, highlighting the necessity for studying bone physiological processes in humans rather than solely relying on mouse models. Our findings provide novel insights into the cellular heterogeneity and potential biological functions of human primary osteoblasts at the single-cell level.
Osteoblasts, which account for 4-6% of resident cells in bone, are derived from multipotent skeletal stem cells (SSCs) through the activation of signaling pathways regulated by osterix (OSX) and Runt-related transcription factor 2 (Runx-2) . In addition, recent studies have demonstrated that periosteal cells  and growth plate chondrocytes  are also significant sources of osteoblasts. Osteoblasts are primarily known for their bone-building functions, including deposition of calcium phosphate crystals (e.g. hydroxyapatite), production of bone matrix constituents (e.g. type I collagen), as well as their ability to secrete a number of important proteins for bone metabolic processes such as integrin-binding sialoprotein (IBSP), secreted phosphoprotein 1 (SPP1), and bone gamma-carboxyglutamic acid-containing protein (BGLAP) . In general, osteoblasts play a crucial role in the mineralization of the bone matrix .
Many osteoblasts ultimately differentiate into osteocytes, which become embedded in the bone matrix to form a communication network for the regulation of bone formation and resorption. The osteoblasts and osteocytes regulate the differentiation of osteoclasts, which are primarily involved in bone resorption activities . For instance, osteoblasts can promote osteoclast proliferation by producing macrophage-colony stimulating factor (M-CSF) . Osteoblasts can also produce the receptor activator of nuclear factor kappa-B ligand (RANKL) and osteoprotegerin (OPG) to modulate osteoclast proliferation through the RANKL/RANK/OPG pathway . On the other hand, osteoclasts may also regulate bone formation by osteoblasts . The complex dynamics between the major bone cells control the delicate balance of bone formation/resorption that is critical for maintaining bone health.
While osteoblasts are typically associated with bone remodeling processes, previous studies have demonstrated that they also have the ability to interact with immune cells via the regulation of hematopoietic stem cell (HSC) niches. Specifically, primary osteoblast lineage cells synthesize granulocyte colony-stimulating factor (G-CSF), granulocyte macrophage CSF (GM-CSF), IL-1, IL-6, lymphotoxin, TGFß, TNFα, leukemia inhibitory factor (LIF), and stem cell factor (SCF); all of which play crucial roles in hematopoiesis . Additionally, osteoblasts are responsible for modulating angiogenesis, which supports the high degree of vascularization in the bone marrow environment to provide sufficient oxygen for bone metabolism .
Cellular heterogeneity is an essential characteristic of various cell populations. Although every cell shares almost the same genome, each cell acquires a unique identity and thus specific functional capabilities through molecular coding across the DNA, RNA, and protein conversions . Therefore, even for the same known classical cell types, cells may be further classified into distinct subpopulations due to systematic differences in their gene expression profiles. Emerging evidence from in vitro studies in mice has revealed notable cell-to-cell heterogeneity within the osteoblast cell population. For instance, one study performed the soft agarose cloning technique on rat osteoblastic cells and detected diverse gene expression patterns in osteoblast cells at different stages of cellular differentiation . Another study showed that the expression levels of osteoblastic specific markers including osteopontin, bone sialoprotein, and osteocalcin, varied in the mature osteoblasts of mice with different cellular morphology, suggesting that even these terminally differentiated osteoblasts were composed of multiple subgroups instead of a single unique cell group . While these early efforts revealed the existence of osteoblast heterogeneity, the functional differences between distinct osteoblast subtypes were not well characterized.
The recent development of the state-of-the-art single-cell RNA sequencing (scRNA-seq) technology is expected to provide the most powerful approach to study the nature and characteristics of cell-to-cell heterogeneity. Compared with the conventional bulk RNA-seq approaches, scRNA-seq can reveal complex and rare cell populations, track the trajectories of distinct cell lineages in development, and identify novel regulatory relationships between genes . Accumulating evidence from a few early scRNA-seq studies in vivo in mice has demonstrated the existence of several osteoblast subtypes. Baryawno et al.  identified pre- and mature osteoblast subpopulations based on their transcriptional profiles, while Tikhonova et al.  classified the osteoblasts into three subgroups and revealed significant changes in subgroup proportions under hematopoietic stress conditions induced by chemotherapy treatment. Although our understanding of osteoblast heterogeneity has evolved substantially based on these studies, the cell subtype characteristics of primary osteoblasts in vivo in humans has not been successfully explored. Studying the cellular heterogeneity in mice, or even the cultured cells from humans (though not yet existing), although useful, may not be ideal for studying human disease etiology. This is because the cell identity may vary between mice and humans, and cell culturing may systematically alter the gene expression of the studied cells .
In this study, we successfully performed the first unbiased examination of the in vivo cellular landscape of freshly isolated osteoblasts via scRNA-seq from one 31-year-old male with osteoarthritis and osteopenia. We identified three distinct cell subtypes along with their differentiation relationships based on the transcriptional profiling of 5,329 osteoblast cells from the femur head of one human subject. We then compared the most differentially expressed genes (DEGs) of each cluster with known cell characterizing markers. Further, we identified distinct functional characteristics of each cell subpopulation suggesting their involvement in bone metabolism, angiogenesis modulation, as well as hematopoietic stem cell (HSC) niche regulation. The discovery of osteoblast subtypes is well beyond the scope of current gene expression studies for bone health and represents an important and necessary step to provide novel insights into bone physiological processes at the refined single-cell level.
We applied an established protocol for in vivo human osteoblast isolation  to obtain the alkaline phosphatase (ALPL)high/CD45/34/31low 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 ). After quality control (Supplementary Figure 1B ), we obtained 8,557 cells, with an average of 2,365 and median of 2,260 genes detected per cell (Supplementary Table 1 ). A recently developed dimension reduction technique for scRNA-seq analysis, uniform manifold approximation and projection (UMAP) , was applied to project the gene expression profiles on a two-dimensional panel for visualization of cellular heterogeneity (Supplementary Figure 1C ). Compared with t-SNE , 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 , 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 , 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 ).
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) , 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 ). Based on the identified marker genes and gene ontology (GO) enrichment analysis of the differentially expressed genes in each cluster (Supplementary Figure 1D–1F ), the contaminant cells were classified as; 1) two nucleated erythrocyte groups (C3 and C4, expressing hemoglobin coding genes HBB and HBA1) ; 2) one smooth muscle cell group (C5, expressing transgelin coding gene, TAGLN , and CNN1, which is specific to differentiated mature smooth muscle cells) ; 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 ALPLhigh/CD45/34/31low 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 . 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 ), 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 ).
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  (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)  (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 . We identified three cell subtypes of osteoblasts (Figure 1B–1D) based on their transcriptional profiles and by comparing classical osteoblastic markers  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)  and vascular cell adhesion molecule 1 (encoded by VCAM1) ; 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.  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) , and demonstrated that IGFBP4 could stimulate adult skeletal growth in males . Meanwhile, tenascin, a glycoprotein encoded by TNC that modulates osteoblast mineralization via matrix vesicles , 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 .
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 , 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 ) 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 ). 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 ). 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  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 , 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 . 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 , 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 ). 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 . 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 . 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 .
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 , 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 . It has been reported that osteoblasts lining the bone surface can generate an extracellular environment supporting non-skeletal hematopoietic cells . 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 , 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) . 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 . Since Aguila et al.  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 . 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 . 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) . 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 . 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 . This conflicts with the previous in vitro findings in mice , 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 . 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 . While CPE plays an important role in RANKL-induced osteoclast differentiation , 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  on mouse osteoblasts in vivo (Figure 5A). After correcting for the potential batch effects between independent experiments using canonical correlation analysis (CCA) , 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 . The result showed that human osteoblast cluster O2 overlapped with mouse osteoblast subtype M2, which represents the mature osteoblasts in mice  (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  (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  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.
While it is now well-appreciated that human osteoblast heterogeneity may be distinguished by differentiational stages , the full spectrum of cells that comprise the osteoblast population, especially in vivo in humans, has remained elusive. In this study, for the first time, we classified the freshly isolated primary osteoblasts from human bone (without any in vitro culturing) into three subpopulations based on systematic differences in gene expression profiles and revealed their distinct functional roles in bone metabolism as well as in the regulation of angiogenesis and hemopoiesis. Further, in contrast to the results proposed by Tarkkonen et al. , our results demonstrate systematic differences in the transcription profiles between humans and mice, emphasizing the importance of in vivo studies in humans.
Here, we highlight some of the key findings. First, our results indicate that the osteoblast isolation technique typically used in the field , FACS isolation based on ALPL, was not sufficient for specific isolation of osteoblasts since the isolated cell population included approximately 40% contamination by other cell types. By comparing the transcriptomic profiles between osteoblasts and contamination cells, we hypothesized that the combination of ALPL and COL1A1 would reduce contamination in osteoblast selection. In addition to known osteoblast subtypes, pre- and mature osteoblasts, we also identified one rare osteoblast subpopulation which highly expressed NR4A1 and NR4A2. Based on the gene expression patterns and the inferred osteoblastic lineage trajectory, we found that: 1) Although both pre-osteoblasts (LEPRhigh/VCAM1high) and undetermined osteoblasts (NR4A1high/NR4A2high) are in the differentiation lineage ordering, preosteblasts are primarily responsible for ECM organization during bone formation processes as well as inducing hematopoiesis and modulating angiogenesis, while undetermined osteoblasts are mainly involved in the regulation of osteoblastogenesis and inhibition of hematopoiesis; 2) mature osteoblasts (IBSPhigh/BGALPhigh) arise at the terminal stages of cellular differentiation and are crucial for bone mineralization. We used pre- and mature osteoblasts to annotate clusters O1 and O2, which are vague concepts largely suggested from in vitro literature rather than in vivo studies, mainly based on the expression levels of the gene markers related to specific differentiation stages in the osteoblastic lineage (e.g., LEPR in early stages, BGLAP in terminal stages, etc.). Consistent with the transcriptional characteristics of these two clusters, the result of the pseudotime analysis showed that while cluster O1 cells are enriched in the early stages of the osteoblastic lineage, cluster O2 cells are distributed in the terminal stages.
Despite the novelty of this scRNA-seq study in freshly isolated human osteoblasts, an important limitation is that all the cells were collected from the femur head of one 31-year-old Chinese male subject with osteoarthritis and osteopenia. This might introduce some bias in profiling the gene expression patterns of osteoblasts. However, since no study has demonstrated the disease-associated changes in cell compositions of human osteoblasts at the single-cell resolution, there is no evidence to suggest the existence of the undetermined osteoblasts is due to the bias stemming from the patient’s condition. Further, if the existence of cluster O3 is due to the osteoarthritis, specific markers NR4A1 and NR4A2, would have important roles in osteoarthritis pathogenesis. While previous studies have reported NR4A1  and NR4A2  to be regulators of inflammatory responses in chondrocytes, no study has proposed the relationship between osteoarthritis and NR4A1 produced by osteoblasts. In addition, scRNA-seq studies in both normal mouse osteoblasts in vivo performed by Tikhonova et al.  and transcriptional profiles of mouse osteoblasts in vitro  showed the presence of NR4A1 in normal osteoblasts. Thus, we speculate that cluster O3 may also exist in the normal osteoblastic lineage.
Another limitation is that, due to limited cell numbers, the differentiation relationships between pre-osteoblasts and undetermined osteoblasts have not been clarified here. Additionally, substantial technical differences between the cross-species comparison might exist which might be difficult to be normalized with CCA method, though there is no evidence supporting this. Future studies based on a larger sample size are needed to uncover the disease associated changes in cell subtype compositions, as this will have significant implications for the development of novel therapeutics. Despite these potential limitations, our results provide the first necessary and valuable insights into the cellular heterogeneity of osteoblasts along with a comprehensive and systematic understanding of the cell-intrinsic and cell-specific mechanisms which may underlie bone metabolism and associated disorders.
Materials and Methods
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 . 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.
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 (g/cm2) 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  and the BMD reference established for Chinese population , 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  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) . 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 . 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)  dimension reduction. After data visualization, we applied an unbiased graph-based clustering method  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 . We then used the Monocle 2 v2.8.0  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 . 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. . Next, we integrated the scRNA-seq data of mouse and human osteoblasts by canonical correlation analysis (CCA) using Seurat R package . 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.
Murine mesenchymal stem cells (MSCs) (1.0 × 104 per well) from STEMCELL (Seattle, WA, United States) were plated in 48-well plates and cultured in MesenCultTM 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.
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 iTaqTM 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′,
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.
YG wrote the main manuscript text and conducted major analysis; JY and LC collected the human sample and corresponding clinical information; XL, YC, and CZ performed the experiments; JG, HS, YPW, YPC and HWD did language proofreading; ZW, LJT, and YB prepared supplementary information and validated the results; the study was conceived, designed, initiated, directed, and supervised by HS, HMX, and HWD. All authors participated in the discussions of the project and reviewed and/or revised the manuscript.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
This research was benefited by grants from the National Institutes of Health (R01AR069055, U19AG055373, P20GM109036, R01AG061917), National Natural Science Foundation of China (Grant No. 81902277), National Key R&D Program of China (Grant No. 2017YFC1001100), Natural Science Foundation of Hunan Province (S2019JJQNJJ2093), Changsha Science and Technology project (kq1907153), Central South University (Grant Nos. 164990007, 2018zzts886), and Xiangya Clinical Big Data Project (xyyydsj9).
- 1. Caplan AI. Mesenchymal stem cells. J Orthop Res. 1991; 9:641–50. https://doi.org/10.1002/jor.1100090504 [PubMed]
- 2. Wang T, Zhang X, Bikle DD. Osteogenic Differentiation of Periosteal Cells During Fracture Healing. J Cell Physiol. 2017; 232:913–21. https://doi.org/10.1002/jcp.25641 [PubMed]
- 3. Hallett SA, Ono W, Ono N. Growth Plate Chondrocytes: Skeletal Development, Growth and Beyond. Int J Mol Sci. 2019; 20:6009. https://doi.org/10.3390/ijms20236009 [PubMed]
- 4. Neve A, Corrado A, Cantatore FP. Osteoblast physiology in normal and pathological conditions. Cell Tissue Res. 2011; 343:289–302. https://doi.org/10.1007/s00441-010-1086-1 [PubMed]
- 5. Boskey AL. Biomineralization: Conflicts, challenges, and opportunities. J Cell Biochem. 1998 (Suppl 30-31); 72:83–91. https://doi.org/10.1002/(SICI)1097-4644(1998)72:30/31+<83::AID-JCB12>3.0.CO;2-F [PubMed]
- 6. Rodan GA, Martin TJ. Role of osteoblasts in hormonal control of bone resorption--a hypothesis. Calcif Tissue Int. 1981; 33:349–51. https://doi.org/10.1007/BF02409454 [PubMed]
- 7. Lacey DL, Boyle WJ, Simonet WS, Kostenuik PJ, Dougall WC, Sullivan JK, San Martin J, Dansey R. Bench to bedside: elucidation of the OPG-RANK-RANKL pathway and the development of denosumab. Nat Rev Drug Discov. 2012; 11:401–19. https://doi.org/10.1038/nrd3705 [PubMed]
- 8. Lacey DL, Timms E, Tan HL, Kelley MJ, Dunstan CR, Burgess T, Elliott R, Colombero A, Elliott G, Scully S, Hsu H, Sullivan J, Hawkins N, et al. Osteoprotegerin ligand is a cytokine that regulates osteoclast differentiation and activation. Cell. 1998; 93:165–76. https://doi.org/10.1016/s0092-8674(00)81569-x [PubMed]
- 9. Katagiri T, Takahashi N. Regulatory mechanisms of osteoblast and osteoclast differentiation. Oral Dis. 2002; 8:147–59. https://doi.org/10.1034/j.1601-0825.2002.01829.x [PubMed]
- 10. Neiva K, Sun YX, Taichman RS. The role of osteoblasts in regulating hematopoietic stem cell activity and tumor metastasis. Braz J Med Biol Res. 2005; 38:1449–54. https://doi.org/10.1590/s0100-879x2005001000001 [PubMed]
- 11. Schipani E, Wu C, Rankin EB, Giaccia AJ. Regulation of Bone Marrow Angiogenesis by Osteoblasts during Bone Development and Homeostasis. Front Endocrinol (Lausanne). 2013; 4:85. https://doi.org/10.3389/fendo.2013.00085 [PubMed]
- 12. Altschuler SJ, Wu LF. Cellular heterogeneity: do differences make a difference? Cell. 2010; 141:559–63. https://doi.org/10.1016/j.cell.2010.04.033 [PubMed]
- 13. Guenther HL, Hofstetter W, Stutzer A, Mühlbauer R, Fleisch H. Evidence for heterogeneity of the osteoblastic phenotype determined with clonal rat bone cells established from transforming growth factor-beta-induced cell colonies grown anchorage independently in semisolid medium. Endocrinology. 1989; 125:2092–102. https://doi.org/10.1210/endo-125-4-2092 [PubMed]
- 14. Liu F, Malaval L, Aubin JE. The mature osteoblast phenotype is characterized by extensive plasticity. Exp Cell Res. 1997; 232:97–105. https://doi.org/10.1006/excr.1997.3501 [PubMed]
- 15. Hwang B, Lee JH, Bang D. Single-cell RNA sequencing technologies and bioinformatics pipelines. Exp Mol Med. 2018; 50:1–14. https://doi.org/10.1038/s12276-018-0071-8 [PubMed]
- 16. Baryawno N, Przybylski D, Kowalczyk MS, Kfoury Y, Severe N, Gustafsson K, Kokkaliaris KD, Mercier F, Tabaka M, Hofree M, Dionne D, Papazian A, Lee D, et al. A Cellular Taxonomy of the Bone Marrow Stroma in Homeostasis and Leukemia. Cell. 2019; 177:1915–32.e16. https://doi.org/10.1016/j.cell.2019.04.040 [PubMed]
- 17. Tikhonova AN, Dolgalev I, Hu H, Sivaraj KK, Hoxha E, Cuesta-Domínguez Á, Pinho S, Akhmetzyanova I, Gao J, Witkowski M, Guillamot M, Gutkin MC, Zhang Y, et al. The bone marrow microenvironment at single-cell resolution. Nature. 2019; 569:222–28. https://doi.org/10.1038/s41586-019-1104-8 [PubMed]
- 18. Jilka RL. The relevance of mouse models for investigating age-related bone loss in humans. J Gerontol A Biol Sci Med Sci. 2013; 68:1209–17. https://doi.org/10.1093/gerona/glt046 [PubMed]
- 19. Fujita K, Roforth MM, Atkinson EJ, Peterson JM, Drake MT, McCready LK, Farr JN, Monroe DG, Khosla S. Isolation and characterization of human osteoblasts from needle biopsies without in vitro culture. Osteoporos Int. 2014; 25:887–95. https://doi.org/10.1007/s00198-013-2529-9 [PubMed]
- 20. Roforth MM, Fujita K, McGregor UI, Kirmani S, McCready LK, Peterson JM, Drake MT, Monroe DG, Khosla S. Effects of age on bone mRNA levels of sclerostin and other genes relevant to bone metabolism in humans. Bone. 2014; 59:1–6. https://doi.org/10.1016/j.bone.2013.10.019 [PubMed]
- 21. Fujita K, Roforth MM, Demaray S, McGregor U, Kirmani S, McCready LK, Peterson JM, Drake MT, Monroe DG, Khosla S. Effects of estrogen on bone mRNA levels of sclerostin and other genes relevant to bone metabolism in postmenopausal women. J Clin Endocrinol Metab. 2014; 99:E81–88. https://doi.org/10.1210/jc.2013-3249 [PubMed]
- 22. Gong Y, Yang J, Li X, Zhou C, Chen Y, Wang Z, Qiu X, Liu Y, Zhang H, Greenbaum J, Cheng L, Hu Y, Xie J, et al. A systematic dissection of human primary osteoblasts in vivo at single-cell resolution. bioRxiv. 2020; 2020.05.12.091975. https://doi.org/10.1101/2020.05.12.091975
- 23. Becht E, McInnes L, Healy J, Dutertre CA, Kwok IW, Ng LG, Ginhoux F, Newell EW. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol. 2019; 37:38–44. https://doi.org/10.1038/nbt.4314 [PubMed]
- 24. Kobak D, Berens P. The art of using t-SNE for single-cell transcriptomics. Nat Commun. 2019; 10:5416. https://doi.org/10.1038/s41467-019-13056-x [PubMed]
- 25. Payen VL, Lavergne A, Alevra Sarika N, Colonval M, Karim L, Deckers M, Najimi M, Coppieters W, Charloteaux B, Sokal EM, El Taghdouini A. Single-cell RNA sequencing of human liver reveals hepatic stellate cell heterogeneity. JHEP Rep. 2021; 3:100278. https://doi.org/10.1016/j.jhepr.2021.100278 [PubMed]
- 26. Zou X, Chen K, Zou J, Han P, Hao J, Han Z. Single-cell RNA-seq data analysis on the receptor ACE2 expression reveals the potential risk of different human organs vulnerable to 2019-nCoV infection. Front Med. 2020; 14:185–92. https://doi.org/10.1007/s11684-020-0754-0 [PubMed]
- 27. Yost KE, Satpathy AT, Wells DK, Qi Y, Wang C, Kageyama R, McNamara KL, Granja JM, Sarin KY, Brown RA, Gupta RK, Curtis C, Bucktrout SL, et al. Clonal replacement of tumor-specific T cells following PD-1 blockade. Nat Med. 2019; 25:1251–59. https://doi.org/10.1038/s41591-019-0522-3 [PubMed]
- 28. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018; 36:411–20. https://doi.org/10.1038/nbt.4096 [PubMed]
- 29. Meyer MB, Benkusky NA, Pike JW. The RUNX2 cistrome in osteoblasts: characterization, down-regulation following differentiation, and relationship to gene expression. J Biol Chem. 2014; 289:16016–31. https://doi.org/10.1074/jbc.M114.552216 [PubMed]
- 30. Capellera-Garcia S, Pulecio J, Dhulipala K, Siva K, Rayon-Estrada V, Singbrant S, Sommarin MN, Walkley CR, Soneji S, Karlsson G, Raya Á, Sankaran VG, Flygare J. Defining the Minimal Factors Required for Erythropoiesis through Direct Lineage Conversion. Cell Rep. 2016; 15:2550–62. https://doi.org/10.1016/j.celrep.2016.05.027 [PubMed]
- 31. Tawfik O, Rao D, Nothnick WB, Graham A, Mau B, Fan F. Transgelin, a Novel Marker of Smooth Muscle Differentiation, Effectively Distinguishes Endometrial Stromal Tumors from Uterine Smooth Muscle Tumors. Int J Gynecol Obstet Reprod Med Res. 2014; 1:26–31. [PubMed]
- 32. Liu R, Jin JP. Calponin isoforms CNN1, CNN2 and CNN3: Regulators for actin cytoskeleton functions in smooth muscle and non-muscle cells. Gene. 2016; 585:143–53. https://doi.org/10.1016/j.gene.2016.02.040 [PubMed]
- 33. Francis A, Bosio E, Stone SF, Fatovich DM, Arendts G, MacDonald SP, Burrows S, Brown SG. Markers Involved in Innate Immunity and Neutrophil Activation are Elevated during Acute Human Anaphylaxis: Validation of a Microarray Study. J Innate Immun. 2019; 11:63–73. https://doi.org/10.1159/000492301 [PubMed]
- 34. Tidwell T, Wechsler J, Nayak RC, Trump L, Salipante SJ, Cheng JC, Donadieu J, Glaubach T, Corey SJ, Grimes HL, Lutzko C, Cancelas JA, Horwitz MS. Neutropenia-associated ELANE mutations disrupting translation initiation produce novel neutrophil elastase isoforms. Blood. 2014; 123:562–69. https://doi.org/10.1182/blood-2013-07-513242 [PubMed]
- 35. Limas CJ, Cohn JN. Alkaline phosphatase in vascular smooth muscle. Nat New Biol. 1973; 245:53–55. https://doi.org/10.1038/newbio245053a0 [PubMed]
- 36. Liu B, Lu Y, Wang Y, Ge L, Zhai N, Han J. A protocol for isolation and identification and comparative characterization of primary osteoblasts from mouse and rat calvaria. Cell Tissue Bank. 2019; 20:173–82. https://doi.org/10.1007/s10561-019-09751-0 [PubMed]
- 37. Yu W, Zhong L, Yao L, Wei Y, Gui T, Li Z, Kim H, Holdreith N, Jiang X, Tong W, Dyment N, Liu XS, Yang S, et al. Bone marrow adipogenic lineage precursors promote osteoclastogenesis in bone remodeling and pathologic bone loss. J Clin Invest. 2021; 131:e140214. https://doi.org/10.1172/JCI140214 [PubMed]
- 38. Zhong L, Yao L, Tower RJ, Wei Y, Miao Z, Park J, Shrestha R, Wang L, Yu W, Holdreith N, Huang X, Zhang Y, Tong W, et al. Single cell transcriptomics identifies a unique adipose lineage cell population that regulates bone marrow environment. Elife. 2020; 9:e54695. https://doi.org/10.7554/eLife.54695 [PubMed]
- 39. Baccin C, Al-Sabah J, Velten L, Helbling PM, Grünschläger F, Hernández-Malmierca P, Nombela-Arrieta C, Steinmetz LM, Trumpp A, Haas S. Combined single-cell and spatial transcriptomics reveal the molecular, cellular and spatial bone marrow niche organization. Nat Cell Biol. 2020; 22:38–48. https://doi.org/10.1038/s41556-019-0439-6 [PubMed]
- 40. Serowoky MA, Arata CE, Crump JG, Mariani FV. Skeletal stem cells: insights into maintaining and regenerating the skeleton. Development. 2020; 147:dev179325. https://doi.org/10.1242/dev.179325 [PubMed]
- 41. Matsushita Y, Nagata M, Kozloff KM, Welch JD, Mizuhashi K, Tokavanich N, Hallett SA, Link DC, Nagasawa T, Ono W, Ono N. A Wnt-mediated transformation of the bone marrow stromal cell identity orchestrates skeletal regeneration. Nat Commun. 2020; 11:332. https://doi.org/10.1038/s41467-019-14029-w [PubMed]
- 42. Rutkovskiy A, Stensløkken KO, Vaage IJ. Osteoblast Differentiation at a Glance. Med Sci Monit Basic Res. 2016; 22:95–106. https://doi.org/10.12659/msmbr.901142 [PubMed]
- 43. Zhou BO, Yue R, Murphy MM, Peyer JG, Morrison SJ. Leptin-receptor-expressing mesenchymal stromal cells represent the main source of bone formed by adult bone marrow. Cell Stem Cell. 2014; 15:154–68. https://doi.org/10.1016/j.stem.2014.06.008 [PubMed]
- 44. Akiyama K, You YO, Yamaza T, Chen C, Tang L, Jin Y, Chen XD, Gronthos S, Shi S. Characterization of bone marrow derived mesenchymal stem cells in suspension. Stem Cell Res Ther. 2012; 3:40. https://doi.org/10.1186/scrt131 [PubMed]
- 45. Hanagata N, Li X, Morita H, Takemura T, Li J, Minowa T. Characterization of the osteoblast-specific transmembrane protein IFITM5 and analysis of IFITM5-deficient mice. J Bone Miner Metab. 2011; 29:279–90. https://doi.org/10.1007/s00774-010-0221-0 [PubMed]
- 46. Komori T. Regulation of bone development and extracellular matrix protein genes by RUNX2. Cell Tissue Res. 2010; 339:189–95. https://doi.org/10.1007/s00441-009-0832-8 [PubMed]
- 47. Kwan Tat S, Pelletier JP, Lajeunesse D, Fahmi H, Lavigne M, Martel-Pelletier J. The differential expression of osteoprotegerin (OPG) and receptor activator of nuclear factor kappaB ligand (RANKL) in human osteoarthritic subchondral bone osteoblasts is an indicator of the metabolic state of these disease cells. Clin Exp Rheumatol. 2008; 26:295–304. [PubMed]
- 48. Akeson G, Malemud CJ. A Role for Soluble IL-6 Receptor in Osteoarthritis. J Funct Morphol Kinesiol. 2017; 2:27. https://doi.org/10.3390/jfmk2030027 [PubMed]
- 49. Sakao K, Takahashi KA, Arai Y, Saito M, Honjo K, Hiraoka N, Asada H, Shin-Ya M, Imanishi J, Mazda O, Kubo T. Osteoblasts derived from osteophytes produce interleukin-6, interleukin-8, and matrix metalloproteinase-13 in osteoarthritis. J Bone Miner Metab. 2009; 27:412–23. https://doi.org/10.1007/s00774-009-0058-6 [PubMed]
- 50. Xi G, Rosen CJ, Clemmons DR. IGF-I and IGFBP-2 Stimulate AMPK Activation and Autophagy, Which Are Required for Osteoblast Differentiation. Endocrinology. 2016; 157:268–81. https://doi.org/10.1210/en.2015-1690 [PubMed]
- 51. Maridas DE, DeMambro VE, Le PT, Nagano K, Baron R, Mohan S, Rosen CJ. IGFBP-4 regulates adult skeletal growth in a sex-specific manner. J Endocrinol. 2017; 233:131–44. https://doi.org/10.1530/JOE-16-0673 [PubMed]
- 52. Li C, Cui Y, Luan J, Zhou X, Li H, Wang H, Shi L, Han J. Tenascin C affects mineralization of SaOS2 osteoblast-like cells through matrix vesicles. Drug Discov Ther. 2016; 10:82–87. https://doi.org/10.5582/ddt.2016.01009 [PubMed]
- 53. Guntur AR, Rosen CJ, Naski MC. N-cadherin adherens junctions mediate osteogenesis through PI3K signaling. Bone. 2012; 50:54–62. https://doi.org/10.1016/j.bone.2011.09.036 [PubMed]
- 54. Tetradis S, Bezouglaia O, Tsingotjidou A, Vila A. Regulation of the nuclear orphan receptor Nur77 in bone by parathyroid hormone. Biochem Biophys Res Commun. 2001; 281:913–16. https://doi.org/10.1006/bbrc.2001.4459 [PubMed]
- 55. Lammi J, Huppunen J, Aarnisalo P. Regulation of the osteopontin gene by the orphan nuclear receptor NURR1 in osteoblasts. Mol Endocrinol. 2004; 18:1546–57. https://doi.org/10.1210/me.2003-0247 [PubMed]
- 56. Lee MK, Choi H, Gil M, Nikodem VM. Regulation of osteoblast differentiation by Nurr1 in MC3T3-E1 cell line and mouse calvarial osteoblasts. J Cell Biochem. 2006; 99:986–94. https://doi.org/10.1002/jcb.20990 [PubMed]
- 57. Pirih FQ, Tang A, Ozkurt IC, Nervina JM, Tetradis S. Nuclear orphan receptor Nurr1 directly transactivates the osteocalcin gene in osteoblasts. J Biol Chem. 2004; 279:53167–74. https://doi.org/10.1074/jbc.M405677200 [PubMed]
- 58. Zeitouni S, Ford BS, Harris SM, Whitney MJ, Gregory CA, Prockop DJ. Pharmaceutical induction of ApoE secretion by multipotent mesenchymal stromal cells (MSCs). BMC Biotechnol. 2008; 8:75. https://doi.org/10.1186/1472-6750-8-75 [PubMed]
- 59. Rosset EM, Bradshaw AD. SPARC/osteonectin in mineralized tissue. Matrix Biol. 2016; 52:78–87. https://doi.org/10.1016/j.matbio.2016.02.001 [PubMed]
- 60. Dufourcq P, Descamps B, Tojais NF, Leroux L, Oses P, Daret D, Moreau C, Lamazière JM, Couffinhal T, Duplàa C. Secreted frizzled-related protein-1 enhances mesenchymal stem cell function in angiogenesis and contributes to neovessel maturation. Stem Cells. 2008; 26:2991–3001. https://doi.org/10.1634/stemcells.2008-0372 [PubMed]
- 61. Choudhuri R, Zhang HT, Donnini S, Ziche M, Bicknell R. An angiogenic role for the neurokines midkine and pleiotrophin in tumorigenesis. Cancer Res. 1997; 57:1814–19. [PubMed]
- 62. Lawler JW, Slayter HS, Coligan JE. Isolation and characterization of a high molecular weight glycoprotein from human blood platelets. J Biol Chem. 1978; 253:8609–16. https://doi.org/10.1016/S0021-9258(17)34336-3 [PubMed]
- 63. Karar J, Maity A. PI3K/AKT/mTOR Pathway in Angiogenesis. Front Mol Neurosci. 2011; 4:51. https://doi.org/10.3389/fnmol.2011.00051 [PubMed]
- 64. Kang DW, Choi KY, Min DS. Functional regulation of phospholipase D expression in cancer and inflammation. J Biol Chem. 2014; 289:22575–82. https://doi.org/10.1074/jbc.R114.569822 [PubMed]
- 65. Sarwar M, Du XJ, Dschietzig TB, Summers RJ. The actions of relaxin on the human cardiovascular system. Br J Pharmacol. 2017; 174:933–49. https://doi.org/10.1111/bph.13523 [PubMed]
- 66. Morrison SJ, Scadden DT. The bone marrow niche for haematopoietic stem cells. Nature. 2014; 505:327–34. https://doi.org/10.1038/nature12984 [PubMed]
- 67. Zumkeller W, Burdach S. The insulin-like growth factor system in normal and malignant hematopoietic cells. Blood. 1999; 94:3653–57. https://doi.org/10.1182/blood.V94.11.3653 [PubMed]
- 68. Cecchini MG, Hofstetter W, Halasy J, Wetterwald A, Felix R. Role of CSF-1 in bone and bone marrow development. Mol Reprod Dev. 1997; 46:75–83. https://doi.org/10.1002/(SICI)1098-2795(199701)46:1<75::AID-MRD12>3.0.CO;2-2 [PubMed]
- 69. Zhang CC, Lodish HF. Cytokines regulating hematopoietic stem cell function. Curr Opin Hematol. 2008; 15:307–11. https://doi.org/10.1097/MOH.0b013e3283007db5 [PubMed]
- 70. Aguila HL, Mun SH, Kalinowski J, Adams DJ, Lorenzo JA, Lee SK. Osteoblast-specific overexpression of human interleukin-7 rescues the bone mass phenotype of interleukin-7-deficient female mice. J Bone Miner Res. 2012; 27:1030–42. https://doi.org/10.1002/jbmr.1553 [PubMed]
- 71. Chen CY, Fuh LJ, Huang CC, Hsu CJ, Su CM, Liu SC, Lin YM, Tang CH. Enhancement of CCL2 expression and monocyte migration by CCN1 in osteoblasts through inhibiting miR-518a-5p: implication of rheumatoid arthritis therapy. Sci Rep. 2017; 7:421. https://doi.org/10.1038/s41598-017-00513-0 [PubMed]
- 72. Gagnon J, Ramanathan S, Leblanc C, Cloutier A, McDonald PP, Ilangumaran S. IL-6, in synergy with IL-7 or IL-15, stimulates TCR-independent proliferation and functional differentiation of CD8+ T lymphocytes. J Immunol. 2008; 180:7958–68. https://doi.org/10.4049/jimmunol.180.12.7958 [PubMed]
- 73. Park IK, Trotta R, Yu J, Caligiuri MA. Axl/Gas6 pathway positively regulates FLT3 activation in human natural killer cell development. Eur J Immunol. 2013; 43:2750–55. https://doi.org/10.1002/eji.201243116 [PubMed]
- 74. Vignudelli T, Selmi T, Martello A, Parenti S, Grande A, Gemelli C, Zanocco-Marani T, Ferrari S. ZFP36L1 negatively regulates erythroid differentiation of CD34+ hematopoietic stem cells by interfering with the Stat5b pathway. Mol Biol Cell. 2010; 21:3340–51. https://doi.org/10.1091/mbc.E10-01-0040 [PubMed]
- 75. Doumas S, Kolokotronis A, Stefanopoulos P. Anti-inflammatory and antimicrobial roles of secretory leukocyte protease inhibitor. Infect Immun. 2005; 73:1271–74. https://doi.org/10.1128/IAI.73.3.1271-1274.2005 [PubMed]
- 76. Su AI, Cooke MP, Ching KA, Hakak Y, Walker JR, Wiltshire T, Orth AP, Vega RG, Sapinoso LM, Moqrich A, Patapoutian A, Hampton GM, Schultz PG, Hogenesch JB. Large-scale analysis of the human and mouse transcriptomes. Proc Natl Acad Sci USA. 2002; 99:4465–70. https://doi.org/10.1073/pnas.012025199 [PubMed]
- 77. Cremasco V, Decker CE, Stumpo D, Blackshear PJ, Nakayama KI, Nakayama K, Lupu TS, Graham DB, Novack DV, Faccio R. Protein kinase C-delta deficiency perturbs bone homeostasis by selective uncoupling of cathepsin K secretion and ruffled border formation in osteoclasts. J Bone Miner Res. 2012; 27:2452–63. https://doi.org/10.1002/jbmr.1701 [PubMed]
- 78. Kim HJ, Hong J, Yoon HJ, Yoon YR, Kim SY. Carboxypeptidase E is a novel modulator of RANKL-induced osteoclast differentiation. Mol Cells. 2014; 37:685–90. https://doi.org/10.14348/molcells.2014.0179 [PubMed]
- 79. Tarkkonen K, Hieta R, Kytölä V, Nykter M, Kiviranta R. Comparative analysis of osteoblast gene expression profiles and Runx2 genomic occupancy of mouse and human osteoblasts in vitro. Gene. 2017; 626:119–31. https://doi.org/10.1016/j.gene.2017.05.028 [PubMed]
- 80. Xiong Y, Ran J, Xu L, Tong Z, Adel Abdo MS, Ma C, Xu K, He Y, Wu Z, Chen Z, Hu P, Jiang L, Bao J, et al. Reactivation of NR4A1 Restrains Chondrocyte Inflammation and Ameliorates Osteoarthritis in Rats. Front Cell Dev Biol. 2020; 8:158. https://doi.org/10.3389/fcell.2020.00158 [PubMed]
- 81. McCoy JM, Walkenhorst DE, McCauley KS, Elaasar H, Everett JR, Mix KS. Orphan nuclear receptor NR4A2 induces transcription of the immunomodulatory peptide hormone prolactin. J Inflamm (Lond). 2015; 12:13. https://doi.org/10.1186/s12950-015-0059-2 [PubMed]
- 82. Xie H, Sun M, Liao XB, Yuan LQ, Sheng ZF, Meng JC, Wang D, Yu ZY, Zhang LY, Zhou HD, Luo XH, Li H, Wu XP, et al. Estrogen receptor α36 mediates a bone-sparing effect of 17β-estrodiol in postmenopausal women. J Bone Miner Res. 2011; 26:156–68. https://doi.org/10.1002/jbmr.169 [PubMed]
- 83. Deng FY, Lei SF, Li MX, Jiang C, Dvornyk V, Deng HW. Genetic determination and correlation of body mass index and bone mineral density at the spine and hip in Chinese Han ethnicity. Osteoporos Int. 2006; 17:119–24. https://doi.org/10.1007/s00198-005-1930-4 [PubMed]
- 84. Liang X, Wu C, Zhao H, Liu L, Du Y, Li P, Wen Y, Zhao Y, Ding M, Cheng B, Cheng S, Ma M, Zhang L, et al. Assessing the genetic correlations between early growth parameters and bone mineral density: A polygenic risk score analysis. Bone. 2018; 116:301–06. https://doi.org/10.1016/j.bone.2018.08.021 [PubMed]
- 85. Kanis JA, Melton LJ 3rd, Christiansen C, Johnston CC, Khaltaev N. The diagnosis of osteoporosis. J Bone Miner Res. 1994; 9:1137–41. https://doi.org/10.1002/jbmr.5650090802 [PubMed]
- 86. Wu XP, Liao EY, Zhang H, Shan PF, Cao XZ, Liu SP. Establishment of BMD reference plots and determination of peak BMD at multiple skeletal regions in mainland Chinese women and the diagnosis of osteoporosis. Osteoporos Int. 2004; 15:71–79. https://doi.org/10.1007/s00198-003-1517-x [PubMed]
- 87. Wu CL, Dicks A, Steward N, Tang R, Katz DB, Choi YR, Guilak F. Single cell transcriptomic analysis of human pluripotent stem cell chondrogenesis. Nat Commun. 2021; 12:362. https://doi.org/10.1038/s41467-020-20598-y [PubMed]
- 88. De Micheli AJ, Spector JA, Elemento O, Cosgrove BD. A reference single-cell transcriptomic atlas of human skeletal muscle tissue reveals bifurcated muscle stem cell populations. Skelet Muscle. 2020; 10:19. https://doi.org/10.1186/s13395-020-00236-3 [PubMed]
- 89. Hildreth AD, Ma F, Wong YY, Sun R, Pellegrini M, O’Sullivan TE. Single-cell sequencing of human white adipose tissue identifies new cell states in health and obesity. Nat Immunol. 2021; 22:639–53. https://doi.org/10.1038/s41590-021-00922-4 [PubMed]
- 90. Blondel VD, Guillaume JL, Lambiotte R, Lefebvre E. Fast unfolding of communities in large networks. J Stat Mech. 2008; 2008:P10008. https://doi.org/10.1088/1742-5468/2008/10/P10008
- 91. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–87. https://doi.org/10.1089/omi.2011.0118 [PubMed]
- 92. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, Lennon NJ, Livak KJ, Mikkelsen TS, Rinn JL. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014; 32:381–86. https://doi.org/10.1038/nbt.2859 [PubMed]
- 93. Stuart T, Satija R. Integrative single-cell analysis. Nat Rev Genet. 2019; 20:257–72. https://doi.org/10.1038/s41576-019-0093-7 [PubMed]