Dcx expression defines a subpopulation of Gdf5 + cells with chondrogenic potentials in E12.5 mouse embryonic limbs.

Growth differentiation factor 5 (Gdf5) and doublecortin (Dcx) genes are both expressed in joint interzone cells during synovial joint development. In this study, we re-analyzed the single cell RNA-sequencing data (Gene Expression Omnibus GSE151985) generated from Gdf5 + cells of mouse knee joints at embryonic stages of E12.5, E13.5, E14.5, and E15.5, with a new focus on Dcx. We found that Dcx expression was enriched in clusters of Gdf5 + cells, with high expression levels of pro-chondrogenic genes including sex determining region Y-box transcription factor 5 (Sox5), Sox6, Sox9, Gdf5, versican, matrilin 4, collagen type II α 1 chain (Col2a1), Col9a1, Col9a2, and Col9a3 at E12.5. Dcx + and Dcx - cells had differential gene expression profiles. The up-regulated genes in Dcx + vs. Dcx - cells at E12.5 and E13.5 were enriched in chondrocyte differentiation and cartilage development, whereas those genes up-regulated at E14.5 and E15.5 were enriched in RNA splicing, protein stability, cell proliferation, and cell growth. Gene expression profiles in Dcx + cells showed rapid daily changes from E12.5 to E15.5, with limited number of genes shared across the time period. Expression of Gdf5, Sox5, Sox6, melanoma inhibitory activity, noggin, odd-skipped related transcription factor 2, matrilin 4, and versican was positively correlated with Dcx expression. Our results demonstrate that Dcx expression defines a subpopulation of Gdf5 + cells with chondrogenic potentials in E12.5 mouse embryonic limbs.

analysis of variance

collagen

doublecortin

growth differentiation factor 5

Gene Set Enrichment Analysis

kinase insert domain receptor

matrilin

melanoma inhibitory activity

noggin

odd-skipped related transcription factor 2

platelet derived growth factor receptor α

super-cluster

single cell RNA-sequencing

Sex Determining Region Y-Box transcription factor

-Distributed Stochastic Neighbor Embedding

versican

wingless-type MMTV integration site family member

Introduction

It has been well established that embryonic articular joint development starts from appearance of joint interzone – a densely packed region of flattened cells in the cartilaginous anlage [1]. The interzone contains a central intermediate layer and two chondrogenous layers that eventually differentiate into articular chondrocytes. In mouse embryos, forelimb buds appear at embryonic stage of 9.5 days postcoitus (i.e., E9.5) and mesenchymal condensation occurs at E12.5 [2]. Interzones appear at the presumptive sites of articular joints at E13.5 [3]. At E15.5, joint cavitation appears, marking formation of a synovial joint [3]. Growth differentiation factor 5 (Gdf5) is first expressed in the forelimb at E11.5 at the presumptive regions of the shoulder and elbow and its expression decreases at E15.5 [4]. Gdf5-Cre reporter mice demonstrate that Gdf5+ interzone cells differentiate into articular chondrocytes, synovial lining, and intrajoint ligaments [[5], [6], [7], [8]]. Bian et al. recently isolated single Gdf5+ cells from mouse knee joints at E12.5, E13.5, E14.5, and E15.5, and performed single cell RNA-sequencing (scRNA-seq) analysis. They found that Gdf5-lineage enriched (GLE) cells contain three super-clusters (SC): SC1 comprises a mixture of progenitor cells of mesenchymal character and chondroprogenitors; SC2 contains a mixture of interzone cells and transient chondrocytes; and SC3 largely consists of fibroblast-related cells. Their work further demonstrates the heterogeneity of Gdf5+ cells at the transcriptional level [9]. How the clusters of Gdf5+ cells commit to different lineages remains to be determined.

Doublecortin (Dcx) is originally identified in migrating and differentiating neurons [10]. In mouse limbs, Dcx is initially expressed in the common mesenchymal precursors at E9.5 through E12.5. Dcx expression is sequestered in the joint interzones at E13.5 and then in immature articular chondrocytes, but Dcx is not expressed in the synovium, ligaments/tendons, and menisci [11,12]. Thus, Dcx appears to be a unique marker of articular chondrocyte lineage [7,[13], [14], [15], [16], [17], [18], [19]]. It is hypothesized that Dcx+ cells are a subpopulation of Gdf5+ cells in the joint interzones. The present study is to test this hypothesis through re-analyzing the scRNA-seq data publicly deposited by Bian et al. [9].

