Distinct gene regulatory dynamics drive skeletogenic cell fate convergence during vertebrate embryogenesis

Abstract

Cell type repertoires have expanded extensively in metazoan animals, with some clade-specific cells being crucial to evolutionary success. A prime example are the skeletogenic cells of vertebrates. Depending on anatomical location, these cells originate from three different precursor lineages, yet they converge developmentally towards similar cellular phenotypes. Furthermore, their 'skeletogenic competency' arose at distinct evolutionary timepoints, thus questioning to what extent different skeletal body parts rely on truly homologous cell types. Here, we investigate how lineage-specific molecular properties are integrated at the gene regulatory level, to allow for skeletogenic cell fate convergence. Using single-cell functional genomics, we find that distinct transcription factor profiles are inherited from the three precursor states and incorporated at lineage-specific enhancer elements. This lineage-specific regulatory logic suggests that these regionalized skeletogenic cells are distinct cell types, rendering them amenable to individualized selection, to define adaptive morphologies and biomaterial properties in different parts of the vertebrate skeleton.

Fig. 1. A convergent transcriptomic signature in skeletogenic cells of different embryonic origins.

a Divergent diversifying versus convergent cell fate decisions. Precursor cells (P) usually differentiate into functionally distinct cell types (C_1-3). During vertebrate skeletogenesis, however, three distinct embryonic lineages—the neural crest (NC, magenta), the somitic mesoderm (SOM, blue), and the lateral plate mesoderm (LPM, yellow)—give rise to functionally similar skeletogenic cells (S). b Early embryonic origins and eventual anatomical locations of skeletogenic cells sampled in this study. c Sampling scheme to assess the transcriptional dynamics of skeletogenic convergence across anatomical locations (color-coded) and embryonic stages (Hamburger–Hamilton, HH). d–f tSNE representations of the three scRNA-seq datasets with ‘'broad'’ cell type annotations, with the mesenchymal populations used for ‘'fine'’ re-clustering highlighted in color. g Unsupervised hierarchical clustering and heatmap representation of pairwise Spearman’s rank correlation coefficients of highly variable genes in pseudobulk transcriptomes of '‘broad'’ cell type clusters. Anatomical origins of pseudobulks are indicated by color code. Mesenchymal cell populations across embryonic origins cluster together (red dotted square). h–k '‘Fine'’ re-clustering of mesenchymal populations. h Dot plot of chondrogenic marker genes (black) and chondrogenic modules (red) expression in '‘fine'’ clusters identified across the three embryonic origins. The number of genes contained in the two modules is indicated in brackets. EC early chondrogenic, PC precursor, IM intramembranous ossification, LS late skeletal. tSNE representations of re-clustered mesenchymal cells of nasal (i), somite (j), and limb (k) origins, with ‘'fine'’ cluster annotations and superimposed streamline plots of scVelo vector fields. Ec clusters are highlighted in red.

Fig. 2. Distinct trans- and cis-regulatory modalities underlie the convergent specification of skeletogenic cells.

