Dissecting human embryonic skeletal stem cell ontogeny by single-cell transcriptomic and functional analyses

Abstract

Human skeletal stem cells (SSCs) have been discovered in fetal and adult long bones. However, the spatiotemporal ontogeny of human embryonic SSCs during early skeletogenesis remains elusive. Here we map the transcriptional landscape of human limb buds and embryonic long bones at single-cell resolution to address this fundamental question. We found remarkable heterogeneity within human limb bud mesenchyme and epithelium, and aligned them along the proximal-distal and anterior-posterior axes using known marker genes. Osteo-chondrogenic progenitors first appeared in the core limb bud mesenchyme, which give rise to multiple populations of stem/progenitor cells in embryonic long bones undergoing endochondral ossification. Importantly, a perichondrial embryonic skeletal stem/progenitor cell (eSSPC) subset was identified, which could self-renew and generate the osteochondral lineage cells, but not adipocytes or hematopoietic stroma. eSSPCs are marked by the adhesion molecule CADM1 and highly enriched with FOXP1/2 transcriptional network. Interestingly, neural crest-derived cells with similar phenotypic markers and transcriptional networks were also found in the sagittal suture of human embryonic calvaria. Taken together, this study revealed the cellular heterogeneity and lineage hierarchy during human embryonic skeletogenesis, and identified distinct skeletal stem/progenitor cells that orchestrate endochondral and intramembranous ossification.

Fig. 1. Integrated analysis of human limb buds and embryonic long bones.

a Representative stereoscope images (left) and H&E images (right) of 5 WPC human limb bud and 8 WPC human long bone (n = 2 embryos). Scale bars, 100 μm. b Sampling workflow and experimental scheme. Human embryonic cells from 5 WPC limb buds and 8 WPC long bones were sorted and subjected to droplet-based scRNA-seq. c Distribution of 35,570 cells from limb buds and long bones. In total, 16 subsets were visualized by uniform manifold approximation and projection (UMAP). d Dot plots showing the expression of curated feature genes in 16 subsets. Dot size represented the proportion of cells expressing specific gene in the indicated subset and color bar represented the gene expression levels. e Proportion of cells from 5 WPC limb buds and 8 WPC long bones in each subset. f Developmental trajectory inferred by RNA velocity and visualized on the UMAP projection. g Partition-based graph abstraction (PAGA) showing the connectivity among subsets in f. The mean expression of representative genes (Mesenchymal, PRRX1; Chondrogenic, SOX9; Osteogenic, RUNX2) in each subset was shown in abstracted graph. Line thickness indicated the strength of connectivity. Color bar represents the gene expression levels.

Fig. 2. Characterization of human limb bud mesenchyme and epithelium.

a UAMP visualization of the ten subsets in 5 WPC limb buds. b Hierarchical clustering of the mesenchymal and epithelial subsets using top 50 principal components (PCs). c The inferred relationships among the mesenchymal and epithelial subsets in PAGA layout. d Stacked bar charts showing the cell cycle distributions in the mesenchymal subsets. e Enriched GO terms of differentially expressed genes (DEGs) in the mesenchymal subsets. f Heatmap showing expression of curated HOX genes scaled across the mesenchymal subsets. Hox genes were clustered into two branches based on hierarchical clustering of the rows, as indicated in green and purple. g Visualization of the mesenchymal subsets (left) with UMAP plots showing the expression of curated PD and AP marker genes (right; Proximal, MEIS2; Distal, HOXD13; Anterior, IRX3; Posterior, SHH). h GSVA analysis of pathway enrichment in the proximal and core mesenchyme (LBM3/OCP) and distal mesenchyme (LBM1/2). T values for each pathway were shown (two-sided unpaired limma-moderated t-test). i Heatmap showing the area under the curve (AUC) score of regulons enriched in the mesenchymal subsets (left). Z-score (row scaling) was computed. Representative regulons were shown on the right. The number of predicted target genes for each regulon was shown in the parenthesis. Hierarchical clustering on columns indicated correlation between cell subsets. Binary activities of representative regulons were shown by UMAP plots (right).

Fig. 3. Characterization of the osteochondral lineage in human long bones identified embryonic SSCs.

a UMAP visualization of seven OCLC subsets in 8 WPC human long bones. b Enriched GO terms of DEGs among the seven OCLC subsets. c Developmental trajectory of seven OCLC subsets inferred by RNA velocity and visualized on the UMAP projection. d Diffusion map visualization of the osteogenic and chondrogenic trajectories simulated by Slingshot across LBDMC, eSSPC, osteoprogenitor, chondroblast and chondrocyte subsets. The corresponding diffusion pseudotime was indicated in the upper right frame. e Heatmap of gene expressions (smoothed over 20 adjacent cells) in LBDMC, eSSPC, osteoprogenitor, chondroblast and chondrocyte subsets ordered by pseudotime of osteogenesis and chondrogenesis in d. Top 200 genes were selected according to the P values of GVM test and representative genes were shown. Shared genes in the two trajectories were indicated in dashed box. f Heatmap showing the AUC score of regulons enriched in human OCLC subsets. Z-score (row scaling) was calculated. Representative regulons were shown on the right. The number of predicted target genes for each regulon was shown in the parenthesis. g Binary activities of FOXP1 and FOXP2 regulons were shown by UMAP plots. h The FOXP1 and FOXP2 regulon networks in OCLC subsets. Line thickness indicated the level of GENIE3 weights. Dot size indicated the number of enriched TF motifs. i, j Immunofluorescent images of FOXP1 (i) and FOXP2 (j) expression in 8 WPC human femur. FOXP1/2⁺ cells were detected in the perichondrium (I) and inside POC (II). Merged and single-channel images of FOXP1/2 (red) and DAPI (blue) were shown (n = 2 embryos). Scale bars in snapshot images, 200 μm; scale bars in magnified images, 50 μm.