Materials and methods

Data source and identification of the target cells. Raw data were generated using scRNA-seq analysis of Gdf5+ cells isolated from mouse (FVB/NJ crossing with C57BL/6J background) embryonic knee joints at E12.5, E13.5, E14.5, and E15.5, and were deposited to Gene Expression Omnibus under accession number GSE151985 by Dr. Patrick Cahan’s lab at Johns Hopkins School of Medicine; Gdf5+ cells expressed enhanced yellow fluorescent protein and were isolated using fluorescence-activated cell sorting (FACS) (see details as described previously [9]). R (version 4.0.5, R Foundation for Statistical Computing, Vienna, Austria)-based Seurat (version 4.0.2) [20] was used to analyze the raw data. We excluded potential cell doublets through removing the cells of the top 5% of total count and excluded potential low-quality libraries that had fewer than 500 genes or >5% of mitochondrially encoded genes. SingleR (version 1.4.1) [21] was used to annotate the scRNA-seq data into 18 main cell types, such as fibroblasts, endothelial cells, macrophages, etc, based on mouse cell type reference database [22]. Fibroblasts were included while the other cell types were excluded from further analysis. Seurat’s t-Distributed Stochastic Neighbor Embedding (tSNE) function was used to cluster the fibroblasts. Like in the previous study [9], we excluded the clusters with >20% cells expressing myogenic differentiation 1 (Myod1, known to be expressed mainly in myoblasts), premelanosome protein (Pmel, known to be expressed mainly in melanocytes), sex determining region Y-box transcription factor 10 (Sox10, known to be expressed mainly in neural crest cells), and/or twist basic helix-loop-helix transcription factor 2 (Twist2) and iroquois homeobox 1 (Irx1) (both genes were known to be expressed mainly in dermis cells). We eventually included a total of 8514 cells in our analysis (Supplementary Table S1), which was 1185 cells more than what were used for downstream analyses by Bian et al. [9]. We speculated that the difference might be due to the difference in algorithms of cell annotation used by Bian et al. [9] and SingleR [21] used here in our analysis.

Analysis of the scRNA-seq data. Software packages used to generate each figure are summarized in Supplementary Fig. S1. Briefly, Seurat was used to generate tSNE plots and visualize cell clusters, and its Dotplot function was used to visualize average expression levels of Dcx and percentages of the cells with Dcx expression. FindMarkers function of MAST package (version 1.16.0) [23] was used to identify 2000 highly variable genes, and their average expression was ranked from the highest (top) to the lowest (bottom) levels in the clusters with the highest Dcx expression. Seurat’s DoHeatmap function was used to generate heatmaps in Fig. S2 and FeaturePlot function was used to separate Dcx+ and Dcx- cells. FindMarkers function was used to identify the genes differentially expressed between Dcx+ and Dcx- cells. Average log2 (Fold Change, threshold −2.09 to 1.15) was calculated and Student’s t-test was used to calculate p values (<0.05 was considered significant). The ggplot function of ggpubr (version 0.4.0 created by Alboukadel Kassambara) was used to generate volcano plots based on average log2(Fold Change) and –log10(p value). A set of significantly upregulated genes were analyzed in Gene Set Enrichment Analysis (GSEA) using enrichR package (version 3.0 created by Wajid Jawaid). In order to compare the gene expression levels across the four embryonic stages, we integrated the datasets of four embryonic stages using Seurat’s integrateData function. Quality control procedures were performed as described above to exclude the low quality libraries and unwanted cell types, resulting in 9319 target cells. These cells were separated into Dcx+ and Dcx- using Seurat’s FeaturePlot function. We selected 139 genes that had been linked to joint development or shown in the previous analysis [9]. The genes were ranked based on their average expression levels in Dcx+ cells at E12.5. The superheat package (version 1.0.0 created by Rebecca Barter and Bin Yu) was used to generate heatmaps in Fig. 3, Fig. 4A, Fig. S4, and Fig. S5, because Seurat’s DoHeatmap function does not allow manual selection of individual genes. Eight genes, i.e., Gdf5, Sox5, Sox6, melanoma inhibitory activity (Mia), noggin (Nog), odd-skipped related transcription factor 2 (Osr2), matrilin 4 (Matn4), and versican (Vcan), showed differential expression patterns similar to Dcx. Matn4 is mainly expressed in articular chondrocytes [24], in contrast to Matn1 that is mainly expressed in endochondral chondrocytes but not in articular chondrocytes [25]. Kendall’s cor.test function in ggpubr package was used to perform correlation between each of the eight genes and Dcx.

