Interspecies Single-Cell RNA-Seq Analysis Reveals the Novel Trajectory of Osteoclast Differentiation and Therapeutic Targets

Abstract

Bone turnover is finely tuned by cells in the bone milieu, including osteoblasts, osteoclasts, and osteocytes. Osteoclasts are multinucleated giant cells with a bone-resorbing function that play a critical role in regulating skeletal homeostasis. Osteoclast differentiation is characterized by dramatic changes in morphology and gene expression following receptor activator of nuclear factor-kappa-Β ligand (RANKL) stimulation. We performed single-cell RNA-sequencing analyses of human and murine osteoclast-lineage cells (OLCs) and found that OLCs in the mitotic phase do not differentiate into mature osteoclasts. We also identified a guanosine triphosphatase (GTPase) family member, RAB38, as a highly expressed molecule in both human and murine osteoclast clusters; RAB38 gene expression is associated with dynamic changes in histone modification and transcriptional regulation. Silencing Rab38 expression by using short hairpin RNA (shRNA) inhibited osteoclast differentiation and maturation. In summary, we established an integrated fate map of human and murine osteoclastogenesis; this will help identify therapeutic targets in bone diseases. © 2022 The Authors. JBMR Plus published by Wiley Periodicals LLC on behalf of American Society for Bone and Mineral Research.

Keywords: ChIP‐SEQUENCING; HISTONE MODIFICATION; OSTEOCLAST; RAB38; SINGLE‐CELL RNA‐SEQUENCING.

Fig. 1

In vitro human and murine osteoclast differentiation following stimulation with recombinant M‐CSF and RANKL. TRAP staining in mature human (A) and murine (B) osteoclasts (scale bar = 200 μm). The expression of differentiation factors in human (C) and murine (D) osteoclasts. The values were measured in triplicate and compared against day 0 (***p < 0.001). M‐CSF = macrophage colony‐stimulating factor; RANKL = receptor activator of nuclear factor‐kappa‐Β ligand; TRAP = tartrate‐resistant acid phosphatase.

Fig. 2

Transcriptional profiling of human osteoclastogenesis by scRNA‐seq. In vitro cultured human mature osteoclasts generated by stimulation with recombinant M‐CSF and RANKL were analyzed by scRNA‐seq. (A) Feature plot (UMAP) demonstrating a total of 11 clusters during human osteoclastogenesis (n = 4707 cells). (B) Cluster annotation and cell number of each cluster in the human OLC data set. (C) Cell‐based trajectory analysis in human osteoclastogenesis using Monocle 3 in UMAP. (D) Phase analysis of cell‐cycle status (red, G1; green, G2/M; blue, S phase). (E) Heat map of the top 25 differentially expressed genes in each OLC cluster. Osteoclastic markers CTSK and ATP6V0D2 are highly expressed in cluster 6. (F,G) Feature plot (F) and violin plot (G) of marker genes. CD14 and CXCR4 are progenitor markers. TNFRSF11A, NFATc1, CTSK, and ATP6V0D2 are osteoclast‐associated genes. (H) Dot plot of GO terms enriched in each cluster related to “dendritic cell,” “macrophage,” “osteoclast,” and “bone.” GO = gene ontology; M‐CSF = macrophage colony‐stimulating factor; RANKL = receptor activator of nuclear factor‐kappa‐Β ligand.

Fig. 3

Transcriptional profiling of murine osteoclastogenesis by scRNA‐seq. In vitro cultured murine mature osteoclasts generated by stimulation with recombinant M‐CSF and RANKL were analyzed by scRNA‐seq. (A) Cell‐based trajectory analysis in murine osteoclastogenesis. (B) Cluster annotation and cell number of each cluster in the murine OLC data set. (C) Pseudo‐time analysis shows a cell‐based trajectory in murine osteoclastogenesis. (D) Phase analysis of cell‐cycle status (red, G1; green, G2/M; blue, S phase). (E) Heat map of the top 25 differentially expressed genes in each OLC cluster. Osteoclastic markers Acp5, Ocstamp, Ctsk, and Atp6v0d2 are highly expressed in cluster 6. (F,G) Feature plot (F) and violin plot (G) of marker genes. Cx3cr1 and Ccr2 are progenitor markers. Tnfrsf11a, Nfatc1, Ctsk, and Atp6v0d2 are osteoclast‐associated genes. (H) Dot plot of gene ontology terms enriched in each cluster related to “dendritic cell,” “macrophage,” “osteoclast,” and “bone.” M‐CSF = macrophage colony‐stimulating factor; RANKL = receptor activator of nuclear factor‐kappa‐Β ligand.

Fig. 4

Analysis of integrated data from human and humanized murine osteoclast differentiation. Human and humanized mouse scRNA‐seq data were integrated. (A) Feature plot (UMAP) demonstrating a total of 13 clusters during osteoclastogenesis (human, n = 4285; humanized mouse, n = 4417; total, n = 8702 cells). Cell‐based trajectory analysis of the integrated data set. Cluster annotation and cell number of each cluster in the human and humanized murine OLC data set. (B) Feature plots and cell‐based trajectory in the reconstructed human and murine OLC data set. Arrows indicate the branch into OLC. (C) Phase analysis of cell‐cycle status in the reconstructed human and murine OLC data set (red, G1; green, G2/M; blue, S phase). (D) Heat map of the top 25 differentially expressed genes in each OLC cluster. Osteoclastic markers ACP5, CTSK, ATP6V0D2, OCSTAMP, DCSTAMP, OSCAR, and NFATC1 are highly expressed in clusters 5 and 6. (E) Feature plot of marker genes. CD14 and CXCR4 are progenitor markers. TNFRSF11A, NFATc1, CTSK, and ATP6V0D2 are osteoclast‐associated genes. (F) Dot plot of GO terms enriched in each cluster related to “dendritic cell,” “macrophage,” “osteoclast,” and “bone.” (G) Expression of detectable RAB family genes in the integrated data set. Note that RAB38 is predominantly expressed in osteoclastic clusters 5 and 6.