a Unsupervised hierarchical clustering and heatmap representation of pairwise Spearman’s rank correlation coefficients of highly variable genes in pseudobulk transcriptomes of ‘'fine'’ mesenchymal clusters. Anatomical origins of pseudobulks are indicated by color code. EC cells (highlighted in red) across embryonic origins cluster together (red dotted square). b Unsupervised hierarchical clustering and heatmap representation of pairwise Spearman’s rank correlation coefficients of expressed transcription factors in pseudobulk transcriptomes of ‘'fine'’ mesenchymal clusters. Anatomical origins of pseudobulks are indicated by color code. All pseudobulks, including EC cells (highlighted in red), cluster by embryonic origins. c Dot plot of embryonic origin-enriched transcription factors' expression in EC cells. d Sampling scheme to assess chromatin accessibility signatures of skeletogenic cells across anatomical locations (color coded) and embryonic stages (Hamburger–Hamilton, HH). tSNE representations of integrated single-cell chromatin accessibility (e) and single-cell transcriptome (f) data across the three anatomical locations. The anatomical origins of cells are indicated by symbols, with the mesenchymal populations used for '‘fine'’ re-clustering highlighted in color. g Unsupervised hierarchical clustering and heatmap representation of pairwise Spearman’s rank correlation coefficients of differentially accessible peaks in pseudobulk chromatin accessibility data of '‘fine'’ mesenchymal clusters. Anatomical origins of pseudobulks are indicated by color code. All pseudobulks, including EC cells (highlighted in red), cluster by embryonic origins. h Coverage plots and heatmap representations of the top 500 differentially accessible peaks (DAPs) in EC cells of neural crest, somatic, and LPM origin. i Promoter depletion and intronic/intergenic elements enrichment of DAPs in EC cells, relative to the consensus peak set. j Boxplots of pairwise Spearman’s rank correlation coefficients of DAPs in EC cells across embryonic origins, calculated using promotor–proximal elements (n = 1023) or an equal number of randomly sampled distal intronic/intergenic elements. Boxplots of 83 pairwise Spearman’s rank correlation coefficients, center line = median, box limits = quartiles, and whiskers = 1.5× interquartile range.

Fig. 3. Trans- and cis-regulatory dynamics of skeletogenic convergence across three embryonic lineages.

tSNE representations of single-cell transcriptomes and co-embedded single-cell chromatin accessibilities (insets) for nasal (a), somite (b), and limb (c) origins. Superimposed on the single-cell transcriptomes are the pseudotime trajectories identified by slingshot, with the chondrogenic trajectories used for further analyses highlighted in red. Pseudotime progression is visualized by heatmaps on scRNA and scATAC data. '‘Fine'’ cluster annotations as in Fig. 1, with EC clusters highlighted in red. d–l Binned pseudobulk dynamics along chondrogenic pseudotime trajectories. d–f Z-score scaled expression dynamics along the chondrogenic pseudotime trajectories for chondrogenic modules (red) and common differentially expressed genes identified in all three embryonic origins (black). g–i Log-transformed expression dynamics of embryonic origin-specific transcription factors identified in nasal (magenta dotted line), somite (blue dotted line), and limb (yellow dotted line) samples. j–l Average normalized peak accessibility dynamics of embryonic origin-specific differentially accessible peaks (DAPs) identified in nasal (magenta dotted line), somite (blue dotted line), and limb (yellow dotted line) samples. Select DAP-adjacent genes are indicated on the right.

Fig. 4. Embryonic origin-specific transcription factor binding motif activities and protein interaction profiles.

a Final number and transcription factor (TF) family distribution of non-redundant de novo identified binding motifs across the three embryonic origins. b Differential cluster-specific motif activities in mesenchymal populations across embryonic origins. EC clusters are highlighted in red. c Embryonic origin-specific TF motif activities (left) and mRNA expression profiles (right), enriched in EC cells. d Venn diagram displaying the overlap of chondrocyte-enriched motif activities across embryonic origins. e Motif activity heatmaps of the ten commonly chondrocyte-enriched TF motifs across embryonic origins. The embryonic origin in which the respective motif was identified is indicated by color-coded circles on the right. EC clusters are highlighted in red. f Venn diagram displaying the overlap of unique chondrocyte-enriched motif activities and commonly enriched chondrogenic genes. g Enriched expression of SOX9, SOX5, and FOXP1 in EC cells of the three embryonic origins. Position weight matrices of a TF motif identified in mesenchymal cells of all embryonic origins (SOX9, h) and a TF motif identified in non-mesenchymal cells (skin) of all embryonic origins (TAF1B, i). Pairwise motif similarity scores are displayed on the right, with wider/darker lines indicating higher similarity. j Boxplots of pairwise motif similarity scores for TFs identified in multiple embryonic origins. Motifs were binned according to which general cell population they were identified in, i.e., mesenchymal (n = 204) versus non-mesenchymal (n = 19). Boxplot center line = median, box limits = quartiles, and whiskers = 1.5× inter-quartile range. k Principal component analysis of replicate SOX9 and FOXP1 RIME (Rapid immunoprecipitation mass spectrometry of endogenous proteins) experiments in NC- (N, magenta), somitic mesoderm- (S, blue), and LPM-derived (L, yellow) tissues. Numbers in brackets indicate the percentage of total variance explained by PC1 and PC2. l UpSet plot of significantly differentially abundant (DA) transcription factors and chromatin modifiers identified pairwise across pairs of embryonic origins for SOX9 and FOXP1 RIME experiments. m Averaged log2-normalized intensities for SOX9-specific (left) and FOXP1-specific (right) DA transcription factors and chromatin modifiers across embryonic origins.