Fig. 3

Gene expression profiles in Dcx+ cells present rapid changes from E12.5 to E15.5. (A, B, C, and D) Heatmaps of top 25 and bottom 25 genes expressed in Dcx+ cells. Genes shared across the four embryonic stages are highlighted in red font. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 4

Expression of a set of genes is positively correlated with Dcx expression. (A) Heatmap of 32 genes from the integrated datasets. (B to I) Kendall’s correlation analyses. τ, correlation coefficient.

Results and discussion

Dcx expression is enriched in a subpopulation of Gdf5+ cells

In order to demonstrate the heterogeneity of Gdf5+ cells, Seurat was used to analyze 2435 Gdf5+ cells of the knee joints at E12.5. A tSNE plot showed that the cells were assigned into 7 clusters (Fig. 1A). Cluster #2 had the highest percentage of cells (12.4%) that expressed Dcx (Fig. 1B), the average expression of which was significantly higher than the other six clusters (p = 0.001 based on analysis of variance -ANOVA). This indicates that Dcx expression is enriched in Cluster #2 of Gdf5+ cells at E12.5. At E13.5, tSNE plot showed that 4072 Gdf5+ cells were assigned into 10 clusters (Fig. 1C). Cluster #6 had the highest percentage of cells (5.1%) that expressed Dcx (Fig. 1D), the average expression of which was significantly higher than the other nine clusters (p = 1.99e-07 based on ANOVA). This indicates that Dcx expression is enriched in Cluster #6 of Gdf5+ cells at E13.5. At E14.5, tSNE plot showed that 965 Gdf5+ cells were assigned into 6 clusters (Fig. 1E). Cluster #0 had the highest percentage of cells (4.2%) that expressed Dcx (Fig. 1F), but Cluster #2 (4.0%), #3 (3.6%), and #4 (3.8%) also had relatively higher percentages of Dcx+ cells than Cluster #1 (1.6%) or #5 (1.1%). However, the average expression of Dcx was not statistically significant among the six clusters (p = 0.366 based on ANOVA). This indicates that Dcx expression is enriched in multiple clusters of Gdf5+ cells at E14.5. At E15.5, tSNE plot showed that 1042 Gdf5+ cells were assigned into 8 clusters (Fig. 1G). Cluster #2 had the highest percentage of cells (30.4%) that expressed Dcx (Fig. 1H), the average expression of which was significantly higher than the other seven clusters (p = 2e-16 based on ANOVA). This indicates that Dcx expression is enriched in Cluster #2 of Gdf5+ cells at E15.5.

Fig. 1

Dcx expression is enriched in a subpopulation of Gdf5+ cells. (A, C, E, and G) Cells were assigned to clusters using t-Distributed Stochastic Neighbor Embedding (tSNE) plots. (B, D, F, and H) Dotplots represent average expression levels of Dcx and percentages of cells with Dcx expression.

