A single-cell transcriptomic atlas of the developing chicken limb

Background Through precise implementation of distinct cell type specification programs, differentially regulated in both space and time, complex patterns emerge during organogenesis. Thanks to its easy experimental accessibility, the developing chicken limb has long served as a paradigm to study vertebrate pattern formation. Through decades’ worth of research, we now have a firm grasp on the molecular mechanisms driving limb formation at the tissue-level. However, to elucidate the dynamic interplay between transcriptional cell type specification programs and pattern formation at its relevant cellular scale, we lack appropriately resolved molecular data at the genome-wide level. Here, making use of droplet-based single-cell RNA-sequencing, we catalogue the developmental emergence of distinct tissue types and their transcriptome dynamics in the distal chicken limb, the so-called autopod, at cellular resolution. Results Using single-cell RNA-sequencing technology, we sequenced a total of 17,628 cells coming from three key developmental stages of chicken autopod patterning. Overall, we identified 23 cell populations with distinct transcriptional profiles. Amongst them were small, albeit essential populations like the apical ectodermal ridge, demonstrating the ability to detect even rare cell types. Moreover, we uncovered the existence of molecularly distinct sub-populations within previously defined compartments of the developing limb, some of which have important signaling functions during autopod pattern formation. Finally, we inferred gene co-expression modules that coincide with distinct tissue types across developmental time, and used them to track patterning-relevant cell populations of the forming digits. Conclusions We provide a comprehensive functional genomics resource to study the molecular effectors of chicken limb patterning at cellular resolution. Our single-cell transcriptomic atlas captures all major cell populations of the developing autopod, and highlights the transcriptional complexity in many of its components. Finally, integrating our data-set with other single-cell transcriptomics resources will enable researchers to assess molecular similarities in orthologous cell types across the major tetrapod clades, and provide an extensive candidate gene list to functionally test cell-type-specific drivers of limb morphological diversification. Electronic supplementary material The online version of this article (10.1186/s12864-019-5802-2) contains supplementary material, which is available to authorized users.


3
Background 50 Embryonic pattern formation relies on the tight coordination of numerous developmental 51 processes, across multiple scales of complexity. From seemingly homogenous progenitor 52 populations, different cell types get specified and arranged in intricate patterns, to give rise to 53 functional tissues and organs. As progenitors mostly share a common genome, this 54 phenotypic specialization relies on the precise execution of distinct gene regulatory networks, 55 to enable cell type specification and ensuing pattern formation [1][2][3]. Slight deviations in 56 these processes contribute to morphological variations within natural populations. More 57 profound aberrations, however, can cause malformations and ultimately result in death of the 58 embryo. To buffer such fragile balance, many cell type specification and pattering processes 59 rely on complex feedback mechanisms, through tightly interconnected molecular loops 60 between spatially distinct signaling centers [4-6] Hence, integration of multiple signaling 61 pathways across space and time define a molecular coordinate grid to instruct organogenesis 62 at the tissue level. Ultimately, however, these multifaceted signaling inputs have to be 63 incorporated at the cellular level, via cell type-specifying gene regulatory networks, as 64 progenitor cells undergo spatially and temporally defined cell fate decisions to contribute to 65 proper pattern formation. 66 Tetrapod limb development has long served as a model to study the genetic and molecular 67 underpinnings of vertebrate pattern formation. Due to its non-essentiality for embryo survival, 68 many fetuses carrying mutations that affect limb development make it to full term. 69 Accordingly, human geneticists have been able to accumulate an impressive catalogue of 70 candidate genes for limb patterning [7][8][9]. Combined with the easy accessibility of the limb in 71 chicken embryos, and molecular genetic tools in the mouse, decades of experimental work 72 have resulted in an in-depth understanding of many of the molecular mechanisms driving limb   Digit elongation then relies on a specialized distal progenitor population, which supports 91 outgrowth of individual digit bones, the phalanges [22,23]. Digit-specific phalanx-formulas, 92 and their stereotypic connection patterns via synovial joints, are established by signals 93 emanating from the posterior interdigit mesenchyme [24,25]. 94 In this study, capitalizing on the power of droplet-based single-cell RNA-sequencing, we 95 resolve the underlying transcriptional dynamics of autopod tissue formation and pattern 96 emergence at single-cell resolution, across three stages of chicken hindlimb development. In 97 total, we present transcriptomic data for 17,628 cells, allowing us to identify all major tissue 98 types of the developing limb, as well as a substantial amount of molecular heterogeneity 99 therein. Through weighted correlation network analysis, we define distinct gene co-expression 100 modules that track corresponding tissue types across developmental time. Finally, we focus 101 on the molecular make-up of cell populations involved in digit pattern formation and, hence, 102 putative drivers of morphological diversification in the autopod.