Fig. 5

Analysis of integrated data from murinized human and murine osteoclast differentiation. Murinized human and murine scRNA‐seq data were integrated. (A) Feature plot (UMAP) demonstrating a total of 13 clusters during osteoclastogenesis (murinized human, n = 4285; mouse, n = 4388; total, n = 8673 cells). Cell‐based trajectory analysis in the integrated data set. Clusters and annotation in the murinized human and murine OLC data set. (B) Feature plots and cell‐based trajectory in the reconstructed human and murine OLC data set. Arrows indicate the branch into OLC. (C) Phase analysis of cell‐cycle status in the reconstructed human and mouse OLC data set (red, G1; green, G2/M; blue, S phase). (D) Heat map analysis of the top 25 differentially expressed genes in each OLC cluster. Osteoclastic markers Acp5, Ctsk, Atpv0d2, Nfatc1, Ocstamp, and Oscar are highly expressed in clusters 5 and 6. (E) Feature plot of marker genes. Cd14, Cxcr4, Cx3cr1, and Ccr2 are progenitor markers. Tnfrsf11A, Nfatc1, Ctsk, and Atp6v0d2 are osteoclast‐associated genes. (F) Dot plot of GO terms enriched in each cluster related to “dendritic cell,” “macrophage,” “osteoclast,” and “bone.” (G) Expression of detectable Rab family genes in the integrated data set. Note that Rab38 is predominantly expressed in osteoclastic clusters 5 and 6.

Fig. 6

ChIP‐seq profiling of human and murine osteoclast precursors and mature osteoclasts. Data represent the ChIP‐seq results for the active chromatin marker H3K4me3 (histone H3 lysine 4 trimethylation), the repressive chromatin marker H3K27me3 (histone H3 lysine 27 trimethylation), the active enhancer marker H3K27ac (histone H3 lysine 27 acetylation), and the osteoclast‐specific marker NfatC1. (A) Histone methylation and acetylation status and the binding of NFATc1 at the Ctsk and Rab38 genes, as analyzed by ChIP‐seq from day 0 to day 14 following RANKL stimulation in human osteoclast precursors and mature osteoclasts. The histone modification patterns changed from H3K4me3/H3K27me3 bivalent to H3K4me3 monovalent at both Ctsk and Rab38 following RANKL stimulation. The histone acetylation H3K27ac and the binding of NFATc1 increased at both the Ctsk and Rab38 gene regions. Bold bars indicate +/− 1 Kb from TSS. (B) Histone methylation and acetylation status and the binding of NFATc1 at the Ctsk and Rab38 genes were analyzed by ChIP‐seq at day 0 or day 2 and day 3 following RANKL stimulation in murine osteoclast precursors and mature osteoclasts. The level of H3K4me3 increased at the Ctsk gene. The histone modification patterns changed from H3K4me3/H3K27me3 bivalent to the H3K4me3 monovalent at the Rab38 gene following RANKL stimulation. The histone acetylation H3K27ac increased at the Rab38 gene, and the binding of NFATc1 increased at both the Ctsk and Rab38 gene regions. Bold bars indicate +/− 1Kb from TSS. (C,D) Gene expression was analyzed by RNA‐seq in human and murine osteoclast precursors and mature osteoclasts. Ctsk and Rab38 gene expression was analyzed by RNA‐seq at day 0 and day 14 following RANKL stimulation to induce human osteoclast differentiation. Gene expression of Ctsk and Rab38 was analyzed by RNA‐seq at day 0 and day 3 following RANKL stimulation during murine osteoclast differentiation. ChIP‐seq = chromatin immunoprecipitation‐sequencing; RANKL = receptor activator of nuclear factor‐kappa‐Β ligand; TSS = transcription start site.

Fig. 7

Effect of Rab38 knockdown by retroviral shRNA on murine osteoclast differentiation. (A) Expression of Rab38 genes in mouse osteoclast differentiation measured by quantitative real‐time PCR. (B) Expression of Rab38 protein in mature osteoclasts (green; Rab38, blue; nucleus; scale bar indicates 50 μm). (C,H) Expression of Rab38 gene under conditions of Rab38 knockdown by retroviral shRNA (C) or by lentiviral Cas9/gRNA (H) and the controls. (D,I) Counts of TRAP‐positive cells (≥3 nuclei). (E,F,J) TRAP staining (E,F,J) and cytoskeletal staining by rhodamine‐phalloidin (F,J) of osteoclasts under the respective conditions of Rab38 knockdown. (G,K) Actin ring formation rate. Ten cells in each condition were measured. (L) Schematic representation of the regulation of osteoclastogenesis by the Nfatc1–Rab38 axis. The values were measured in triplicate and compared against control (*p < 0.05, ***p < 0.001).