To demonstrate how the clusters with the highest Dcx expression were defined at transcriptional levels, we ranked levels of gene expression in the clusters and generated heatmaps to show the relative levels of those genes ranked at top 20 or bottom 20. The top ranked 20 genes in Cluster #2 at E12.5 included Sox5, Sox6, Sox9, Gdf5, Vcan, Matn4, collagen type II α 1 chain (Col2a1), Col9a1, Col9a2, and Col9a3 that are known to be involved in chondrocyte differentiation (Supplementary Fig. S2A). In contrast, the bottom ranked 20 genes in Cluster #2 at E12.5 included Col1a2 and Col3a1 that are related to fibroblasts (Supplementary Fig. S2A). Although the top ranked 20 genes were highly expressed in Cluster #2, they were expressed at low levels in the other six clusters. In contrast, the bottom ranked 20 genes in Cluster #2 were expressed at relatively high levels in the rest clusters (Supplementary Fig. S2A). Since the expression of top 20 genes showed more dramatic difference compared to the rest clusters, we expanded the heatmap to include the top 50 genes in Cluster #2 and found a similar contrast between Cluster #2 and the rest clusters (Supplementary Fig. S2B). These findings suggest that Cluster #2 of Gdf5+ cells at E12.5 had a distinct transcriptional profile that is pro-chondrogenic. We did similar analyses for E13.5 (Supplementary Figs. S2C–D), E14.5 (Supplementary Figs. S2E–F), and E15.5 (Supplementary Figs. S2G–H), and found that the clusters with the highest Dcx expression showed distinct transcriptional profiles compared to the other clusters.

Dcx+ cells have distinct gene expression profiles compared to Dcx- cells

Since the clusters contained both Dcx+ and Dcx- cells, we separated Dcx+ and Dcx- cells and compared their transcriptional profiles. A volcano plot showed that many genes including Sox9, Col2a1, Col9a2, Col9a3, and Vcan were significantly upregulated in Dcx+ cells compared to Dcx- cells at E12.5 (Fig. 2A, p < 0.05). In contrast, some genes, including Col1a1, phospholipid phosphatase 3 (Plpp3), dickkopf WNT signaling pathway inhibitor 2 (Dkk2), Lumican (Lum), and cysteine dioxygenase type 1 (Cdo1), were significantly downregulated in Dcx+ cells compared to Dcx- cells (Fig. 2A, p < 0.05). Similarly, we found significantly upregulated and downregulated genes in Dcx+ cells compared to Dcx- cells at E13.5 (Fig. 2B), E14.5 (Fig. 2C), and E15.5 (Fig. 2D). Of note, most of the significantly up- or down-regulated genes had about 1 fold change, which might be due to the fact that the cells were all Gdf5+ cells, thus their differences in gene expression were subtle, rather than dramatically different. These findings suggest that Dcx+ cells had distinct gene expression profiles compared to Dcx- cells. To explore the functions of the upregulated genes, we performed GSEA of the significantly upregulated genes at each embryonic stage. We found that at E12.5, the upregulated genes were enriched in chondrocyte differentiation, skeletal system development, embryonic morphogenesis, and cartilage development (Fig. 2E). These functions are highly related to chondrogenesis, suggesting that the upregulated genes in Dcx+ cells are involved the development of synovial joints. GSEA for E13.5 showed that the upregulated genes were enriched in mesenchymal cell proliferation, chondrocyte differentiation, and cartilage development (Supplementary Fig. S3A). GSEA for E14.5 showed that the upregulated genes were enriched in RNA splicing, mRNA processing, and protein stability (Supplementary Fig. S3B). GSEA for E15.5 showed that the upregulated genes were enriched in cell growth, epithelial cell proliferation, cell cycle arrest, and skeletal system development (Supplementary Fig. S3C). We noted that the upregulated genes at E12.5 and E13.5 were enriched mainly in chondrocyte differentiation and cartilage development, whereas the upregulated genes at E14.5 and E15.5 were enriched mainly in RNA splicing, protein stability, cell growth and proliferation. It is well known that E12.5 to E13.5 is the time when joint interzones emerge from the condensed mesenchyme. The joint interzones are destined to give rise to the joint, whereas the non-interzone mesenchymal cells become transient chondrocytes that eventually will be replaced by the bone. More importantly, the interzones maintain Dcx expression, whereas the non-interzone mesenchymal cells turn off Dcx expression at E13.5 [11]. If we assess our findings under the aforementioned context, we may reasonably speculate that the upregulated genes in Dcx+ cells at E12.5 and E13.5 may play critical roles in determining the choice of cell fates, i.e., to become articular chondrocyte lineage. In contrast, at E14.5 and E15.5, the cell fates have already been determined, thus the upregulated genes in Dcx+ cells serve to maintain cell proliferation and cartilage growth.

Fig. 2