103
Collectively, we present a comprehensive genomics resource that for the first time reveals the 104 transcriptome dynamics of the developing chicken foot at cellular level. Our study identifies 105 novel and known marker genes in co-expression modules of patterning-relevant cell 106 populations, thereby providing an extensive catalogue of candidate genes for functional 107 follow-up studies, to elucidate the molecular mechanisms of autopod pattern formation and 108 diversification. particularly for the skeletal apparatus and its associated tissues. Namely, stage HH25 is 117 dominated by overall autopod outgrowth and delineation of the main embryonic axes, at 118 HH29 digit-specific patterns differentiate, and at HH31 digit elongation is phasing out. We 119 designed our tissue sampling strategies accordingly. At HH25, we captured the entire distal 120 part of the growing limb (Fig. 1a), at HH29 we dissected two digits with distinct skeletal 121 formulas, digit 3 and 4, as well as their adjacent interdigit mesenchyme (Fig. 1b), and at 122 HH31 we focused on the tip of digit 4 with its growth-relevant progenitor population ( Fig.   123 1c). We dissociated the micro-dissected tissue pieces using enzymatic digest combined with 124 mechanical shearing and prepared single-cell suspensions for droplet-based high-throughput  On average, we detected 2,879 unique molecular identifiers (UMIs) and 1,081 genes per cell 135 (Additional file 1: Fig. S1b,c).

136
Autopod tissue composition at cellular resolution 137 Using unsupervised graph-based clustering, we identified 5, 10 and 5 clusters at stages HH25, 138 HH29 and HH31, respectively. Projecting these clusters onto stage-specific tSNE (t- well as distinct outlier groups ( Fig. 1 d-f). Based on the expression of known marker genes 142 and gene ontology (GO)-term enrichment analyses, we were able to attribute these broadly S2a-c). At stage HH25, they comprise a largely undifferentiated and proliferating 145 mesenchymal population (red), early skeletal progenitors (blue), muscle cells invading the 146 limb (black), as well as skin (purple) and blood cells (grey) (Fig. 1d,g). We recovered cell 147 populations corresponding to those same five tissue types in our HH29 sample, with the 148 exception that the "blood cluster" was now dominated by white blood cells and not  Although all expected major tissue types were recovered in our primary analyses, smaller cell 165 populations, some well known to be essential for limb outgrowth and patterning, remained 166 elusive. Hence, given our sampling depth, we next examined our data for additional sub-167 structure. Indeed, upon closer inspection using finer-tuned clustering parameters, we did find To gain further insights into the regulatory programs that maintain these transcriptional 189 signatures, and explore their potential biological significance, we tested for the occurrence of  Transcriptionally and spatially distinct sub-populations in the interdigit mesenchyme 230 As expected by developmental stage, interdigit populations were only recovered in samples 231 HH29 and HH31. In total, we identified four associated co-expression modules (Fig. 4a-d).

232
High Orange and Olivegreen module activities were coinciding with the same interdigit sub-233 population (Fig. 4e,f), which was recognizable in both HH29 and HH31 samples and marked 234 by RDH10 expression (Fig. 2e,f). Noticeably, all genes with high membership in module 235 Olivegreen were transcription factors (TFs), while module Orange was enriched for 236 enzymatic activities (Fig. 4a,b). Both, however, scored high for GO-terms related to retinoic  (Fig. 4c,g). Lastly, module Midnightblue showed multiple TFs and its activity was 241 restricted to HH29 sub-cluster 2 (Fig. 4d,h).

