Distinct Gene Regulatory Dynamics Drive Skeletogenic Cell Fate Convergence During Vertebrate Embryogenesis

Distinct Gene Regulatory Dynamics Drive Skeletogenic Cell Fate Convergence During Vertebrate Embryogenesis | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Distinct Gene Regulatory Dynamics Drive Skeletogenic Cell Fate Convergence During Vertebrate Embryogenesis Menghan Wang , Ana Di Pietro-Torres , View ORCID Profile Christian Feregrino , Maëva Luxey , View ORCID Profile Chloé Moreau , View ORCID Profile Sabrina Fischer , View ORCID Profile Antoine Fages , View ORCID Profile Patrick Tschopp doi: https://doi.org/10.1101/2024.03.26.586769 Menghan Wang 1 DUW Zoology, University of Basel , 4051 Basel, Switzerland Find this author on Google Scholar Find this author on PubMed Search for this author on this site Ana Di Pietro-Torres 1 DUW Zoology, University of Basel , 4051 Basel, Switzerland Find this author on Google Scholar Find this author on PubMed Search for this author on this site Christian Feregrino 1 DUW Zoology, University of Basel , 4051 Basel, Switzerland Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Christian Feregrino Maëva Luxey 1 DUW Zoology, University of Basel , 4051 Basel, Switzerland Find this author on Google Scholar Find this author on PubMed Search for this author on this site Chloé Moreau 1 DUW Zoology, University of Basel , 4051 Basel, Switzerland Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Chloé Moreau Sabrina Fischer 1 DUW Zoology, University of Basel , 4051 Basel, Switzerland Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Sabrina Fischer Antoine Fages 1 DUW Zoology, University of Basel , 4051 Basel, Switzerland Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Antoine Fages Patrick Tschopp 1 DUW Zoology, University of Basel , 4051 Basel, Switzerland Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Patrick Tschopp For correspondence: patrick.tschopp{at}unibas.ch Abstract Full Text Info/History Metrics Preview PDF Abstract Cell type repertoires have expanded extensively in metazoan animals, with some clade-specific cells being paramount to their evolutionary success. A prime example are the skeletogenic cells of vertebrates that form the basis of their developing endoskeletons. Depending on anatomical location, these cells originate from three different embryonic precursor lineages – the neural crest, the somites, and the lateral plate mesoderm – yet they converge developmentally towards similar cellular phenotypes. Furthermore, these lineages have gained ‘skeletogenic competency’ at distinct timepoints during vertebrate evolution, thus questioning to what extent different parts of the vertebrate skeleton rely on truly homologous cell types. Here, we investigate how lineage-specific molecular properties of the three precursor pools are integrated at the gene regulatory level, to allow for phenotypic convergence towards a skeletogenic cell fate. Using single-cell transcriptomics and chromatin accessibility profiling along the precursor-to-skeletogenic cell continuum, we examine the gene regulatory dynamics associated with this cell fate convergence. We find that distinct transcription factor profiles are inherited from the three precursor states, and that lineage-specific enhancer elements integrate these different inputs at the cis -regulatory level, to execute a core skeletogenic program. We propose a lineage-specific gene regulatory logic for skeletogenic convergence from three embryonic precursor pools. Early skeletal cells in different body parts thus share only a partial ‘deep homology’. This regulatory uncoupling may render them amenable to individualized selection, to help to define distinct morphologies and biomaterial properties in the different parts of the vertebrate skeleton. Introduction Metazoan bodies are defined by obligate multicellularity and the presence of functionally, morphologically, and molecularly distinct cell types 1 - 3 On a developmental timescale, this requires the implementation of distinct gene regulatory programs, to generate cellular diversity from a single fertilized cell, the zygote, that inherits a common genomic blueprint to its cellular descendants 4 - 6 . Such regulatory programs can also inform us about evolutionary relationships amongst different cell types, both within and across species boundaries. Namely, on an evolutionary timescale, novel cell types are often thought to arise through the genomic re-use of pre-existing regulatory modules, via re-shuffling, neo-functionalization or co-option. Accordingly, similarities or dissimilarities in the gene regularity dynamics underlying their developmental specification can help to resolve potential cell type homologies 1 , 7 . In this context, vertebrate skeletogenesis offers a relevant opportunity to study the gene regulatory dynamics of cell fate specification across both developmental and evolutionary timescales. A cartilaginous and often ossified endoskeleton is a hallmark of the vertebrate clade, with many accessory connective tissue types required to build a functional skeleton 8 . Developmentally, the specification of the cells required to build these tissues initiates with mesenchymal precursors condensing, at the onset of vertebrate skeletogenesis, and then progresses through distinct, cell type-specific differentiation processes 9 - 12 . A peculiarity in this process is that – unlike as for many other cell lineages ( Fig.1A , left) – these multipotent skeletal progenitors can arise convergently, from distinct embryonic sources ( Fig.1A , right). Namely, depending on anatomical location, three distinct embryonic lineages – the cranial neural crest, the somitic sclerotome, and the somatopleure of the lateral plate mesoderm – are the developmental origin of the cranial, axial and appendicular skeleton, respectively 8 ( Fig. 1B ). The early specification of these cells thus needs to integrate discrete molecular states, inherited from their respective embryonic precursor sources, to facilitate a skeletogenic cell fate convergence. Download figure Open in new tab Figure 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 (C1-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 variably 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. ( I - K ) 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. Early chondrogenic clusters are highlighted in red. The distinct embryonic trajectories of skeletogenic cells are also reflecting the evolutionary histories of the different parts of the vertebrate skeleton. The first vertebrate skeletal elements are thought to have originated in the head region, mirroring the prevertebrate presence of cartilage-like structures supporting a feeding apparatus 13 . Subsequently, the ability to form a progenitor cell-based endoskeleton expanded along the primary and secondary body axes, giving rise to structurally supportive yet flexible elements in the axial and appendicular skeletons, respectively 14 . However, these distinct evolutionary and developmental histories may also challenge our understanding to what extent the convergently specified cells of the vertebrate skeleton can be considered truly homologous 15 - 17 . Here, using single-cell functional genomics along the three distinct mesenchymal precursor-to-skeletogenic cell trajectories, we investigate the genome-wide regulatory dynamics in a vertebrate embryo at cellular resolution. We provide evidence that lineage-specific transcription factor profiles are inherited from the respective embryonic origins and that these are integrated at distinct cis -regulatory elements, to canalize developmental cell fate trajectories towards an early skeletogenic convergence point. These distinct dynamics imply a regulatory uncoupling between skeletogenic cells at different anatomical locations. We discuss the resulting implications for cell type homology assessment in the vertebrate skeleton, and the potential of distinct evolutionary trajectories in skeletal cell and tissue properties, upon co-optiondependent convergence of gene regulatory programs across embryonic lineages. Results A convergent transcriptomic signature in skeletogenic cells of different embryonic origins To establish the temporal progression of skeletogenic initiation and maturation across the three anatomical locations and embryonic lineages, we first performed chromogenic in situ hybridizations (ISH) on a developmental time series of chicken embryos. We used cranial sagittal cryosections covering the frontonasal prominence, and brachial trunk transversal sections including the emerging forelimb buds. We investigated the expression of SOX9 , an early marker of skeletogenic induction 9 , and Aggrecan ( ACAN ), an extracellular matrix protein of more mature skeletal cells 18 ( Fig.S1A,B ). Based on the observed expression dynamics, we devised a sampling strategy following the cell- and tissue-specific transcriptional changes of the three embryonic lineages, from mesenchymal precursors towards the onset and maturation of skeletogenic tissues, using 10x Chromium single-cell RNA-sequencing (scRNA-seq) profiling ( Fig.1C ). In total, we obtained over 23’000 high quality single-cell transcriptomes. For each of the three anatomical locations, we integrated the three sampled embryonic stages and performed tSNE non-linear dimensionality reduction and graph-based clustering using Seurat 19 ( Fig.1D-F ). These ‘broad’ clusters were annotated with the help of expression profiles of known marker genes ( Fig.S1C,E ). All samples showed comparable transcriptome complexities, and contributed similarly to the different clusters ( Fig.S1F,G ). To assess transcriptomic similarities amongst these clusters, we next generated cluster-based pseudobulks and calculated Spearman’s rank correlation coefficients on differentially expressed genes across all cell types and anatomical locations. Unsupervised hierarchical clustering revealed that cell types originating from the same embryonic lineage, but sampled at different anatomical locations – like, e.g., skin or blood cells – showed highly similar transcriptional profiles, in agreement with their shared developmental history ( Fig.1G ). Intriguingly, however, our analysis revealed that also mesenchymal cells stemming from discrete embryonic lineages – that is, from the neural crest, the somites, or the lateral plate mesoderm – clustered together, indicating a transcriptional convergence amongst them ( Fig.1G , red dotted square). Across our three anatomical sampling sites, we focused on ‘broad’ clusters that likely contained cells transitioning from a mesenchymal precursor state towards a skeletogenic fate ( Fig.1D-F , color coded), and re-clustered them at finer resolutions. Within these ‘fine’ clusters, we again used expression profiles of known marker genes, as well as two previously identified early chondrogenic gene co-expression modules, ‘IMM’ 20 and ‘RED’ 21 , to assess their respective skeletogenic differentiation. In each embryonic lineage, we identified a ‘fine’ cluster enriched for an early skeletogenic, i.e. chondrogenic, signature ( Fig.1H , highlighted in red). For the limb sample, we additionally isolated a cluster showing signs of more mature skeletal cells ( Fig.1H , limb ‘fine’ cluster 5). Using scVelo 22 , we approximated cell fate transition trajectories in silico and projected the predicted vector fields using streamline plots on tSNE representations of our mesenchymal samples. For all three anatomical locations, scVelo predicted trajectories connecting the mesen-chymal precursors – as identified by known marker genes – to early chondrogenic populations ( Fig.1I-K ). Using ISH and single-cell transcriptomics, our analyses revealed distinct temporal dynamics, but converging molecular signatures amongst mesenchymal cells with skeletogenic potential across the three embryonic lineages. Furthermore, our scVelo analysis suggested that our single-cell transcriptomics data captured the entire specification spectrum, from uncommitted mesenchymal precursor cell to early chondrocyte. Distinct trans - and cis -regulatory modalities underlie the convergent specification of skeletogenic cells To detail the transcriptional signatures underlying the switch from uncommitted mesenchymal precursor cell to early chondrocyte, we focused our analyses on the ‘fine’ mesenchymal clusters with skeletogenic potential. We re-assessed transcriptional similarities of these sub-populations using differentially expressed genes and Spearman’s rank correlation coefficients of pseudobulk transcriptomes. Again, we found high transcriptional similarities amongst early skeletogenic cells at the three anatomical locations, as indicated by them clustering together ( Fig.2A , red dotted square). Download figure Open in new tab Figure 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 variably genes in pseudobulk transcriptomes of ‘fine’ mesenchymal clusters. Anatomical origins of pseudobulks are indicated by color code. Early chondrogenic 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 early chondrogenic cells (highlighted in red), cluster by embryonic origins. ( C ) Dot plot of embryonic originenriched transcription factors expression in early chondrogenic cells. ( D ) Sampling scheme to assess chromatin accessibility signatures of skeletogenic cells across anatomical locations (color coded) and embryonic stages (Hamburger-Hamilton, HH). ( E - F ) tSNE representations of integrated single-cell chromatin accessibility ( E ) and single-cell transcriptome ( F ) data across the three anatomical locations. Anatomical origins of cells are indicated by symbols, with the mesenchymal populations used for ‘fine’ reclustering 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 early chondrogenic cells (highlighted in red), cluster by embryonic origins. ( H ) Coverage plots and heatmap representations of the top500 differentially accessible peaks (DAPs) in early chondrogenic cells of neural crest, somitic and lateral plate mesoderm origin. ( I ) Promoter depletion and intronic/intergenic elements enrichment of DAPs in early chondrogenic cells, relative to the consensus peak set. ( J ) Boxplots of pairwise Spearman’s rank correlation coefficients of DAPs in early chondrogenic cells across embryonic origins, calculated using promotor-proximal elements or an equal number of randomly sampled distal intronic/intergenic elements. To investigate potential upstream regulatory inputs controlling these similar transcriptomes, we next focused our attention on transcription factors. Spearman’s rank correlation coefficients on transcription factor expression profiles, however, re-clustered the mesenchymal sub-populations strictly by embryonic origins, irrespective of their skeletogenic differentiation state ( Fig.2B ). This suggested that the mesenchymal precursor populations carried over a lineage-specific repertoire of expressed transcription factors, while undergoing skeletogenic induction. Indeed, when looking at transcriptional regulators enriched in chondrogenic cells of the three anatomical locations, we find clear evidence of a lineage-specific heritage of expressed transcription factors, many of which have known developmental functions in their respective anatomical locations ( Fig.2C ). Thus, counterintuitively, this indicated that – across anatomical locations and embryonic lineages – overall similar transcriptional signatures were generated with distinct upstream trans -regulatory inputs, i.e. lineage-specific transcription factor expression profiles. We reasoned that these distinct trans -regulatory inputs could potentially be integrated at the cis -regulatory level of key skeletogenic genes, to facilitate transcriptomic convergence. To investigate this possibility, we performed single-cell chromatin accessibility assays (single-cell assay for transposase-accessible chromatin with sequencing, or scATAC-seq 23 ) to identify genomic elements with potential regulatory activity. We followed a similar sampling scheme as for our scRNA-seq approach, although – reasoning that such cis -regulatory recoding would occur during early stages of skeletogenic induction – we excluded the late time points ( Fig.2D ). In total, we obtained over 36’000 cells with high quality chromatin accessibility profiles and, using MACS2 24 , we identified a consensus set of 678,707 peaks across the three anatomical locations. We integrated the two embryonic stages per sampling site, performed non-linear dimensionality reduction and identified potential cell type-specific clusters ( Fig.S2A,C ). We annotated these clusters with the help of our scRNA-seq data, using label transferring and non-negative least squares (NNLS) regression, as well as visual inspection of scATAC-seq ‘marker peaks’ (see Material and Methods, and Fig.S2D). At all three anatomical locations, overall similar cell type repertoires were recovered as in our scRNA-seq sampling. We then combined all samples across the three anatomical locations, for both scRNA-seq and scATAC-seq data sets, and performed anchor-based integration. Mesenchymal cells in our scATAC-seq data visually appeared to intermingle less than in our scRNA-seq data ( Fig.2E,F ). This may indicate the presence of embryonic origin-specific states in chromatin accessibility, as opposed to the convergent signatures at the transcriptomic level ( Fig.1G ). To investigate potential lineage-specific chromatin accessibilities, we re-clustered the mesenchymal scATAC-seq cell populations with skeletogenic potential and annotated the resulting ‘fine’ clusters using our scRNA-seq data ( Fig.S2E,G ). We then performed unsupervised hierarchical clustering on Spearman’s rank correlation coefficients of differentially accessible peaks (DAPs) in pseudobulks of these ‘fine’ mesenchymal clusters. Akin to our scRNA-seq analyses of transcription factor expression profiles, this resulted in a strict embryonic origin-dependent clustering of mesenchymal populations, including early chondrogenic cells ( Fig.2B,G ). Major cell type classifications, however, still clustered according to their ‘broad’ annotations, with similar chromatin accessibility signatures echoing their shared embryonic origins ( Fig.S2H ). This further implied that mesenchymal and skeletogenic cells across the three anatomical locations carry distinct chromatin accessibility profiles, reflecting their discrete embryonic origins and lineage histories. Indeed, coverage plots across the three anatomical locations revealed that a substantial fraction of DAPs within the respective chondrogenic populations are distinct and hence embryonic origin-specific ( Fig.2H ). We classified these DAPs according to their genomic location, with respect to the transcription start sites of neighboring genes. Interestingly, we found that promoter-proximal peaks were depleted in our DAP set, relative to the consensus gene set, while more distal elements – intronic and intergenic – appeared enriched ( Fig.2I ). This implied that the majority of differences in chromatin accessibilities, between chondrogenic cells from different embryonic lineages, originated at distally located peaks. Indeed, DAPs at promoter-proximal peaks showed on average higher Spearman’s rank correlation coefficients across embryonic origins than distal intronic/intergenic ones ( Fig.2J ). This suggested that the similar transcriptional profiles observed at the RNA level originate from similar promoter repertoires. Distal peaks, however, where one would anticipate putative long-range enhancer elements to be located, showed much lower similarity in accessibilities amongst the different embryonic origins ( Fig.2J ). Collectively, our combined scRNA-seq and scATAC-seq approach revealed the presence of distinct trans - and cis -regulatory signatures in skeletogenic cells at the three anatomical locations, as evidenced by embryonic origin-specific transcription factor expression profiles and the presence of discrete chromatin accessibility signatures at distal locations. Trans - and cis -regulatory dynamics of skeletogenic convergence across the three embryonic lineages To follow the trans - and cis -regulatory changes underlying these convergent cell fate transitions, we next investigated scRNA-seq and scATAC-seq dynamics along chondrogenic pseudotime trajectories in the three anatomical locations. Using tSNE embeddings of our scRNA-seq data, we constructed minimum spanning trees using slingshot 25 , with preset start points corresponding to the clusters with naïve mesenchymal expression signatures ( Fig.3A-C ). We observed overall similar trajectory predictions as for our scVelo analysis ( Fig.1H-J ). For each anatomical location, we were able to retrieve a single trajectory either ending in an early chondrogenic cluster ( Fig.3A,B ), or traversing it towards more mature skeletal cells ( Fig.3C ). We transferred the pseudotime values of our chondrogenic scRNA-seq trajectories to our scATAC-seq data, to follow the accompanying chromatin accessibility dynamics ( Fig.3A-C , insets). Download figure Open in new tab Figure 3. Trans - and cis -regulatory dynamics of skeletogenic convergence across three embryonic lineages. ( A - C ) 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 visualized by heatmaps on scRNA and scATAC data. ‘Fine’ cluster annotations as in Figure 1 , with early chondrogenic 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. We then binned the respective scRNA-seq pseudotimes into equidistant pseudobulks and used TrAGEDy 26 to align the chondrogenic gene expression dynamics along the respective trajectories of the three embryonic origins. We found overall higher similarities towards the ends of these pair-wise comparisons, indicating chondrogenic convergence at the transcriptional level ( Fig.S3A,C ). Both nasal and somite trajectories failed to align to the last section of our limb chondrogenic pseudotime. This corroborated the notion that we had recovered more mature skeletal cells in our limb samples, compared to nasal and somite ( Fig.S3A,F , Fig.1H ). Accordingly, we excluded the corresponding limb pseudotime bins containing these cells from further analyses. We first checked the expression dynamics of the two chondrogenic gene co-expression modules ‘IMM’ and ‘RED’. Both modules increased in their expression along the three pseudotime trajectories, in nasal, somite and limb samples ( Fig.3D-F , red). Next, we looked at the expression dynamics of a core set of common chondrogenic genes, found to be enriched in chondrocytes across the three embryonic lineages ( Fig.S4A ). All genes showed an increase in expression along the respective trajectories ( Fig.3D-F , black). Contained within this shared set of genes were known chondrogenic regulators and extracellular matrix proteins, as well as a selection of ribosomal proteins ( Fig.3D-F , Fig.S4A ). Additionally, using tradeSeq 27 , we identified lineage-specific expression dynamics of known and novel regulators activated in chondrogenesis across the three anatomical locations ( Fig.S4B,D ). In all three lineages, the chondrogenic wave appeared to be preceded by an increase in expression of origin-specific transcription factors ( Fig.3G-I , Fig.2C ). Finally, we found evidence for distinct chromatin accessibility dynamics. Chondrocyte-specific DAPs showed increased accessibilities along the respective chondrogenic scATAC pseudotime trajectories, but in an embryonic origin-specific manner ( Fig.3J-L ). Using integrative scRNA-seq and scATAC-seq pseudotime analyses, we detailed the emergence of common transcriptional signatures and distinct trans - and cis -regulatory profiles during skeletogenic convergence. Namely, while increased expression of chondrogenic modules and a core set of differentially expressed genes was shared across the three embryonic origins, these were accompanied by lineage-specific transcription factor expression and chromatin accessibility dynamics. Embryonic origin-specific activities and specificities of transcription factor binding motifs To investigate the interplay of distinct transcription factor profiles and origin-specific chromatin accessibilities, we next evaluated cell type-specific activities of transcription factor binding motifs. Given the scarcity of publicly available and experimentally validated binding motifs for chicken transcription factors, we decided to define our own set of DNA position weight matrices. Briefly, we used Homer 28 to identify enriched de novo motifs in our scATAC-seq data in a cluster-by-cluster manner, across the three anatomically distinct samples. We annotated these de novo motifs with candidate transcription factors using public repositories, and selected the best matches based on motif similarity and the correlation of motif activities with scRNA-seq transcription factors expression profiles (see Material and Methods). In total, we identified and annotated 1373 de novo motifs across the three anatomical locations ( Fig.S5A ), with 540 non-redundant ones used for the further analyses ( Fig.4A ). Motifs for members of the homeobox, C2H2 zinc fingers and basic helix-loop-helix (bHLH) protein transcription factor families were amongst the most frequently identified ones ( Fig.4A ). Furthermore, our limb mesenchyme de novo motif for SOX9 matched a chick limb ChIP-seq 29 validated motif more closely than publicly available position weight matrices ( Fig.S5B ). Encouraged by this, we continued all subsequent analyses with our de novo motifs only. Download figure Open in new tab Figure 4. Embryonic origin-specific de novo prediction and activity profiles of transcription factor binding motifs. ( A ) Final numbers and transcription factor (TF) family distribution of non-redundant de novo binding motifs across the three embryonic origins. ( B ) Differential cluster-specific motif activities in mesenchymal populations across embryonic origins. ‘Early chondrogenic’ clusters are highlighted in red. ( C ) Embryonic origin-specific TF motif activities (left) and mRNA expression profiles (right), enriched in early chondrogenic 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. ‘Early chondrogenic’ 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, FOXP1 in early chondrogenic cells of the three embryonic origins. ( H, I ) 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 versus non-mesenchymal. Using this custom set of position weight matrices, we next conducted differential motif activity analyses. At both ‘broad’ and ‘fine’ cluster resolutions, we identified motif activity signatures that were predictive for specific cell types, including early chondrocytes and mesenchymal precursor cells ( Fig.4B and Fig.S5C,D ). Moreover, we identified motif activities in chondrogenic cells with high embryonic origin-specificity. These were mirrored by expression signatures of the corresponding transcription factors. This emphasized the presence of distinct trans -regulatory inputs in chondrogenic cells of the three embryonic lineages, at both RNA and motif activity levels ( Fig.4C ). Of all motifs whose activities were enriched in chondrogenic cells, only ten were shared across the three anatomical locations ( Fig.4D,E , Fig.S5E ). Comparing them to the core chondrogenic genes identified in our differential expression analyses revealed only three genes showing consistent chondrocyte enrichment at both RNA expression and motif activity levels: SOX9, SOX5 and FOXP1 ( Fig.4F,G ). With the length range of our Homer de novo motif search (8-22 bp) we were able to investigate the occurrence of potential co-binding patterns of multiple transcription factors. Sequences bound by multimeric protein complexes are increasingly recognized as an integral aspect of DNA’s regulatory grammar, to provide robustness or diversify activity patterns in a combinatorial manner 30 - 32 . Indeed, many of our motifs showed a bimodal distribution of nucleotide enrichment in their position weight matrices, indicative of dimers binding to them. For example, the SOX9 motifs identified in somite and limb mesenchymal cells indicated a homodimer-like binding, in agreement with previous findings 29 ( Fig.4H and Fig.S5B ). Additionally, in nasal samples, our analyses predicted a motif showing a SOX9-like monomer signature at its 3’ end, but a potential heterodimeric binding partner at its 5’ end ( Fig.4H ). To investigate this phenomenon more systematically – that is, the same transcription factor having dissimilar binding motifs predicted in the three different embryonic lineages – we trimmed the extremities of our position weight matrices based on minimal nucleotide enrichment scores and calculated motifs similarities across anatomical locations. In our analysis, we split the motifs based on whether their position weight matrices were identified in mesenchymal cells of different embryonic origins, or in non-mesenchymal populations of the same embryonic lineage (e.g. skin, Fig.4I ). On average, position weight matrices for transcription factors identified in mesenchymal cells showed lower similarity scores than motifs identified in non-mesenchymal cell types ( Fig.4J ). This indicated that binding motifs for the same transcription factors were less conserved in mesenchymal cells of different embryonic origins, potentially due to lineage-specific differences in their cells’ chromatin environment, or distinct co-factors binding the motifs in a heterodimeric manner ( Fig.4H-J ). Collectively, we identified de novo transcription factor binding motifs, some of which showed common cell type-specific activities while others were embryonic lineage-restricted. Furthermore, we find partially diverging sequence logos for the same transcription factors in mesenchymal cells at different anatomical locations. This substantiates our previous findings of lineage-specific trans -regulatory inputs during the skeletogenic convergence of different embryonic precursor lineages at the level of motif architectures and activities. Lineage-specific regulatory architectures and enhancer activities in skeletogenic cells We used ArchR 33 on our mesenchymal scATAC-seq data for peak-to-gene links analyses, to connect putative distal enhancer elements to target genes and test for lineage-specific activities. Based on a minimal correlation coefficient and adjusted p-value cutoffs, we identified over 28’000 presumptive peak-to-gene links ( Fig.5A-C and Fig.S6A,C ). Within these peak-to-gene links, target genes were generally predicted to be contacted by less than three putative cis -regulatory elements (CRE), and the majority of CREs was interacting with only one target gene ( Fig.S6D,F ). Using hierarchical k-means (hkmeans) clustering, we sorted our peak-to-gene links according to their aggregate activity profiles and plotted z-score normalized chromatin accessibilities and imputed target gene expression levels ( Fig.5A-C and Fig.S6G,I ). For both somite and limb samples, and to some extent for nasal cells, we identified clusters enriched for early chondrogenic cells ( Fig.5A-C ), with corresponding enrichments for skeletogenesis-related terms ( Fig.S6G,I ). Download figure Open in new tab Figure 5. Lineage-specific regulatory architectures and enhancer activities in skeletogenic cells. ( A - C ) hkmeans-clustered peak-to-gene link heatmaps displaying Z-score normalized peak accessibilities (left, ochre) and imputed target gene expression levels (right, grey) in single cells (columns). Cell clusters dominated by early chondrogenic cells are highlighted in red (top). The top40 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 (left) and overlap of target genes contained within these peak-to-gene links (right). ( E ) 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. ( F ) Pseudobulk aggregate genome accessibility tracks at the SOX9 locus. For each embryonic origin, tracks for precursor populations (light colors) and early chondrogenic cells (dark colors) are displayed. Identified peaks and neighboring genes are indicated. Below, embryonic origin-specific peak-to-gene link plots. ( G ) Pseudobulk aggregate genome accessibility tracks of peaks tested in enhancer reporter assays in H . Coordinates of candidate enhancers are highlighted by grey dotted boxes. Precursor population tracks in light colors, early chondrogenic tracks in dark colors. ( H ) 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 lateral plate mesoderm (LPM). Embryonic origin-specificity of GFP expression in SOX9-positive chondrogenic condensations is indicated by green squares. The candidate ‘nasal’ enhancer showed weak GFP signal in the limb mesenchyme as well (green dotted square, asterisk). Of all peak-to-gene links, only 41 CRE coordinates were detected in all three anatomical locations using our stringent cutoffs. However, that number increased to 607, if we were only focusing on the overlap of predicted target genes therein ( Fig.5D ). This suggested that many commonly expressed genes were contacted by lineage-specific enhancer elements. Indeed, looking at the correlation scores of our core chondrogenic gene set ( Fig.S4A , excluding ribosomal proteins) revealed largely lineage-specific peak-to-gene link activities ( Fig.5E ). We further explored this in more detail at the SOX9 locus, a well-known regulator of skeletogenic induction 9 . While a high density of consensus peaks was present in the genomic landscape flanking that gene, the overall pseudobulk chromatin accessibility profiles looked distinct from one anatomical location to another, both in mesenchymal precursors (bright colors) as well as at the early chondrogenic stage (dark colors) ( Fig.5F , top). Importantly, link plots at the locus revealed that many of the peak-to-gene links were in fact predicted to be embryonic origin-specific ( Fig.5F , bottom). To evaluate putative enhancer functions of these peaks and test for their lineage-specificity, we first used differential chromatin accessibility analysis to define chondrocyte-enriched DAPs at the SOX9 locus. For all three anatomical locations, we identified peaks with high chromatin accessibilities in early chondrogenic cells, relative to mesenchymal precursors and chondrocytes of other embryonic origins ( Fig.5G ). We isolated the corresponding sequences from genomic DNA and cloned them into reporter plasmids, upstream of a minimal promoter driving green fluorescent protein (GFP) expression. We electroporated the resulting constructs into the precursor populations of the respective skeletogenic lineages in ovo . As electroporation control, we included a plasmid containing a strong constitutively active promoter driving tdTomato expression. Post-electroporation, we let the embryos develop further for two days, harvested the targeted tissues and processed them for histology. We performed immunohistochemistry against endogenous SOX9 protein, to determine the location of chondrogenic condensations, and against tdTomato and GFP, to evaluate electroporation efficiency and enhancer activity, respectively. Within electroporated condensations, we scored the presence or absence of GFP signal, to indicate chondrogenic enhancer activity. The predicted somite and limb enhancers showed specific GFP reporter activity, restricted to the respective embryonic lineages they were identified in ( Fig.5H , green boxes). The nasal candidate enhancer reporter drove strong GFP signal in cranial neural crest-derived chondrogenic tissue ( Fig.5H , green box). It additionally did so at reduced levels in limb condensations as well ( Fig.5H , green dotted box and asterisk), suggesting only partial lineage-specificity of this element. Overall, our peak-to-gene link analyses and enhancer reporter assays uncovered the presence of lineage-specific CREs contacting an overlapping set of target genes across the three anatomical locations. This suggested a partial lineage dependency in cis -regulatory interactions of target genes in early chondrogenic cells of different embryonic origins, mirroring their previously established specificity at the trans -regulatory level. Discussion Cell type specification in animals relies on the execution of distinct gene regulatory programs during embryonic and post-embryonic development. Consequently, cell type evolution depends on the origination of new regulatory modalities, to drive innovations in cellular form and function. Here we have documented the gene regulatory dynamics underlying the embryonic specification of skeletogenic cells in vertebrates, an iconic cell type central to their evolutionary success. Following their first evolutionary appearance at the base of vertebrates, additional embryonic lineages have subsequently acquired the ability to form skeletal cell types. This makes the underlying gene regulatory programs interesting from both evolutionary and developmental perspectives. Namely, how the three embryonic lineages have acquired skeletogenic competency at a gene regulatory level during vertebrate evolution, and, developmentally, how these lineages are transcriptionally recoded during the earliest steps of skeletogenesis, to converge from distinct molecular profiles of their origin populations towards a similar skeletogenic phenotype ( Fig.1A ). Transcriptional recoding and lineagememory in mesenchymal cells Using scRNA-seq profiling, we have followed the early transcriptional dynamics of vertebrate skeletogenesis in three distinct embryonic lineages, at three different anatomical locations: cranial neural crest (frontonasal prominence), somitic sclerotome (axial skeleton) and lateral plate mesoderm (appendicular skeleton) ( Fig.1B ). We find shared transcriptional signatures between both mesenchymal and skeletal cells of the three embryonic origins ( Fig.1G ). This aligns with previous observations delineating a ‘mesenchymal signature’ as one of the major transcriptional programs present in primary cells from different tissue types, across developmental lineages 34 . Transitioning towards such a mesenchymal signature may transcriptionally prime different precursor lineages to increase their cell fate plasticity, akin to what is observed in various types of metastatic cancers, while still maintaining partially distinct transcription and chromatin memories of their embryonic origins 35 , 36 . Indeed, skeletogenic cells of the three lineages investigated here all start out as embryonic epithelia and undergo an epithelial-to-mesenchymal transition (EMT), before migrating to their respective anatomical locations in the periphery 37 - 39 . There, they form mesenchymal condensations and appear to initiate a core skeletogenic program which – at a global scale – shows similarity in its overall gene regulatory state 20 , 40 , 41 ( Fig.1H ). We also noted a shared increase in ribosomal protein transcription, potentially to meet the translational needs for increased extracellular matrix protein production and secretion 42 ( Fig.S4A ). Importantly, however, our cellularly resolved analyses along the mesenchymal precursor-to-skeletal progenitor trajectories also uncovered embryonic lineage-specific trans - and cis -regulatory dynamics during the formation of the vertebrate skeleton. Distinct trans - and cis -regulatory dynamics along skeletogenic differentiation trajectories of different embryonic origins We find evidence for a lineage-specific heritage in transcription factor expression profiles that carries over to skeletogenic cells at different anatomical locations of the vertebrate skeleton ( Fig.2B,C ). This implies that the convergent specification of functionally analogous and transcriptionally similar skeletal cell types can be induced by distinct upstream trans -regulatory inputs, even across germ layers. Similarly distinct inputs with convergence were identified between different neuronal lineages in the developing visual system of Drosophila 43 . Here, using scATAC-seq data, we add an additional layer of regulatory information, and demonstrate that distinct chromatin accessibility signatures are accompanying these specific trans -regulatory inputs ( Fig.2G,H ). These lineage-specific cis -regulatory signatures separate unequally between gene-proximal and -distal sites, with promoter accessibilities exhibiting overall higher similarities in skeletal cells across embryonic origins than distal sites with putative enhancer functions ( Fig.2I,J ). Using pseudotime analyses we followed the trans - and cis -regulatory dynamics along the skeletogenic differentiation trajectories of the three embryonic origins and identify both known and novel signatures of this transcriptional and phenotypic cell fate convergence ( Fig.3D-I , Fig.S4A,D ). Intriguingly, the expression levels of lineage-specific trans regulators along our pseudotemporal orderings peak right before the onset of the core skeletogenic program ( Fig.3D-I ), thus raising the possibility of lineage-specific transition states in which mesenchymal precursors are transcriptionally recoded towards a similar skeletogenic cell fate 44 . Thereafter, known skeletal regulators become transcribed in the respective convergence trajectories. We provide evidence that the activation of these genes relies on non-overlapping set of enhancer elements, with lineage-specific chromatin accessibility profiles and promoter-enhancer interactions ( Fig.3J-L , Fig.5E,F ). Using in vivo reporter assays, we demonstrate the lineage-specific activity of distal enhancer elements at a core chondrogenic factor ( Fig.5G,H ). Human genetics and molecular studies using in vivo and in vitro models have previously documented the presence of large regulatory landscapes at skeletogenic genes 45 - 47 . We argue that this regulatory strategy enables distinct up-stream trans inputs to activate a common downstream program in the different embryonic lineages, through integration at the level of cis elements ( Fig.6 ). In this scenario, the implementation of distinct long-range enhancer actions would transcriptionally recode distinct precursor signatures towards a similar skeletogenic cell fate. Download figure Open in new tab Figure 6. 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 precursors 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. Across the different embryonic lineages, it is unlikely that this convergence involves only a single skeletogenic ‘master regulator’. Indeed, SOX9 – a prime candidate for such function – acts in many non-skeletogenic lineages 48 , and its overexpression is unable to fully reprogram cells with skeletogenic potential towards a chondrocyte fate by itself 49 . Rather, an entire battery of transcriptional regulators – some shared, many lineage-specific – seems to drive the lineage-specific convergence towards a skeletogenic phenotype, with further distinctions present at the level of distal cis elements ( Fig.3D-L , Fig.2H-J ). In fact, at early stages of this convergence process, we find shared trans - and cis -evidence for only three previously identified factors, SOX9, SOX5 and FOXP1 50 ( Fig.4F,G ). This apparent regulatory complexity appears further potentiated by differences in transcription factor binding motifs, in mesenchymal cells across the embryonic origins ( Fig.4H-J ). Whether due to lineage-specific differences in local chromatin environment or the presence of distinct heteromeric binding partners, such expanded ‘regulatory grammar’ at the level of motif diversity is known to increase the combinatorial flexibility and complexity in cell fate decision and patterning processes 30 - 32 . Origin-specific evolutionary trajectories of skeletogenic cells The distinct regulatory modalities that specify the skeletogenic cells in different parts of the vertebrate embryo also hold implications for how we treat them in an evolutionary comparative context. Indeed, thanks to ever more detailed molecular descriptions of cell types across the tree of life, the emerging ‘regulatory blueprints’ have helped to delineate what actually defines distinct cell types, and disentangle putative evolutionary histories and developmental relationships amongst them 1 , 2 . In vertebrates, different anatomical regions have acquired the potential to form an endoskeleton at distinct evolutionary timepoints. It is generally believed that an acellular, cartilage-like support of feeding structures preceded vertebrates 13 . A possible incorporation of the underlying gene regulatory network into skeletogenic precursor cells would then have paved the way for the emergence of the vertebrate cranium 51 - 53 . As additional developmental lineages gained skeletogenic competency, repeated co-option of such network would appear as the most parsimonious solution to developmentally specify a shared cellular phenotype from distinct embryonic sources. However, the partial regulatory independence we document here suggests network recruitment with modification, and thus only an incomplete ‘deep homology’ amongst early skeletogenic cells of different embryonic origins 17 . While these different regulatory strategies might originally have been a necessity – namely, to absorb and integrate the distinct trans - and cis -regulatory modalities of the evolutionarily novel skeletogenic precursor lineages – ultimately, they might have proven beneficial and even adaptive. By relying on partially distinct specification networks, independent changes in effector gene expression – e.g., to affect extracellular matrix composition and ensuing tissue properties 42 , 54 – or cellular growth dynamics become possible in the different parts of the skeleton 55 - 57 . Indeed, naturally occurring genetic variation as well as induced targeted mutations in some of the embryonic origin-specific regulators we document here show anatomical location-restricted effects 57 - 59 . Hence, despite seemingly similar cellular phenotypes, the distinct regulatory strategies at work in the three embryonic lineages may help to make different parts of the vertebrate skeleton become independent targets of evolutionary selection, with distinct biomaterial and patterning properties resulting therefrom. Materials and Methods Tissue Collection Embryos were dissected in ice-cold PBS and staged according to Hamburger and Hamilton 60 (HH). Embryonic tissue was either processed for single-cell functional genomics experiments, or fixed in 4% PFA at 4C. RNA In Situ Hybridization Cranial or brachial tissue samples ranging from HH17 to HH24 were dehydrated and cryo-embedded side-by-side in OCT, to allow for staining of entire time series on single slides. Sectioning was performed on a Leica CM3050S cryostat and RNA in situ hybridization against SOX9 and ACAN was performed using standard protocols 61 . Brightfield images were acquired on an Olympus FLUOVIEW FV3000 and globally processed for color balance and brightness using Adobe Photoshop . Single-Cell RNA-Sequencing (scRNA-seq) Data Collection We sampled the frontonasal prominence at stages HH15, HH18 and HH22, the dorsal part of the brachial region at stages HH12, HH15 and HH20, and entire forelimbs at stages HH21, HH24 and HH27. Tissue dissociation, 10x Genomics Chromium 3’ Kit library preparation and sequencing was performed as reported previously 21 . Data was processed and mapped with CellRanger ( 10x Genomics ), using our in-house improved GRCg6a genome annotation with elongated 3’ UTRs 62 . scRNA-seq Data Processing Unique Molecular Identifier (UMI) count matrices were filtered for quality based on a cell’s total and relative UMI counts (i.e., >4*mean and median + 3*MAD (median absolute deviation) & >0.1, except if UMI count >median), and relation of UMI count/genes detected (i.e., <0.15 & UMI count <2/3) 62 . In total, 8468 cells remained for frontonasal (2702, 4558 and 1208 cells for HH15, HH18 and HH22), 7777 cells for somite (3093, 2993 and 1691 cells for HH12, HH15 and HH20), and 10350 cells for forelimb (2987, 5293 and 2070 cells for HH21, HH24 and HH27). scRNA-seq Data Normalization, Dimensionality Reduction and Clustering Using the R package Seurat (v4) 19 , UMI counts were normalized by sequencing depth and log transformed. A cell cycle score was calculated using SCRAN 63 . Variations of sequencing depth, mitochondrial UMI percentage, and the difference in S and G2M cycle scores were regressed out, using SCTransform from Seurat 19 , 62 , 64 . Genes with a higher value of standardized variance than the sum of median and median absolute deviation were considered as ‘highly variable’. These steps were carried out independently for the three stages of the three embryonic origins. Using Seurat , we integrated samples from the same embryonic origin, and used principal component analysis (PCA) on highly variable genes, followed by tSNE and FFT-accelerated Interpolation-based tSNE algorithms 65 for nonlinear dimensionality reduction on the first 19 (nasal), 21 (somite) and 19 (limb) principal components. Using Seurat functions, we performed Leiden graph-based clustering 66 on all cells with a resolution of 0.2 (=‘broad clustering’). A second round of clustering was conducted on select mesenchyme populations, with resolutions of 0.4 (somite, limb) and 0.5 (nasal) (=‘fine clustering’). Cell type assignments of clusters were based on visual inspection of known marker gene expression patterns, and the activity of the two previously identified early chondrogenic gene expression modules ‘IMM’ 20 and ‘RED’ 21 using the Seurat function ‘AddModuleScore’. Differential Expression Analysis Differential expression analysis was based on a logistic regression framework 67 using Seurat , with cell cycle differences and embryonic stages as latent variables. Genes expressed in at least 10% of the cells and showing differences with an adjusted p-value 0.5 (‘broad’) or >0.25 (‘fine’) were considered as significantly differentially expressed. To minimize batch effects, differential expression analysis of chondrocytes from different embryonic origins was performed on pseudobulk counts using the R package muscat 68 . scRNA-seq Data Integration Across Embryonic Origins We filtered out potential doublets using the R package doubletFinder 69 and removed clusters enriched for mitochondrial counts. The resulting UMI count matrix was divided by size factor and logtransformed using SCRAN 63 . The top variable genes (getTopHVGs, SCRAN ) identified in at least two samples were kept for downstream analyses. Using Seurat , we then integrated the count matrices using anchors in canonical correlation analysis (CCA) reduction, to compute batch corrected matrices of the three embryonic origins. To calculate co-embedding projections, the PCA dimension were reduced sample-wide. Anchors for integration were identified using ‘FindIntegrationAnchors’ in reciprocal PCA reductions. We used ‘IntegrateEmbeddings’ to integrate PCA reduction, followed by tSNE calculations (‘RunTSNE’). Correlation analyses were performed on ‘pseudobulk’ average gene expression values ( Seurat function ‘AverageExpression’) in each cluster. scRNA-seq Pseudotime Analyses We generated spliced/unspliced count matrices of our selected mesenchymal populations using velocyto 70 and assessed the directional transcriptional dynamics of highly variable genes with sufficient spliced/unspliced counts in scVelo with the default parameters 22 . We visualized the recovered dynamics on tSNE projections of the three embryonic origins 65 . We then used these tSNE embeddings as input space, and constructed a minimum spanning tree with preset start cluster in the R package slingshot 25 . Alignment of embryonic origin-specific pseudotimes was performed with TrAGEDy 26 using 40 interpolated points along the respective chondrogenic trajectories. Module expression dissimilarities were calculated by Spearman correlation (1-ρ) and optimal alignment was identified by dynamic time warping with default settings. Using the R package tradeSeq 27 , we detected temporally differentially expressed genes along the respective chondrogenic trajectories. Single-Cell ATAC-Sequencing (scATAC-seq) Data Collection Tissue dissociation was performed as previously described 21 . Cell concentration and viability were assessed using the Nexcelom Cellometer K2 and ∼1*10 6 cells were used to perform nuclei isolation following the 10x Genomics protocol. Nuclei suspensions were loaded onto a Next GEM chip H, and transposition, nuclei partitioning and library preparation were performed according to the 10X Genomics ATAC User Guide. scATAC libraries were quantified on an Agilent 2100 Bioanalyzer system (Agilent) and sequenced on a NovaSeq 6000 system (Illumina). We used CellRanger ATAC v1 . 2 . 0 ( 10x Genomics ) for read processing and quantification, and mapped the fragments to the chicken ENSEMBL genome Gallus_gallus-6.0 71 with our in-house improved GRCg6a genome annotation 62 . PEAK_MERGE_DISTANCE was changed to 50, with all other parameters at default settings. scATAC-seq Data Pre-Processing . We removed doublets with ArchR (v1.0.1) 33 and selected high quality cells in Signac (v1.1.1) 72 using the following thresholds: total number of fragments in peaks ranging from 1000 to 100000, fraction of reads in peaks >15%, nucleosome signal 2. Using these criteria, we ended up with 11527 cells for frontonasal (6171 and 5356 cells for HH15 and HH18), 14106 cells for somite (11232 and 2874 cells for HH12 and HH15), and 10982 cells for forelimb (4453 and 6529 cells for HH21 and HH24). scATAC-seq Data Merging, Dimensionality Reduction and Clustering We merged samples from same embryonic origins and summed fragment counts in 5kb genomic tiling windows located on autosomes and chromosome Z (208680 tiles in total). We performed latent semantic indexing (LSI) dimension reduction on a term frequency-inverse document frequency (TF-IDF) normalized matrix with the top 75% of tiles (top 0.1% tiles are removed, putative repetitive elements or alignment errors) using Signac and removed batch effects on LSI components using Harmony 73 . We used tSNE and FFT-accelerated Interpolation-based t-SNE algorithm 65 to carry out non-linear dimensionality reduction with LSI dimensions 2:30, and performed Leiden graph-based clustering in Seurat on all cells with resolutions of 0.4 (nasal, somite) and 0.6 (limb) (=‘broad clustering’), and a second round of clustering on select mesenchyme populations with resolution of 0.4 (=‘fine clustering’). We annotated cell types for both ‘broad’ and ‘fine’ clusters using scATAC-seq gene activity matrices (promoters and gene bodies) and scRNA-seq expression data, with a combination of label transfers in Seurat and non-negative least squares (NNLS) regression on cluster specific genes 74 , as well as manual inspection of peaks at known marker genes. Peak Calling and Differential Accessibility Analysis We identified peaks using MACS2 (version 2.2.7.1) 24 with parameters “—nomodel --shift 100 --extsize 200 --keep-dup all --call-summits” on pseudobulks of each cluster, for each embryonic origin, respectively. Peaks used summits as center and were extended to a width of 501bp. We merged peaks from different clusters of the same embryonic origin and, for overlapping peaks, kept only the most significant one, using adapted code from ArchR . To get a consensus peak set, we merged peaks from three embryonic origins and removed redundant and/or overlapping peaks using the same logic. Differential accessibility analysis was performed in Seurat , using the total number of fragments and embryonic stages as latent variables. Peaks accessed in at least 10% of the cells and showing differences with adjusted p-value less than 0.05 and a log fold change larger than 0.25 were considered as significantly differentially accessed. Peak-centered heatmaps of differential accessible peaks were visualized with deepTools2 75 . scATAC-seq Data Integration Across Embryonic Origins The top variable peaks (getTopHVGs, SCRAN ) identified in at least two samples were kept for downstream analyses. To calculate a co-embedding projection, first we performed reciprocal LSI dimensional reduction to find anchors ( Seurat function ‘FindIntegrationAnchors’) and constructed transformation matrices between each query cell and anchor. We computed the integration matrices based on the original LSI matrix with dimensions from 2 to 30 and the transformation matrix using the Seurat functions ‘IntegrateEmbeddings’ and ‘runTSNE’ on the integrated LSI dimensions. To remove the batch effects among peak matrices after merging, we binarized the matrix based on presence/absence of counts. Correlation analyses were performed on ‘pseudobulk’ average count values in each cluster using the function ‘AverageExpression’. scATAC-seq Pseudotime Analyses We transferred pseudotime values from our scRNA analyses using the Seurat function ‘TransferData’. First, we integrated scRNA-seq expression matrices and scATAC-seq gene activity matrices and performed CCA dimensional reduction to find anchors. We then constructed a transformation matrix between each query cell and each ancho ( Seurat function ‘FindTransferAnchors’) and computed the transferred scATAC-seq pseudotimes based on the original scRNA-seq pseudotimes and the transformation matrices. For all three embryonic origins, we restricted this transfer to only chondrogenesis related cell type clusters. de novo Motif Enrichment Analysis and Annotation We performed de novo motif enrichment analysis for each cluster, using Homer 28 ‘findMotifsGenome.pl’ with -mset vertebrates -size -250,250 -fdr 5 and motif length between 8-22 bp ( Homer p-value < 1e-11), using highly accessible peaks for each cluster. We obtained candidate TF annotations for this set of de novo motifs with the help of three databases ( Homer vertebrates, JASPAR20 vertebrates 76 , CisBP v2 chicken 77 ) using Homer and STAMP 78 . We selected the best matches based on scRNA-seq expression levels of the predicted TFs, and Spearman correlations between motif activity and gene expression of candidate TFs in scATAC and scRNA aggregates (default k=50, n=400; for small clusters, k=20, n=200). For motifs with similarity scores >0.8, only the one with the lowest p-value was retained. Additionally, we checked for paralog TFs expression along our pseudotime trajectories and calculated its expression correlation to motif activity. We combined annotated de novo motifs from each embryonic origin and calculated motif similarity scores using PWMEnrich 79 . For our final set of annotated de novo motifs, we computed a per-cell motif deviation scores using chromVAR 80 and conducted analysis of differential motif activity using Seurat . Peak-to-Gene Link Analysis We generated imputed pseudoexpression data for each scATAC cell based on scRNA-seq data, using the ArchR function ‘addGeneIntegrationMatrix’. 500 cell aggregates were generated with, with 100 cells per aggregate. We then computed the Pearson correlation between peak accessibility and pseudoexpression in mesenchymal aggregates using the ArchR function ‘addPeak2GeneLinks’. Clustering of peak-to-gene links was calculated by the hkmeans method in factoextra 81 . Functional enrichment analyses of peak-to-gene link clusters was conducted in rGREAT 82 . Enhancer Reporter Assays Genomic regions of candidate enhancers were amplified by PCR from chicken genomic DNA and cloned upstream of a minimal promoter driving GFP expression. Fertilized chicken eggs were incubated at 38.5 °C in a humified incubator to stage HH14 for lateral plate mesoderm and somite electroporation and stage HH7 for neural crest electroporation. DNA solutions containing our enhancer reporter constructs and a constitutively expressing tdTomato co-electroporation control were injected and electroporated into epithelial precursor populations of the cranial neural crest, the brachial somites or the lateral plate mesoderm at forelimb levels 39 , 83 , 84 . Embryos were harvested two days post-electroporation and tissue was processed for immunohistochemistry 85 . Author Contributions This study was conceived and designed by MW, ADPT, CF and PT. Single-cell functional genomics data was generated by CF, ML, CM and AF. Bioinformatics analyses were conducted by MW, CF and PT. Enhancer reporter assays were performed by ADPT and SF. PT wrote the paper, with feedback from all other authors. Competing Interest Statement The authors declare no competing interest. Supplementary Figures Download figure Open in new tab Figure S1. In situ dynamics of skeletogenic induction across three embryonic lineages and quality metrics of single-cell transcriptomics experiments. ( A, B ) Chromogenic in situ hybridizations against SOX9 and ACAN on cranial sagittal ( A ) and brachial transversal ( B ) cryosections spanning embryonic stages HH17 to HH24. First appearance of signal in the respective embryonic lineages is marked by an arrow (NC: neural crest, SOM: somite, LPM: lateral plate mesoderm). ( C - E ) Dot plots of select marker genes expression used for ‘broad’ cluster cell type annotations. ( F ) Violin plots of number of genes detected in the nine scRNA-seq samples, color coded by anatomical origin and embryonic stage. ( G ) Stacked bar plots indicating stage-wise contributions to the total number of cells in ‘broad’ clusters, color coded by anatomical origin and embryonic stage. Download figure Open in new tab Figure S2. The single-cell chromatin accessibility landscape of skeletogenic induction across three embryonic lineages. ( A - C ) tSNE representations of the three scATAC-seq datasets and ‘broad’ cell type annotations, with mesenchymal populations used for ‘fine’ re-clustering highlighted in color. ( D ) Pseudobulk aggregate genome accessibility tracks of peaks near ‘broad’ cell type-specific marker genes. ( E - G ) tSNE representations of re-clustered mesenchymal cells of nasal ( E ), somite ( F ) and limb ( G ) origins, with ‘fine’ cluster annotations and normalized scores based on non-negative least squares (NNLS) regression, between gene activity from ATAC clusters and gene expression from RNA clusters. High values are indicative of cell type correspondences between scATAC and scRNA clusters. Early chondrogenic clusters are highlighted in red. ( H ) Unsupervised hierarchical clustering and heatmap representation of pairwise Spearman’s rank correlation coefficients of differentially accessible peaks in pseudobulk chromatin accessibility data of ‘broad’ and ‘fine’ mesenchymal clusters. Anatomical origins of pseudobulks are indicated by color code. All mesenchymal pseudobulks cluster by embryonic origins, including early chondrogenic cells (highlighted in red) and somite-derived muscles in the limb sample (arrow). Download figure Open in new tab Figure S3. Skeletogenic convergence across three embryonic lineages. ( A - F ) TrAGEDy alignments of embryonic origin-specific chondrogenic pseudotime trajectories. ( A - C ) Pairwise expression dissimilarities of the chondrogenic module ‘IMM’ were calculated by Spearman correlation across all pseudotimes and embryonic origins, and are displayed as 1-ρ. Dynamic time warping identifies skeletogenic convergence across all three embryonic lineages, except for the ‘late chondrocytes’ population in the limb sample that is not accounted for in our nasal and somite data. ( D - F ) TrAGEDy alignment plots of interpolated points of chondrogenesis trajectories and cell type contributions per pseudotime bin, displayed as stacked bar plots. Pseudotime bins where most cells are classified as ‘early chondrogenic’ are highlighted in red. Download figure Open in new tab Figure S4. Core skeletogenic genes and embryonic origin-specific transcriptional dynamics. ( A ) Venn diagram displaying the overlap of differentially expressed genes in chondrogenic cells compared to other mesenchymal populations of the three embryonic origins. ( B - D ) Z-score scaled RNA levels of dynamically expressed genes along the chondrogenic trajectories of nasal ( B ), somite ( C ) and limb ( D ) mesenchymal and chondrogenic cells. SOX9 expression is highlighted by a red arrow. Genes without annotated symbol are labeled with ‘EG_’, followed by the 5 last digits of their ENSGAL identifiers. Download figure Open in new tab Figure S5. De novo identification and annotation of transcription factor binding motifs. ( A ) Transcription factor (TF) family distribution of all de novo identified binding motifs in all populations across the three embryonic origins. ( B ) Publicly available position weight matrices for SOX9 in ‘jaspar20’, ‘cisbp2’ and ‘homer’ databases. Our de novo position weight matrix resembles more closely the one identified using SOX9 ChIP-seq data in chicken embryonic limbs 1 . ( C ) Differential ‘broad’ cluster-specific motif activities across embryonic origins. ( D ) tSNE representations of reclustered mesenchymal cells, with motif activity heatmaps plotted for a core skeletogenic factor (SOX9, middle) and an embryonic origin-specific precursor factor (ZIC2, MEIS1, LH-2A, right). ‘Early chondrogenic’ clusters are highlighted in red, ‘mesenchymal’ clusters in turquoise (left). ( E ) Position weight matrices 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. Acknowledgments The authors wish to thank I. Adameyko and all members of the lab for insightful discussions. Calculations for single-cell functional genomics analyses were performed at sciCORE ( http://scicore.unibas.ch/ ), scientific computing center at the University of Basel. This work was supported by research funds from the Swiss National Science Foundation [SNSF project grant number 310030_189242], the Swiss 3R Competence Centre [3RCC grant OC-2018-005] and the University of Basel to PT. Footnotes 1 SOX9 ChIP-seq position weight matrix from Yamashita et al. 29 , re-used under CC-BY 4.0 license References 1. ↵ Arendt , D. et al. The origin and evolution of cell types . Nat Rev Genet 17 , 744 – 757 ( 2016 ). OpenUrl CrossRef PubMed 2. ↵ Xia , B. & Yanai , I. A periodic table of cell types . Development 146 , dev169854 – 9 ( 2019 ). OpenUrl Abstract / FREE Full Text 3. ↵ Zeng , H. What is a cell type and how to define it? Cell 185 , 2739 – 2755 ( 2022 ). OpenUrl CrossRef PubMed 4. ↵ Peter , I. S. & Davidson , E. H. Assessing regulatory information in developmental gene regulatory networks . Proceedings of the National Academy of Sciences 114 , 5862 – 5869 ( 2017 ). OpenUrl Abstract / FREE Full Text 5. Stadler , T. , Pybus , O. G. & Stumpf , M. P. H. Phylodynamics for cell biologists . Science 371 , ( 2021 ). 6. ↵ Domcke , S. & Shendure , J. A reference cell tree will serve science better than a reference cell atlas . Cell 186 , 1103 – 1114 ( 2023 ). OpenUrl CrossRef PubMed 7. ↵ Wagner , G. P. Homology, genes, and evolutionary innovation . ( 2014 ). 8. ↵ Hall , B. K. Bones and Cartilage . ( Academic Press , 2014 ). 9. ↵ Akiyama , H. et al. Osteo-chondroprogenitor cells are derived from Sox9 expressing precursors . Proc. Natl. Acad. Sci. U.S.A . 102 , 14665 – 14670 ( 2005 ). OpenUrl Abstract / FREE Full Text 10. Eames , B. F. , Sharpe , P. T. & Helms , J. Hierarchy revealed in the specification of three skeletal fates by Sox9 and Runx2 . Dev. Biol . 274 , 188 – 200 ( 2004 ). OpenUrl CrossRef PubMed Web of Science 11. Pacifici , M. et al. Cellular and molecular mechanisms of synovial joint and articular cartilage formation . Ann. N. Y. Acad. Sci . 1068 , 74 – 86 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 12. ↵ Kozhemyakina , E. , Lassar , A. B. & Zelzer , E. A pathway to bone: signaling molecules and transcription factors involved in chondrocyte development and maturation . Development 142 , 817 – 831 ( 2015 ). OpenUrl Abstract / FREE Full Text 13. ↵ Jandzik , D. et al. Evolution of the new vertebrate head by co-option of an ancient chordate skeletal tissue . Nature 518 , 534 – 537 ( 2015 ). OpenUrl CrossRef 14. ↵ Hirasawa , T. & Kuratani , S. Evolution of the vertebrate skeleton: morphology, embryology, and development . Zoological Letters 1 , 2 – 17 ( 2015 ). OpenUrl 15. ↵ Tarazona , O. A. , Slota , L. A. , Lopez , D. H. , Zhang , G. & Cohn , M. J. The genetic program for cartilage development has deep homology within Bilateria . Nature 533 , 86 – 89 ( 2016 ). OpenUrl CrossRef 16. Brunet , T. & Arendt , D. Animal Evolution: The Hard Problem of Cartilage Origins . Current Biology 26 , 1 – 4 ( 2016 ). OpenUrl CrossRef PubMed 17. ↵ Tschopp , P. & Tabin , C. J. Deep homology in the age of next-generation sequencing . Philosophical Transactions of the Royal Society B: Biological Sciences 372 , 20150475 ( 2017 ). OpenUrl CrossRef PubMed 18. ↵ Kiani , C. , Chen , L. , Wu , Y. J. , Yee , A. J. & Yang , B. B. Structure and function of aggrecan . Cell Research 12 , 19 – 32 ( 2002 ). OpenUrl CrossRef PubMed Web of Science 19. ↵ Hao , Y. et al. Integrated analysis of multimodal single-cell data . Cell 1 – 45 ( 2021 ). 20. ↵ Picos , P. G. & Eames , B. F. Limb mesoderm and head ectomesenchyme both express a core transcriptional program during chondrocyte differentiation . Front. Cell Dev. Biol . ( 2022 ). 21. ↵ Feregrino , C. , Sacher , F. , Parnas , O. & Tschopp , P. A single-cell transcriptomic atlas of the developing chicken limb . BMC Genomics 20 , 401 ( 2019 ). OpenUrl CrossRef 22. ↵ Bergen , V. , Lange , M. , Peidli , S. , Wolf , F. & Theis , F. J. Generalizing RNA velocity to transient cell states through dynamical modeling . Nat Biotechnol 1 – 13 ( 2020 ). 23. ↵ Buenrostro , J. D. et al. Single-cell chromatin accessibility reveals principles of regulatory variation . Nature ( 2015 ). 24. ↵ Zhang , Y. et al. Model-based analysis of ChIP-Seq (MACS) . Genome Biol . 9 , R137 ( 2008 ). OpenUrl CrossRef PubMed 25. ↵ Street , K. et al. Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics . BMC Genomics 19 , 477 ( 2018 ). OpenUrl CrossRef PubMed 26. ↵ Laidlaw , R. F. , Briggs , E. M. , Matthews , K. R. , McCulloch , R. & Otto , T. D. TrAGEDy: Trajectory Alignment of Gene Expression Dynamics . bioRxiv 2022.12.21.521424 ( 2022 ). 27. ↵ Van den Berge , K. et al. Trajectory-based differential expression analysis for single-cell sequencing data . Nat Comms 1 – 13 ( 2020 ). 28. ↵ Heinz , S. et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities . Molecular Cell 38 , 576 – 589 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 29. ↵ Yamashita , S. et al. Comparative analysis demonstrates cell type-specific conservation of SOX9 targets between mouse and chicken . Sci. Rep . 9 , 1 – 15 ( 2019 ). OpenUrl CrossRef PubMed 30. ↵ Jolma , A. et al. DNA-dependent formation of transcription factor pairs alters their binding specificity . Nature 527 , 384 – 388 ( 2015 ). OpenUrl CrossRef PubMed 31. Hojo , H. , Ohba , S. , He , X. , Lai , L. P. & McMahon , A. P. Sp7/Osterix Is Restricted to Bone-Forming Vertebrates where It Acts as a Dlx Co-factor in Osteoblast Specification . Developmental Cell 37 , 238 – 253 ( 2016 ). OpenUrl CrossRef 32. ↵ Kim , S. et al. DNA-guided transcription factor cooperativity shapes face and limb mesenchyme . Cell 1 – 47 ( 2024 ). 33. ↵ Granja , J. M. et al. ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis . Nature Genetics 53 , 403 – 411 ( 2021 ). OpenUrl 34. ↵ Breschi , A. et al. A limited set of transcriptional programs define major cell types . Genome Research 30 , 1047 – 1059 ( 2020 ). OpenUrl Abstract / FREE Full Text 35. ↵ Nieto , M. A. Epithelial plasticity: A common theme in embryonic and cancer cells . Science 342 , 1234850 ( 2013 ). OpenUrl Abstract / FREE Full Text 36. ↵ Haerinck , J. , Goossens , S. & Berx , G. The epithelial–mesenchymal plasticity landscape: principles of design and mechanisms of regulation . Nat Rev Genet 24 , 590 – 609 ( 2023 ). OpenUrl 37. ↵ Sauka-Spengler , T. & Bronner-Fraser , M. Development and evolution of the migratory neural crest: a gene regulatory perspective . Current Opinion in Genetics & Development 16 , 360 – 366 ( 2006 ). OpenUrl 38. Scaal , M. Early development of the vertebral column . Seminars in Cell & Developmental Biology 49 , 83 – 91 ( 2016 ). OpenUrl CrossRef PubMed 39. ↵ Gros , J. & Tabin , C. J. Vertebrate limb bud formation is initiated by localized epithelial-to-mesenchymal transition . Science 343 , 1253 – 1256 ( 2014 ). OpenUrl Abstract / FREE Full Text 40. ↵ Ohba , S. , He , X. , Hojo , H. & McMahon , P. Distinct Transcriptional Programs Underlie Sox9 Regulation of the Mammalian Chondrocyte . CellReports 12 , 229 – 243 ( 2015 ). OpenUrl 41. ↵ Darbellay , F. , Ramisch , A. , Lopez-Delisle , L. , Kosicki , M. , Visel , A. & Andrey , G. Chondrogenic Enhancer Landscape of Limb and Axial Skeleton Development . bioRxiv 1 – 27 ( 2023 ). 42. ↵ Sacher , F. , Feregrino , C. , Tschopp , P. & Ewald , C. Y. Extracellular matrix gene expression signatures as cell type and cell state identifiers . Matrix Biology Plus 100069 ( 2021 ). 43. ↵ Konstantinides , N. et al. Phenotypic Convergence: Distinct Transcription Factors Regulate Common Terminal Features . Cell ( 2018 ). 44. ↵ Sáez , M. , Briscoe , J. & Rand , D. A. Dynamical landscapes of cell fate decisions . Interface Focus 12 , 20220002 ( 2022 ). OpenUrl 45. ↵ Attanasio , C. et al. Fine Tuning of Craniofacial Morphology by Distant-Acting Enhancers . Science 342 , 1241006 – 1241006 ( 2013 ). OpenUrl Abstract / FREE Full Text 46. Yao , B. et al. The SOX9 upstream region prone to chromosomal aberrations causing campomelic dysplasia contains multiple cartilage enhancers . Nucleic Acids Research 43 , 5394 – 5408 ( 2015 ). OpenUrl CrossRef PubMed 47. ↵ Chen , L.-F. et al. Structural elements promote architectural stripe formation and facilitate ultra-long-range gene regulation at a human disease locus . Molecular Cell 1 – 32 ( 2023 ). 48. ↵ Jo , A. et al. The versatile functions of Sox9 in development, stem cells, and human diseases . Genes & Diseases 1 , 149 – 161 ( 2014 ). OpenUrl 49. ↵ Akiyama , H. et al. Misexpression of Sox9 in mouse limb bud mesenchyme induces polydactyly and rescues hypodactyly mice . Matrix Biology 26 , 224 – 233 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 50. ↵ Liu , C.-F. , Samsa , W. E. , Zhou , G. & Lefebvre , V. Transcriptional control of chondrocyte specification and differentiation . Seminars in Cell & Developmental Biology 62 , 34 – 49 ( 2017 ). OpenUrl CrossRef 51. ↵ Zhang , G. , Eames , B. F. & Cohn , M. J. Chapter 2. Evolution of vertebrate cartilage development . Curr. Top. Dev. Biol . 86 , 15 – 42 ( 2009 ). OpenUrl CrossRef PubMed Web of Science 52. Kaucka , M. & Adameyko , I. Evolution and development of the cartilaginous skull: From a lancelet towards a human face . Seminars in Cell & Developmental Biology 91 , 2 – 12 ( 2019 ). OpenUrl 53. ↵ Saunders , L. M. et al. Embryo-scale reverse genetics at single-cell resolution . Nature 1 – 34 ( 2023 ). 54. ↵ Humphrey , J. D. , Dufresne , E. R. & Schwartz , M. A. Mechanotransduction and extracellular matrix homeostasis . Nat Rev Mol Cell Biol 1 – 11 ( 2014 ). 55. ↵ Cretekos , C. J. et al. Regulatory divergence modifies limb length between mammals . Genes & Development 22 , 141 – 151 ( 2008 ). OpenUrl Abstract / FREE Full Text 56. Cooper , K. L. et al. Multiple phases of chondrocyte enlargement underlie differences in skeletal proportions . Nature 495 , 375 – 378 ( 2013 ). OpenUrl CrossRef PubMed 57. ↵ Lamichhaney , S. et al. Evolution of Darwin’s finches and their beaks revealed by genome sequencing . Nature 1 – 16 ( 2015 ). 58. Pineault , K. M. , Song , J. Y. , Kozloff , K. M. , Lucas , D. & Wellik , D. M. Hox11 expressing regional skeletal stem cells are progenitors for osteoblasts, chondrocytes and adipocytes throughout life . Nat Comms 1 – 15 ( 2019 ). 59. ↵ Ushiki , A. et al. Deletion of Pax1 scoliosis-associated regulatory elements leads to a female-biased tail abnormality . CellReports 43 , 113907 ( 2024 ). OpenUrl 60. ↵ Hamburger , V. & Hamilton , H. L. A series of normal stages in the development of the chick embryo . J. Morphol . 88 , 49 – 92 ( 1951 ). OpenUrl CrossRef PubMed Web of Science 61. ↵ McGlinn , E. & Mansfield , J. H. Detection of gene expression in mouse embryos and tissue sections . Methods Mol. Biol . 770 , 259 – 292 ( 2011 ). OpenUrl CrossRef PubMed 62. ↵ Feregrino , C. & Tschopp , P. Assessing evolutionary and developmental transcriptome dynamics in homologous cell types . Dev. Dyn . 251 , 1472 – 1489 ( 2022 ). OpenUrl 63. ↵ Lun , A. T. L. , McCarthy , D. J. & Marioni , J. C. A step-by-step workflow for low-level analysis of single-cell RNA-seq data . F1000Res 5 , 2122 ( 2016 ). OpenUrl 64. ↵ Scialdone , A. et al. Computational assignment of cell-cycle stage from single-cell transcriptome data . Methods 85 , 54 – 61 ( 2015 ). OpenUrl CrossRef PubMed 65. ↵ Linderman , G. C. , Rachh , M. , Hoskins , J. G. , Steinerberger , S. & Kluger , Y. Fast interpolation-based t-SNE for improved visualization of single-cell RNA-seq data . Nat Meth 1 – 10 ( 2019 ). 66. ↵ Traag , V. A. , Waltman , L. & van Eck , N. J. From Louvain to Leiden: guaranteeing well-connected communities . Sci. Rep . 9 , 5233 – 12 ( 2019 ). OpenUrl CrossRef PubMed 67. ↵ Ntranos , V. , Yi , L. , Melsted , P. & Pachter , L. A discriminative learning approach to differential expression analysis for single-cell RNA-seq . Nat Meth 1 – 8 ( 2019 ). 68. ↵ Crowell , H. L. et al. muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data . Nat Comms 1 – 12 ( 2020 ). 69. ↵ McGinnis , C. S. , Murrow , L. M. & Gartner , Z. J. DoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors . Cell Systems 8 , 329 – 337 .e4 ( 2019 ). OpenUrl 70. ↵ La Manno , G. et al. RNA velocity of single cells . Nature 1 – 25 ( 2018 ). 71. ↵ Yates , A. D. et al. Ensembl 2020 . Nucleic Acids Research 48 , D682 – D688 ( 2020 ). OpenUrl CrossRef PubMed 72. ↵ Stuart , T. , Srivastava , A. , Madad , S. , Lareau , C. A. & Satija , R. Single-cell chromatin state analysis with Signac . Nat Meth 18 , 1333 – 1341 ( 2021 ). OpenUrl 73. ↵ Korsunsky , I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony . Nat Meth 1 – 16 ( 2020 ). 74. ↵ Domcke , S. et al. A human cell atlas of fetal chromatin accessibility . Science 370 , eaba7612 – 17 ( 2020 ). OpenUrl Abstract / FREE Full Text 75. ↵ Ramírez , F. et al. deepTools2: a next generation web server for deep-sequencing data analysis . Nucleic Acids Research 44 , W160 – 5 ( 2016 ). OpenUrl CrossRef PubMed 76. ↵ Fornes , O. et al. JASPAR 2020: Update of the open-Access database of transcription factor binding profiles . Nucleic Acids Research 48 , D87 – D92 ( 2020 ). OpenUrl CrossRef PubMed 77. ↵ Weirauch , M. T. et al. Determination and inference of eukaryotic transcription factor sequence specificity . Cell 158 , 1431 – 1443 ( 2014 ). OpenUrl CrossRef PubMed Web of Science 78. ↵ Mahony , S. & Benos , P. V. STAMP: A web tool for exploring DNA-binding motif similarities . Nucleic Acids Research 35 , W253 – 8 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 79. ↵ Stojnić , R. & Diez , D. Package ‘PWMEnrich’ . 1 – 75 ( 2015 ). 80. ↵ Schep , A. N. , Wu , B. , Buenrostro , J. D. & Greenleaf , W. J. chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data . Nat Meth 14 , 975 – 978 ( 2017 ). OpenUrl 81. ↵ Kassambara , A. & Mundt , F. factoextra : Extract and Visualize the Results of Multivariate Data Analyses . 1 – 2 ( 2024 ). 82. ↵ Gu , Z. & Hübschmann , D. rGREAT: an R/bioconductor package for functional enrichment on genomic regions . Bioinformatics 39 , ( 2023 ). 83. ↵ Sauka-Spengler , T. & Barembaum , M. Gain- and Loss-of-Function Approaches in the Chick Embryo . 1 – 20 ( 2008 ). 84. ↵ Scaal , M. , Gros , J. , Lesbros , C. & Marcelle , C. In ovo electroporation of avian somites . Dev. Dyn . 229 , 643 – 650 ( 2004 ). OpenUrl CrossRef PubMed Web of Science 85. ↵ Tschopp , P. et al. A relative shift in cloacal location repositions external genitalia in amniote evolution . Nature 516 , 391 – 394 ( 2014 ). OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted March 29, 2024. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Distinct Gene Regulatory Dynamics Drive Skeletogenic Cell Fate Convergence During Vertebrate Embryogenesis Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Distinct Gene Regulatory Dynamics Drive Skeletogenic Cell Fate Convergence During Vertebrate Embryogenesis Menghan Wang , Ana Di Pietro-Torres , Christian Feregrino , Maëva Luxey , Chloé Moreau , Sabrina Fischer , Antoine Fages , Patrick Tschopp bioRxiv 2024.03.26.586769; doi: https://doi.org/10.1101/2024.03.26.586769 Share This Article: Copy Citation Tools Distinct Gene Regulatory Dynamics Drive Skeletogenic Cell Fate Convergence During Vertebrate Embryogenesis Menghan Wang , Ana Di Pietro-Torres , Christian Feregrino , Maëva Luxey , Chloé Moreau , Sabrina Fischer , Antoine Fages , Patrick Tschopp bioRxiv 2024.03.26.586769; doi: https://doi.org/10.1101/2024.03.26.586769 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Developmental Biology Subject Areas All Articles Animal Behavior and Cognition (7643) Biochemistry (17717) Bioengineering (13910) Bioinformatics (42018) Biophysics (21480) Cancer Biology (18629) Cell Biology (25537) Clinical Trials (138) Developmental Biology (13392) Ecology (19935) Epidemiology (2067) Evolutionary Biology (24356) Genetics (15617) Genomics (22531) Immunology (17755) Microbiology (40438) Molecular Biology (17200) Neuroscience (88706) Paleontology (667) Pathology (2840) Pharmacology and Toxicology (4832) Physiology (7657) Plant Biology (15171) Scientific Communication and Education (2046) Synthetic Biology (4304) Systems Biology (9828) Zoology (2272)