Dcx+ cells have distinct gene expression profiles compared to Dcx- cells. (A, B, C, and D) Volcano plots represent significantly up-regulated and down-regulated genes in Dcx+ cells compared to Dcx- cells. Log2 (Fold Change) > 0.263035 (i.e., fold change >1.2) was considered up-regulated and log2 (Fold Change) < −0.263035 (i.e., fold change <0.83) was considered down-regulated. It was considered significant if –log10 (p value) was >1.30103 (i.e., p value < 0.05). (E) Gene Set Enrichment Analysis of significantly up-regulated genes in Dcx+ cells at E12.5. The most significant functions based on gene ontology are shown.

Gene expression profiles in Dcx+ cells present rapid changes from E12.5 to E15.5

To further investigate the gene expression profiles in Dcx+ cells in comparison to Dcx- cells, we ranked the gene expression levels from the highest to the lowest in Dcx+ cells at individual embryonic stages and generated heatmaps to compare Dcx+ and Dcx- cells. We found that at E12.5, the top 25 genes in Dcx+ cells included Col9a2, Col2a1, Sox9, Vcan, Col9a3, Sox5, Col9a1, and Gdf5 (Fig. 3A), which were similarly highly expressed in Cluster #2 at E12.5 (Supplementary Fig. S2A). The bottom 25 genes in Dcx+ cells included Col1a2 and Col3a1 (Fig. 3A), which were similarly weakly expressed in Cluster #2 at E12.5 (Supplementary Fig. S2A). The results suggest that the gene expression profile in Cluster #2 at E12.5 is likely dominated by Dcx+ cells in the cluster. The top 25 genes in Dcx+ cells were weakly expressed in Dcx- cells, whereas the bottom 25 genes in Dcx+ cells were expressed at higher levels in Dcx- cells (Fig. 3A). Of note, the levels of gene expression in Dcx- cells were not dramatically different (mostly in yellow to green color codes around “0” – the average level). This pattern of differential expression was also observed when comparing the top 50 genes and bottom 50 genes in Dcx+ cells against Dcx- cells (Supplementary Figs. S4A–B). At E13.5, the top 25 and bottom 25 genes in Dcx+ cells were largely different from E12.5 (Fig. 3B), except for a few genes, including CD24 molecule (Cd24a), leucine rich repeats and immunoglobulin like domains 3 (Lrig3), high mobility group nucleosomal binding domain 2 (Hmgn2), integral membrane protein 2A (Itm2a), snail family transcriptional repressor 1 (Snai1), Vcan, zinc finger homeobox 3 (Zfhx3), Col1a2, Col3a1, and Col1a1 (Fig. 3B). The top 50 genes in Dcx+ cells were weakly expressed in Dcx- cells, while the bottom 50 genes in Dcx+ cells were expressed at higher levels in Dcx- cells (Supplementary Figs. S4C–D). At E14.5, the top 25 and bottom 25 genes in Dcx+ cells were completely different from E13.5 (Fig. 3C). Only one of the top 25 genes, calpain 6 (Capn6), was shared with the top 25 genes at E12.5. The top 50 genes in Dcx+ cells were weakly expressed in Dcx- cells, while the bottom 50 genes in Dcx+ cells were expressed at higher levels in Dcx- cells (Supplementary Figs. S4E–F). At E15.5, the top 25 and bottom 25 genes were mostly different from E14.5 (Fig. 3D). A few genes including mab-21 like 2 (Mab21l2), growth arrest specific 2 (Gas2), cytokine receptor like factor 1 (Crlf1), ankyrin repeat domain 37 (Ankrd37), and fibronectin 1 (Fn1) were shared with those of E14.5. Hmgn2, Vcan, Zfhx3, brain abundant membrane attached signal protein 1 (Basp1), and decorin (Dcn) were shared with those of E12.5 and/or E13.5 (Fig. 3A–D). The top 50 genes in Dcx+ cells were weakly expressed in Dcx- cells, while the bottom 50 genes in Dcx+ cells were expressed at higher levels in Dcx- cells (Supplementary Figs. S4G–H). The lists of top and bottom 50 genes with annotated functions are shown in Supplementary Table S2. At E13.5, some genes highly expressed in Dcx+ cells (Fig. 3B), e.g., Itm2a, Top2a, Lpar4, Nexn, and Nrk, were also highly expressed in cluster #6 with dominant Dcx expression (Fig. S2D). Similarly, at E15.5, Foxp1, Lgfbp7, Col26a1, Dpt, Gas6, Angptl1, Cxcl12, Apoe, Gas2, Dbi, Ptn, Mgp, and Vcan genes in Fig. 3D were shown in Cluster #2 with dominant Dcx expression (Fig. S2H). However, at E14.5, such concurrence was not found (comparing Fig.3C and Fig. S2F), which might be due to the fact that Dcx expression at E14.5 was distributed in 4 clusters (Fig. 1F). These findings demonstrate that the gene expression profiles in Dcx+ cells undergo rapid changes during the four-day period from E12.5 to E15.5, implying a quick differentiation of Dcx+ cells.