Fig. 5. Lineage-specific enhancer-promoter interactions of core chondrogenic genes.

a-c hkmeans-clustered peak-to-gene link heatmaps displaying Z-score normalized peak accessibilities (left, ocher) and imputed target gene expression levels (right, gray) in single cells (columns). Cell clusters dominated by EC cells are highlighted in red (top). The top 40 peak-to-gene links are shown. Total numbers of identified peak-to-gene links per cluster are indicated by barplots on the right, cluster numbers (kX) on the left. d Venn diagram displaying the overlap of mesenchymal peak-to-gene link CREs across embryonic origins (top) and the overlap of target genes contained within these peak-to-gene links (bottom). e Average evolutionary conservation (calculated as phastCons scores) across CREs in peak-to-gene links shared amongst all embryonic origins (n = 40), shared by more than one embryonic origin (n = 886), and that are origin-specific ones (nasal n = 4396, somite n = 8452, and limb n = 6975). Boxplot center line = median, box limits = quartiles, and whiskers = 1.5× inter-quartile range. f Hierarchically clustered peak-to-gene link correlation heatmap of links identified at core chondrogenic genes across embryonic origins. Select target genes are indicated on the right.

Fig. 6. Distinct cis- and trans-regulatory dynamics during the lineage-specific activation of the core chondrogenic transcription factor SOX9.

a Pseudobulk aggregate genome accessibility tracks at the SOX9 locus. For each embryonic origin, tracks for PC populations (light colors) and EC cells (dark colors) are displayed. Identified peaks and neighboring genes are indicated below. b Relative transcription factor binding motif abundances in DAPs at the SOX9 locus in cells of nasal (left), somite (middle), and limb origin (right). Genome-wide dynamics of motif activities (c) and the corresponding RNA expression profiles (d) of transcription factors identified in (b), following chondrogenic pseudotime trajectories in nasal (top), somite (middle), and limb cells (bottom). For reference, both motif activities and RNA expression profiles of SOX9 are included. e embryonic origin-specific peak-to-gene link plots at the SOX9 locus. f Pseudobulk aggregate genome accessibility tracks of peaks tested in enhancer reporter assays in (g). Candidate enhancers are highlighted by gray dotted boxes. PC population tracks in light colors, EC tracks in dark colors. g In vivo enhancer reporter assays. Chondrocyte candidate enhancer elements with predicted embryonic origin-specificity were cloned into reporter constructs driving GFP expression and electroporated (EP) together with a tdTomato control plasmid into cells of the cranical neural crest (cNC), the somites, or the forelimb LPM. Embryonic origin-specificity of GFP expression in SOX9-positive chondrogenic condensations is indicated by green squares. The candidate ‘'nasal'’ enhancer showed a weak GFP signal in the limb mesenchyme as well (green dotted square, asterisk).

Fig. 7. Model for the transcriptional convergence of skeletogenic cells of different embryonic origins.

Depending on anatomical location, different embryonic lineages give rise to the skeletal elements of the vertebrate skeleton. These mesenchymal PC express distinct lineage-specific transcription factor profiles, according to their embryonic origin, which are integrated at lineage-specific enhancer elements to activate a shared set of genes belonging to a core transcriptional program of vertebrate skeletal progenitors.