Differentiating human osteoclasts
Osteoclast differentiation is a coordinated process starting with myeloid commitment of hematopoietic stem cells to monocytes that differentiate into macrophages following macrophage colony-stimulating factor (M-CSF) signaling. M-CSF acts on its receptor, colony-stimulating factor 1 receptor (CSF1R), which induces expression of RANK22, a receptor for RANKL. The RANK-RANKL binding activates crucial pathways necessary for macrophages to differentiate into pre-osteoclasts including the nuclear factor κβ (NFκβ) and MAPK signaling pathways23,24. Activation of these pathways increases NFATc1 expression, which drives osteoclast differentiation and is thus a key regulator of osteoclastogenesis25. RANKL signaling then causes pre-osteoclasts to mature and fuse into multinucleated bone resorbing osteoclasts26, thereby fulfilling the differentiation and maturation of osteoclasts. Notably, osteoprotegerin (OPG), which is released from stromal cells and osteoblasts, is a decoy receptor for RANKL that prevents RANK-RANKL binding and thereby inhibits RANKL-mediated signaling in osteoclasts27.
Although these regulatory mechanisms of osteoclastogenesis are well-described, little is known about the overall changes in gene expression that occur during human osteoclastogenesis. Therefore, we first established an overview of the transcriptional reprogramming of human osteoclastogenesis using osteoclasts differentiated from human CD14+-monocytes. These were isolated from peripheral blood from eight anonymous female blood donors aged 18–49 years using concomitant treatment with M-CSF for 9 days (days 0–9) and RANKL for seven days (days 2–9). RNA was harvested at day 0 and days 2, 5, and 9 after induction of differentiation (Fig. 1A). Cell culture media was analyzed for the osteoclast marker tartrate-resistant acid phosphatase (TRAcP) activity at day 7 and 9. Furthermore, resorptive activity was determined by seeding and incubating osteoclasts at day 9 onto bovine bone slices for 72 hours.
We found that cells from all donors differentiated into TRAcP-expressing multinucleated osteoclasts capable of resorbing bone and as previously seen28,29. TRAcP and resorptive activity varied between donors (Fig. 1B). RNA-seq analyses revealed time-dependent downregulation of monocyte-specific genes (e.g., TREM1, SELL, CLEC10A)30 and monocyte-macrophage markers (e.g., ADGRE1)31 and a like-wise time-dependent increase in the expression of osteoclast-specific marker genes (e.g., CTSK, ACP5, MMP9, CA2) from day 0 to day 9 (Fig. 1C). Analyses of the variance in gene expression at different time points demonstrated 10839 genes that were differentially expressed, i.e., the level of gene expression changed between at least two time points across donors. For some of these genes, such as COX11, the donors showed a highly similar temporal expression pattern whereas expression of other genes, e.g., LPCAT7, showed substantial inter-donor variation (Fig. 1D). To avoid donor-dependent gene expression patterns, we only included differentially expressed genes with high similarity across donors (Fig. 1E & F). This resulted in 8446 genes that were differentially expressed during osteoclast differentiation and formed the basis of subsequent analyses.
RNA-seq reveals temporal gene expression patterns with distinct cellular functions and implications for human bone biology
Using k-means clustering, we observed that the 8446 genes could be grouped into eight temporal gene expression patterns (clusters) (Fig. 2A). Among these, two clusters were characterized by a transient decrease or increase in gene expression levels, and six clusters were characterized by early (peak at day 2), middle (peak at day 5), or late (peak at day 9) changes in gene expression levels. Gene ontology (GO) and reactome pathway analyses revealed cluster specific enrichments for distinct biological processes (Fig. 2B-C) such as metabolic reprogramming and mitochondrial activation, which were linked to early upregulated genes (cluster 2) (Fig. 2B). In line with this, we found distinct temporal profiles for metabolic genes involved in glucose metabolism and tricarboxylic acid cycle (TCA)-mediated metabolizing of fatty acids (Fig. 2C). Downregulated genes were on the other side barely linked to metabolic processes but important for cytokine production and immune cell activation. We also found that several genes in the late upregulated cluster 4 were involved in mature osteoclast related processes such as cytoskeleton organization, pH regulation, bone remodelling, cell migration, and cell-cell fusion32–34. Furthermore, cluster 4 contained genes that were upregulated during late stages of osteoclast differentiation known to be important for osteoclast function, such as CTSK, ACP5, DCSTAMP, and CA2 (ref. 28,35–38). Based on this, cluster 4 contained genes characterized as mature osteoclast specific genes.
To link the differentially expressed gene patterns with bone biology, we then investigated if there was an association between the osteoclast gene expression profiles and genes associated with bone mineral density (BMD), bone mineral content (BMC), and/or bone morphology in knockout mouse models. Based on data from the International Mouse Phenotyping Consortium (IMPC)39, knockout phenotypes characterized by increased BMC and BMD (based on Dual-energy X-ray absorptiometry (DXA)) were linked to genes repressed early during osteoclast differentiation, i.e. from monocyte to macrophage (clusters 5, 6, and 7), while knockout phenotypes characterized by decreased BMC and BMD were linked to genes repressed during late osteoclast differentiation (cluster 8) (Fig. 2D). Gene knockout phenotypes characterized by other skeletal defects on X-ray were linked to genes induced in mature osteoclasts (cluster 4) and genes with decreased expression during osteoclastogenesis (clusters 5–8). While disrupting osteoclast function is linked to bone phenotype causing-mutations in mice, these data indicate existence of specific associations between subgroups of bone phenotypes and temporal gene expression patterns during osteoclast differentiation.
Next, we tested whether the differential expressed gene patterns were aberrantly expressed in patients with osteoporosis. Using published microarray data on RNA in iliac crest bone tissue from 27 osteoporotic and 39 non-osteoporotic postmenopausal women40 we found that expression of ACTN2, a gene related to osteoclast fusion41 and a member of cluster 4, was higher in bone from osteoporotic women (Fig. 2E). In a genome-wide context, we found an overlap between genes that were upregulated in samples from osteoporotic women and genes that were induced (clusters 2 and 4) during osteoclast differentiation and vice versa for genes that were repressed in both data sets (clusters 7 and 8) (Fig. 2F). This implicated higher abundance of osteoclast-specific markers and absence of progenitor related genes in bone in osteoporotic women. To further elucidate this, we used published micro array-data from peripheral blood monocytes (PBMs) from 73 non-osteoporotic pre- and postmenopausal Caucasian women grouped into low or high hip BMD assessed using DXA42,43. However, in contrast to clear enrichment patterns observed in the expression data from postmenopausal osteoporotic iliac crest biopsies, we could not find an association between genes that were upregulated (cluster 1–4) or downregulated (cluster 5–8) during osteoclast differentiation and genes with altered expression in PBMs from non-osteoporotic pre- and postmenopausal women with low or high BMD (data not shown).
Finally, we investigated if the osteoclast gene expression profiles were linked to genetic sequence variations reported to be associated with human BMD estimated by heel quantitative ultrasound (eBMD). Using previously published data44, we found that the density of SNPs associated with eBMD was greater nearby genes with differential expression during osteoclast differentiation compared to genomic background (random distribution in the genome). This was especially the case for genes repressed after M-CSF-stimulation (cluster 7) and for late-induced genes representing mature osteoclast (cluster 4) (Fig. 2G). These data indicate that genes linked to monocyte progenitors and mature osteoclasts, i.e., that strongly change expression following differentiation from macrophage stage, are more likely to be affected by human sequence variations associated with alterations in eBMD.
Transcriptional networks drive osteoclast differentiation
We next investigated if the temporal changes in gene expressions were linked through regulatory networks, i.e., that gene induction during late stages of osteoclast differentiation was a consequence of early changes in gene expression. We applied the machine learning algorithm “Integrated System for Motif Activity Response Analysis” (ISMARA)45 to model transcription factor activity and to predict target genes of each transcription factor using gene expression data and motif occurrences in promotor regions. Importantly, the algorithm predicted a strong increase in the motif activity of two key transcription factors of osteoclast differentiation, NFATC1 and JUN46 (Fig. 3A). We also found that more than half of the 682 transcription factors that are part of the ISMARA motif database and our expression data set changed activity during differentiation (Fig. 3B), which included many of the transcription factors annotated in the GO terms “osteoclast differentiation” or “bone resorption”. Besides the changes in activity, ISMARA predicts target genes for a given transcription factor, i.e., genes that are very likely to depend on the presence of that specific factor in the given expression data set. In accordance with the continuous increase in JUN and NFATC1 activity during osteoclast differentiation and their documented role in the maturation of osteoclasts46, we found most of their predicted target genes to be upregulated during the late phase of differentiation, i.e., members of cluster 4 (Fig. 3C). To test this model for target gene prediction, we performed co-expression analysis of transcription factors and their predicted target genes on published single cell RNA-seq (scRNA-seq) data from in vitro differentiated human osteoclasts47. Due to the sparse nature of scRNA-seq data, we first performed deep clustering to estimate the correlation of co-expression of transcription factors and predicted targets at cluster level (Fig. 3D). As exemplified by the predicted target genes for MYB related protein 2 (MYBL2), a transcription factor known to regulate cell cycle genes48, which were enriched in cluster 3 (data not shown), we found a strong overlap between the expression of MYBL2 and the predicted target genes at both single cell (Fig. 3E) and cluster level (Fig. 3F). Importantly, similar associations were only observed among transcription factors with cluster specific expression patterns (Fig. 3G), supporting that expression of target genes follows the expression pattern of the regulating transcription factors.
We next constructed a directed transcriptional network by combining the predicted target genes of all 682 transcription factors part of the ISMARA database using network enrichment analysis49. This demonstrated that members of the early, transiently induced gene cluster (cluster 1) were particularly important for the regulation of genes of the early and middle-induced clusters (cluster 2 and 3), which in turn control the induction of osteoclast-specific genes (cluster 4) (Fig. 3H). Using the transcriptional network approach, we showed that the proportion of transcription factors highly relevant for the regulation of mature osteoclast genes (cluster 4) and / or causal eBMD genes44 were independent of the factor itself changing expression levels or being a causal eBMD gene (Fig. 3I). As an example, we identified transcription factors, including members of the HOXA and MEF2 family, such as HOXA10 and MEF2A, which were predicted to be important for osteoclast and bone function despite not being differentially expressed or associated with eBMD SNPs (Fig. 3I). Taken together, these network analyses demonstrated associations in the stepwise remodelling of the transcriptional networks regulating osteoclast differentiation, i.e., feedforward loops in gene regulation to overcome transitions between states of cellular differentiation in osteoclasts. In addition, our data suggested that network-based analysis can predict transcription factors important for osteoclast function based on the level of transcription factor activity, i.e., factors that contribute to the regulation of osteoclast genes through post-transcriptional mechanisms.
Subgrouping osteoclast transcriptional networks
As it was recently suggested that mature murine osteoclasts can undergo fission into osteomorphs, which are transcriptionally distinct from osteoclasts and macrophages and can recycle back into osteoclasts19, we next questioned whether our gene expression data could be used to subgroup osteoclast transcriptional networks and give insight into the transcriptional regulation of osteomorph-related genes. We focused on a previously published gene signature of 82 genes characterized by higher expression in murine osteomorphs than in osteoclasts and their progenitors19, and for which we could detect the human homologs in our expression data. We found that osteomorph-related genes were enriched among gene clusters with transient activation and repression or gene clusters that were upregulated at day 5 of osteoclast differentiation (cluster 1,3, and 5) (Fig. 4A). Based on our bulk RNAseq data, osteomorph-related genes were enriched among differentially expressed gene profiles with high expression at late-stage differentiation (cluster 3 and 5) but did not follow the same expression pattern as classical mature osteoclast-genes like CTSK, ACP5, and DCSTAMP (cluster 4) (Fig. 4A). Therefore, we speculated if osteomorph-related genes19 could be enriched among genes originally characterized as mature osteoclasts in previously published scRNA-seq data including data from human monocytes and mature osteoclasts47. First, our analyses showed that the distribution of genes expressed in mature osteoclasts were sub-grouped in two pools of cells (left and right arm) (Fig. 4B) of which both arms contained mature osteoclast-specific genes and originated from a group of cells enriched in monocyte specific genes (middle connective piece) (Fig. 4C). Combining published scRNA-seq data from osteomorphs19 with signature genes for the three groups (Fig. 4D), we found that osteomorph-related genes were enriched among the left arm containing mature osteoclast genes (Fig. 4E). However, since the enrichment was only due to the overlap of five genes (Fig. 4F) in the left arm and none in the right arm, the evidence for the left arm representing osteomorphs rather than mature osteoclasts was considered limited. These findings let us to further elucidate other potential differences between the left and right arms. Therefore, we assessed biological functional differences between the two arms and found the left arm-selective genes to be enriched among biological processes such as lipid metabolism, cytokine production and cell migration whereas the right arm-selective genes were enriched among mitochondrial processes, ATP-transport, and translation (Fig. 4G). Notably, left- and right-arm selective genes also showed distinct temporal expression patterns in our bulk RNA-seq profiles (Fig. 4H) and importantly, could be clustered apart using our gene regulatory network (Fig. 4I). Here, the well-known osteoclast transcription factors JUN and NFATc1 showed importance for mature osteoclast genes independent of the subpopulations while other osteoclast-associated transcription factors (GO terms ‘osteoclast differentiation’ and ‘bone resorption’) showed specificity for the genes of the left (e.g., CEBPb, FOS, BCL6) or right arm (e.g., ATF1, NRF1, SIX5). Although it is unknown if these differences were caused by the left arm also containing immature and non-fused mature osteoclasts and the right arm containing mainly fully fused osteoclasts, these analyses demonstrated that subpopulation-selective osteoclast genes have distinct temporal expression patterns, and that osteoclast-subtype specificity is defined by the interaction of a core osteoclast transcriptional network (e.g., NFATc1) with cell-type specific ones.
Differentially expressed GPCRs during human osteoclast differentiation
In the perspective of identifying potential anti-osteoporotic treatment targets, we then speculated if we could use our data set to detect and test previously unrecognized regulators of osteoclast differentiation and activity. To test this, we focused on the role of GPCRs during the transcriptional reprogramming of osteoclast differentiation as GPCRs implicate with bone biology and represent feasible therapeutic targets accounting for approximately 34% of FDA-approved drugs (2017) that target a total of 108 different GPCRs50. An example of this is teriparatide, a bone anabolic anti-osteoporotic drug that targets the parathyroid hormone (PTH) receptor on osteoblasts or their precursors51. Furthermore, 36 GPCRs are known to be associated with altered BMD, bone morphology, and bone-related diseases such as arthritis in humans52. For example, activation of the calcitonin receptor, encoded by CALCR, reduces osteoclast activity53, and polymorphisms in the CALCR gene locus are associated with changes in lumbar spine BMD and fracture risk in postmenopausal Korean women54.
Based on the GPCR database (GPCRdb)55, we identified 135 GPCRs and 42 gene-encoded ligands expressed during osteoclast differentiation (Fig. 5A). Importantly, GPCRs and ligands as well as ligand-receptor pairs were enriched among the differentially expressed genes (Fig. 5B & 5C). Interestingly, when focusing on all human encoded non-olfactory GPCRs, we did not find them to be enriched among the 835 genes predicted to be causally related to eBMD44. However, when subgrouping those GPCRs that were expressed or differentially expressed during osteoclast differentiation, we observed that these were enriched among genes predicted to be casually related to eBMD (Fig. 5D). Among the differentially expressed GPCRs, we found well-described regulators of osteoclast function like CALCR; LGR4, which encodes leucine-rich repeat containing G-protein coupled receptor 4 and acts as a receptor for RANKL thereby inhibiting osteoclast differentiation and activity56; and GIPR, which encodes the receptor of the incretin hormone glucose-dependent insulinotropic polypeptide (GIPR) that inhibits osteoclast differentiation and activity29. Taken together, these data indicate a pivotal role for differentially expressed GPCRs and their ligands in the regulation of osteoclast differentiation and activity, and BMD.
Based on this, we hypothesized that other differentially expressed GPCRs could affect osteoclast differentiation or function. To test this hypothesis, we focused on differentially expressed GPCRs with known downstream G protein-coupling57, which left us with 64 GPCRs (Fig. 5E). Of these, we selected three (marked in red in Fig. 5E), namely C5AR1 (complement 5a receptor 1), SSTR2 (somatostatin receptor 2), and FFAR4 (free fatty acid receptor 4/GPR120), as all receptors had commercially available agonists and the regulatory effects on human osteoclast differentiation and activity had not previously been investigated.
C5AR1 activation impairs osteoclast differentiation due to its progenitor-restricted expression
The expression pattern of C5AR1, a cluster 7 member, was characterized by a strong decrease during osteoclast differentiation (Fig. 6A). C5AR1 expression could be detected in scRNA-seq from human monocytes and osteoclasts47 in which expression was limited to undifferentiated cells (Fig. 6B), indicating that C5AR1 is not a functional receptor on mature osteoclasts. C5AR1 is a class A GPCR, which predominantly couples to inhibitory G proteins (Gi) following ligand binding. This leads to a decrease in intracellular cAMP levels58, and e.g., mediates innate immunity and promotes inflammation. Mice with C5ar1 knockout have increased bone mass, reduced osteoclast formation and increased osteoblast number compared to wild-type mice59. In cellular models, adding supernatant from C5AR1-activated murine osteoblasts enhances murine osteoclastogenesis60, and activating C5AR1 with C5a on human osteoblast-like cells induces interleukin-6 production, which may trigger bone resorption61. Thus, C5a may couple osteoblast and osteoclast activity, and C5AR1 could have osteoclastogenic properties in human cells.
To determine the effects of C5AR1 on human osteoclast differentiation and activity, we treated differentiating osteoclasts and mature osteoclasts respectively with the selective agonist BM221 (ref. 62). We found that activation of C5AR1 with BM221 enhanced osteoclastogenesis by increasing the number of newly formed osteoclasts (Fig. 6C). Importantly, this effect was reversed by adding PMX205, a C5AR1 antagonist63 (Fig. 6C). BM221 did not affect the number of nuclei per osteoclasts or the resorptive activity of mature osteoclasts (Fig. 6D), demonstrating that C5AR1 activation solely affect the quantity of osteoclasts rather than the resorptive potential of the cell.
SSTR2 activation did not change osteoclast differentiation, but decreased osteoclast activity
The expression pattern of SSTR2, a member of cluster 5, was transiently decreased during osteoclast differentiation with high expression at day 0 in some, but not all donors, and then upregulated during late differentiation among all donors (Fig. 6E). Unlike C5AR1, expression of SSTR2 was not detected in scRNA-seq data of human osteoclastogenesis. SSTR2 is a class A GPCR, which couples to inhibitory G proteins (Gi or Go) following ligand binding64 leading to a decrease in intracellular cAMP levels. SSTR2 activity can be modulated by the somatostatin analogue octreotide that is used clinically to control growth of and hormone secretion from neuroendocrine tumors such as insulinomas and growth hormone secreting pituitary tumors65. Furthermore, octreotide has been shown to acutely increase both markers for bone resorption and formation in healthy subjects66.
To determine the function of SSTR2 in human osteoclasts, we first measured cAMP signaling after exposure to either somatostatin or octreotide with or without pre-treatment with the specific SSTR2 antagonist BIM-23627 (ref. 67) using cAMP-Glo assays in mature (day 9) osteoclasts (Fig. 6F). We observed that somatostatin and octreotide dose-dependently decreased cAMP levels, although somatostatin only had effects at the highest concentration investigated (100 nM). Importantly, the effects on cAMP were reversed when cells were pre-treated with BIM-23627. Treatment with somatostatin or octreotide did not affect osteoclast differentiation assessed as total number of mature osteoclasts or the number of nuclei per cell in mature osteoclasts (Fig. 6G). However, both somatostatin (100 nM) and octreotide (10 nM) inhibited the bone resorptive activity of mature osteoclasts, and the effect was abolished when osteoclasts were pre-treated with BIM-23627 (100 nM) (Fig. 6H). Taken together, these findings demonstrate that SSTR2 is a functional differentially expressed GPCR, that has anti-resorptive effects on mature human osteoclasts.
FFAR4 activation impaired osteoclast differentiation and decreased resorptive activity of mature osteoclasts
The expression pattern of FFAR4, a member of cluster 2, peaked at day 2 and remained high throughout differentiation (Fig. 6I). Like SSTR2, expression of FFAR4 was not detected in scRNA-seq data of human osteoclastogenesis. FFAR4 is a class A GPCR, which upon ligand binding stimulates the generation of inositol 1,4,5-trisphosphate (IP3). This induces intracellular calcium mobilization and phosphorylation of ERK via Gq/11-mediated signaling pathways68. Ffar4 expression in mice is strongly upregulated during osteoclastogenesis and negatively regulates osteoclast differentiation, function, and survival69. Treating human CD14+-monocytes with docosahexaenoic acid (DHA), an n-3 fatty acid ligand for several free fatty acid receptors including FFAR4, is reported to inhibit osteoclast formation and activity, and downregulate expression of genes important for osteoclast function like DCSTAMP and NFATC1 (ref. 70). However, as it remained unknown if these inhibitory effects of DHA were mediated via FFAR4 activation, we extended the analyses of the role of FFAR4 agonism on osteoclast differentiation and activity using the selective FFAR4 agonist TUG891 (ref. 71).
To determine the functionality of FFAR4 in human osteoclasts, we first measured IP1 levels, a stable downstream metabolite of IP3, upon exposure to TUG-891 using an HTRF-assay. We observed a significant increase in IP1 for all concentrations investigated (0.1, 1, and 10 nM), with the most potent effect at 10 nM (Fig. 6J). Next, using 10 nM TUG891 we showed that FFAR4 activation impaired osteoclast differentiation by reducing the number of osteoclasts, but not the number of nuclei per osteoclast (Fig. 6K). Finally, treating mature osteoclasts seeded on bone slices with TUG-891 inhibited the resorptive activity (Fig. 6L), demonstrating that FFAR4 activation inhibits osteoclastogenesis by limiting initial maturation of osteoclasts and the resorptive activity of mature human osteoclasts.