Expression of a set of genes is positively correlated with Dcx expression

To track the expression of genes from E12.5 to E15.5, we integrated the datasets and normalized the gene expression levels to the same scale across the four embryonic stages and generated a heatmap of 32 genes that have been studied in joint development [19]. We found that Sox4, Sox9, Sox11, Jun, paired like homeodomain 1 (Pitx1), Col2a1, Vcan, Matn4, and platelet derived growth factor receptor α (Pdgfra) were highly expressed from E12.5 to E15.5 in both Dcx+ and Dcx- cells, whereas scleraxis basic helix-loop-helix transcription factor (Scx), transforming growth factor beta receptor 2 (Tgfbr2), Indian hedgehog (Ihh), leucine rich repeat containing G protein-coupled receptor 5 (Lgr5), wingless-type MMTV integration site family member 4 (Wnt4), Wnt9a, Wnt16, erythroblast transformation-specific transcription factor Erg (Erg), Col22a1, and kinase insert domain receptor (Kdr) were weakly expressed in both Dcx+ and Dcx- cells (Fig. 4A). Wnt4 and Wnt9a expression was almost absent at E12.5 and slightly increased at E15.5, suggesting that Wnt signaling is not involved in interzone specification, but may be involved in joint maintenance. This finding is consistent with a previous study using Wnt4−/−Wnt9a-/- double mutant mice, which revealed that the mice developed joint fusion after joint interzone formation [26]. Sox9, Jun, and Col2a1 were expressed at higher levels in Dcx+ cells than Dcx- cells at E12.5, but this pattern was not kept at E13.5 to E15.5. On the other hand, Tgfbr2 and Matn2 expression was increased in Dcx+ cells only at E15.5 (Fig. 4A). Matn1 levels were slightly higher in Dcx+ cells than Dcx- cells at E13.5 to E14.5, but it was weakly expressed at E12.5 and E15.5. A set of eight genes including Gdf5, Sox5, Sox6, Mia, Nog, Osr2, Matn4, and Vcan were consistently expressed at higher levels in Dcx+ cells than Dcx- cells from E12.5 to E15.5 (Fig. 4A). Gdf5 expression varied from E12.5 to E15.5. Kendall’s correlation analysis showed that expression of these eight genes was positively correlated with Dcx expression (Fig. 4B–I). Heatmaps of other 48 genes (Supplementary Fig. S5A) and 59 genes (Supplementary Fig. S5B) showed some genes that were highly or weakly expressed in both Dcx+ and Dcx- cells. However, there was not any consistent pattern of differential expression between Dcx+ and Dcx- cells across the four embryonic stages from E12.5 to E15.5. These 107 genes were analyzed in the previous study by Bian et al. [9].

In summary, we found that Dcx expression was enriched in clusters of Gdf5+ cells expressing pro-chondrogenic genes. The up-regulated genes in Dcx+ cells were enriched in chondrocyte differentiation and cartilage development at E12.5 and E13.5, whereas those genes up-regulated at E14.5 and E15.5 were enriched in RNA splicing, protein stability, cell proliferation, and cell growth. Gene expression profiles in Dcx+ cells showed rapid daily changes from E12.5 to E15.5. Expression of Gdf5, Sox5, Sox6, Mia, Nog, Osr2, Matn4, and Vcan was positively correlated with Dcx expression from E12.5 to E15.5. These results demonstrate that Dcx expression defines a subpopulation of Gdf5+ cells with chondrogenic potentials in E12.5 mouse embryonic limbs, which provides clues to further investigate Dcx’s roles in joint development.

Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.