242
Since relevant patterning information is contained in the interdigit, posteriorly adjacent to 243 each forming digit, we next wondered whether some of the sub-clustering structure 244 corresponded to spatially distinct interdigit populations along the anterior-posterior axis of the 245 9 autopod. At HH29, we detected three interdigit sub-clusters (Fig. 4i). Using differential 246 expression analyses, we defined marker genes that distinguish the three sub-clusters from 247 each other (Fig. 4j). To assign putative spatial information to our single-cell interdigit 248 transcriptomes, we reanalyzed a bulk RNA-seq dataset covering stages HH29 and HH31 of 249 the developing chicken hindlimb autopod [34]. This dataset is based on dissections of 250 individual digits, together with their posteriorly associated interdigit mesenchyme, and thus 251 provided an opportunity to identify spatially resolved marker genes. We contrasted their 252 transcriptomic data of digit/interdigit III against digit/interdigit IV and found a total of 54 253 genes to be significantly differentially expressed at both developmental time points (Fig. 4k). 254 Comparing the digit/interdigit IV-specific subset of these genes to our differential expression   (Fig. 5i,j). Finally, we identified skeletal progenitor populations at all three time points (Fig. 6a-c). 280 According to the developmental stages we sampled, only cartilage-producing skeletal cells 281 were recovered. In all three samples, we found a cell population resembling early 282 chondrocytes (sub-clusters HH25-4, HH29-15 and HH31-2). At stages HH29 and HH31, a 283 seemingly more mature chondrocyte type emerged (HH29-3, HH31-1), and an additional 284 cartilaginous cluster was evident in the HH29 sample (HH29-17). Concomitantly, we 285 identified two co-expression modules associated with these cell populations, Turquoise and 286 Red (Fig 6d,e). Turquoise is centered on CD24, CHGB and SULF1, whereas module Red  (Fig. 6f). Interestingly, compared to module Turquoise, the activity of 292 module Red was generally more restricted and specifically excluded from sub-cluster HH29-293 17 (Fig. 6g,h). Upon closer inspection, we identified high expression of several known 294 synovial joint markers genes in this population, thus identifying it as the forming 295 interphalangeal joints (Fig. 6i, Additional file 3). 296 Hence, through a combination of differential gene expression and GO-term enrichment 297 analyses, as well as gene co-expression modules, we identified spatially and/or temporally   to form (Fig. 2d, Fig. 3b). As such, it suggests an early lineage priming, without necessarily 333 implying a definite switch in cell fate or clear morphological distinctions. In agreement with 334 12 this, our ZFHX3-containing module Darkgreen appears to be the most basic and least specific 335 of the co-expression modules that coincide with the nsCT population. We detect its activity at 336 all three time points, marking the prospective nsCT as well as parts of the PRRX1-positive 337 mesenchymal progenitor population (Fig. 5c,f). Only later do more mature and restricted 338 nsCT sub-divisions and their corresponding co-expression modules occur, as exemplified by 339 the activity of module Darkgrey and some of its members known to be involved in the 340 formation of periskeletal tissues and tendon attachment sites (Fig. 5a,d)  module Magenta with TFAP2B, WNT5A and high BMP signaling (Fig. 3c-f). Certain module 373 members have been functionally implied in regulating autopod growth and digit elongation 374 [24,[54][55][56]], yet others remain completely unexplored in this context. 375 Moreover, we identify distinct sub-populations of interdigit mesenchyme cells in our HH29 376 and HH31 samples, with four associated gene co-expression modules (Fig. 4a-h). Module 377 Olivegreen contains SNAI and ID genes, known to be expressed in interdigits, and likely 378 relates to the various BMP-driven processes in this tissue [57][58][59][60][61][62]. On the other hand, module

389
On the other hand, genes in module Turquoise do not, for the most part, evoke a classical 390 chondrogenic transcriptional profile (Fig. 6d). Again, this module might rather reflect an early 391 transcriptional priming, only this time towards the skeletogenic lineage. In agreement with 392 this, we only detect low expression levels for the canonical early skeletogenic marker SOX9 in 393 HH25 sub-cluster 4 (Fig. 2d), which itself is specifically enriched for Turquoise activity.

394
Likewise, our synovial joint-like HH29 sub-cluster 17 shows high activity for Turquoise, 395 while excluding the more mature chondrocyte module Red (Fig. 6g-i).

411
Tissue sampling 412 We collected tissue samples from embryonic hind limbs at different developmental stages 413 (Fig. 1,a-c). Limbs were dissected in cold PBS, and chopped coarsely with a razorblade.  Filtering thresholds for mapped data were adapted for each sample, depending on the different 433 library complexities. Cells with an UMI count of more than 4 times the sample mean or less 434 than 20% of the sample median were filtered out, cells with a mitochondrial or ribosomal 435 contribution to UMI count of more than 10% were also filtered out. Using the R package