Fig. 4. Identification of CADM1 as a phenotypic marker of eSSPCs.

a Dot plots showing the expression of differentially expressed cell surface genes (left) and candidate SSC markers (right) in 8 WPC human long bone subsets. Asterisks indicated positive markers that were used to enrich eSSPCs. b Immunofluorescent images of PDPN⁺CADM1⁺ cells in 8 WPC human long bones. Overviews of PDPN⁺CADM1⁺ cells (arrows) in the articular (upper left) and POC (bottom left) regions were shown on the left. PDPN⁺CADM1⁺ cells were found in the inner layer of perichondrium in the articular regions (I) and surrounding POC (II). A few PDPN⁺CADM1⁺ cells were also found inside POC (III). Arrow heads indicated enlarged PDPN⁺CADM1⁺ cells. Merged and single-channel images of DAPI (blue), CADM1 (red) and PDPN (green) were shown (n = 2 embryos). Scale bars in snapshot images, 50 μm; scale bars in magnified images, 5 μm. c Flow cytometry gating strategies for sorting different populations in 8 WPC long bones (n = 3 embryos). d Representative crystal violet staining of CFU-F colonies generated by the sorted populations as indicated in c. Magnified images of the boxed areas were shown on the right. Scale bars, 25 μm. e Quantifications of the number (left) and mean diameter (right) of the CFU-F colonies (n = 3 embryos). The statistical significance of differences was determined using one-way ANOVA with multiple comparison tests (LSD). *P < 0.05; **P < 0.01; ***P < 0.001. Error bars indicated SEM.

Fig. 5. Functional characterizations of eSSPCs in vitro and in vivo.

a Flow cytometry plots showing the maintenance of phenotypic eSSPCs after serially passaging clonally expanded PDGFRA^low/–PDPN⁺CADM1⁺ cells (n = 3 clones). b Representative oil red O (top), alizarin red (middle) and toluidine blue (bottom) staining after adipogenic, osteogenic and chondrogenic differentiation of clonally expanded eSSPCs (PDGFRA^low/–PDPN⁺CADM1⁺). Magnified images of the boxed areas were shown on the right. Scale bars, 200 μm. c qPCR analyses of adipogenic, osteogenic and chondrogenic marker genes in clonally expanded eSSPCs before and after trilineage differentiation in vitro (n = 3 clones). The statistical significance of differences was determined using Wilcoxon signed rank test. *P < 0.05; **P < 0.01. Error bars indicated SEM. d Renal subcapsular transplantation. The work flow for functional characterization of eSSPC in vivo (top). Subcapsular xenografts were dissected and sectioned 8 weeks after transplantation of culture expanded eSSPCs into immunodeficient mice. Bright field (middle), Movat pentachrome staining (bottom left; cartilage, blue; bone and fibrous tissue, yellow) and immunofluorescent staining images (bottom right; DAPI, blue; collagen I (COL1), red; collagen II (COL2), green) were shown (n = 9 grafts from three embryos). Scale bars, 50 μm.

Fig. 6. Characterization of the osteogenic lineages in human embryonic calvaria identified neural crest-derived skeletal progenitors.

a UMAP visualization of 12 subsets in 8 WPC calvarial bones (n = 2 embryos). Inset illustrated the position of calvarial bone. b Violin plots showing the expression of feature genes for each subset. c Heatmap showing the transcriptome correlation between osteogenic subsets in calvarial and OCLC subsets in long bone. Asterisks indicated subsets with correlation coefficients > 0.8. d Dot plots (left) and UMAP plots (right) showing the expression of eSSPC marker genes in 12 subsets of 8 WPC calvarial. e UMAP visualization of the two osteogenic trajectories simulated by Slingshot across NC, mig_NC, NCDC, osteoprogenitor, PMSC1 and PMSC2 subsets (Upper left). Expression UMAP plots of marker genes (NC, FOXC1; Mesoderm, TWIST2; Osteoprogenitor, DLX5). f Heatmap of the gene expressions (smoothed over 20 adjacent cells) in subsets ordered by pseudotime of osteogenesis as in e. Top 200 genes were selected according to the P values of GVM test and representative genes were shown. Shared genes in two trajectories were indicated in dashed box. g Heatmap showing the AUC scores of regulons enriched in the osteogenic subsets. Z-score (row scaling) was computed. Representative regulons were shown on the right. h Binary activities of FOXP1/2/4 regulons were shown by UMAP plots.