Single cell RNA-sequencing and RNA-tomography of the avian embryo extending body axis

Introduction

The generation of somites, which arise in a regular sequence during embryogenesis, is fundamental for creating the vertebrate segmented body plan (Benazeraf and Pourquie, 2013). Pairs of somites form on either side of the neural tube from unsegmented, paraxial mesoderm, and the process of somitogenesis, which involves waves of cycling gene expression, has been studied extensively in chick embryos (Pourquie, 2004). Prospective paraxial mesoderm cells emerge from the primitive streak during gastrulation (Psychoyos and Stern, 1996) and follow a stereotypical migration trajectory towards their destination (Yang et al., 2002; Iimura et al., 2007). As the body axis elongates, bi-potential neuromesodermal progenitors (NMP) located in the tailbud continue to generate paraxial mesoderm and cells of the neural tube (Wilson et al., 2009; Henrique et al., 2015; Wymeersch et al., 2021). The dynamics of this specialised cell population has been mapped in detail in chick embryos (Guillot et al., 2021) and it has been shown that the extension of neural and paraxial mesoderm tissues in the embryonic body is coordinated by mechanical interactions (Xiong et al., 2020).

Somite differentiation proceeds along the posterior-to-anterior axis and serves as a paradigm for the study of cell fate specification. Multiple signals from surrounding tissues are integrated by somite cells to produce the lineages of the musculoskeletal system, including chondrocytes of the axial skeleton and skeletal muscles of the trunk and limbs (Brent and Tabin, 2002; Christ et al., 2007; Christ and Scaal, 2008). Cell fate specification is intimately linked to stereotypic morphological changes resulting in somite compartmentalisation. For example, tracking of GFP-labelled cells showed that the dorsal dermomyotome produces the myotome layer in multiple waves, with the first myocytes specified adjacent to the neural tube (Gros et al., 2004). Live-imaging of cellular rearrangements examined the morphological transformations of somites, from epithelial structures to somites with a mesenchymal sclerotome, located ventrally, and an epaxial myotome abutting the neural tube (McColl et al., 2018). This uncovered differential cell sizes and regions of proliferation as well as a directed movement of dermomyotomal progenitor cells towards the rostro-medial domain of the dermomyotome, where skeletal muscle formation initiates.

To better characterise the regulation of these morphogenetic events and their integration with cell specification and differentiation, we previously used bulk sequencing to assess the dynamic changes of the transcriptome and of chromatin accessibility across presegmented mesoderm and early, maturing and differentiating somites (Mok et al., 2021). Associating differentially accessible chromatin with nearby genes, differentially expressed along the axis, identified candidate cis-regulatory elements (CREs) involved in expression of transcription factors important for somite formation and differentiation. Time-lapse microscopy in accessible chick embryos of fluorescent CRE-reporters revealed their spatio-temporal activity, and mutation analysis uncovered some upstream regulators. Similarly in mice, we examined matched gene expression and open chromatin profiles for newly formed somite pairs across a developmental time series. This provided a high-resolution view of the molecular signatures underlying the conserved maturation programme followed by all somites after segmentation (Ibarra-Soria et al., 2023).

Here we focus on the initial phase of trunk development in chick embryos. This complements recent in vitro organoid models of axis elongation, which are based on the differentiation of mouse or human pluripotent stem cells and use pharmacological activation or inhibition of crucial signalling pathways (Veenvliet and Herrmann, 2021). These protocols generate self-organising structures including gastruloids (Turner et al., 2017; Moris et al., 2020), trunk-like structures (Veenvliet et al., 2020), somitoids (Sanaki-Matsumiya et al., 2022), or axioloids (Yamanaka et al., 2023). These structures often comprise mesoderm, including somites, although the notochord, which is involved in patterning of trunk tissues, is missing. Thus, it is important to reconstruct the molecular profiles and cellular composition in the native tissues.

Here, we use single cell transcriptomics combined with an RNA-tomography based approach, analogous to Tomo-seq (Junker et al., 2014; Kruse et al., 2016) but using a modified G&T-seq approach (Macaulay et al., 2015). This generated a spatio-temporal map of the emerging trunk and identified genes not previously known to be involved in pre-somitic/pre-segmented mesoderm (psm) and somite maturation. Our study complements data in chick embryos, from earlier developmental stages (HH4-HH11) (Vermillion et al., 2018; Williams et al., 2022; Rito et al., 2023), from tailbud (Guillot et al., 2021) and from prospective neural plate, neural plate border and non-neural ectoderm (Trevers et al., 2023). The dataset is relevant for cell type specification during early body formation and may provide insights into the molecular genetics that underlie diseases of the musculoskeletal system.

Results

Single cell profiling of the developing chick embryonic body

To investigate the molecular basis underpinning somitogenesis and axis elongation of the growing chick body, we mapped the transcriptomes of individual cells at embryonic stage HH14 (Hamburger and Hamilton, 1992) (Figure 1A). This stage embryo has 22 somites, including cervical level somites (6–19) and thoracic level somites (20–22). The unsegmented paraxial mesoderm and tailbud comprise prospective somites of the thoracic, lumbar and sacral regions (Weldon and Munsterberg, 2022). Single cell suspensions from the posterior part of five pooled embryos included the extraembryonic region, tailbud, pre-somitic mesoderm and the most recently formed six somites (Figure 1B). Following enzymatic digestion and mechanical dissociation the suspension was processed using the 10X Genomics Chromium. A total of 6158 cells were sequenced with a median of 517 genes and 900 UMIs per cell.

Figure 1. Cell population composition and signatures of the HH14 chicken embryo trunk. (A) Trunk regions, indicated by stippled lines, of five stage HH14 chicken embryos were collected for scRNA-seq using 10X Genomics Chromium. (B) Illustration of developing tissues captured: pre-somitic mesoderm (psm), epithelial somites, maturing somites, lateral plate mesoderm (lpm), surface ectoderm, neural tube and neural crest cells, and notochord. (C) Unsupervised UMAP subdivides cells within the trunk into 10 clusters—lateral plate mesoderm, neural progenitors, early somite, maturing somite, pre-somitic mesoderm, ectoderm, blood progenitors, endothelial progenitors and neural crest. (D) Heatmap of the top 10 genes significantly enriched in each cluster; representative genes are shown. (E) UMAPs show log normalised counts of a representative gene for each cluster. Colour intensity is proportional to expression level of each gene. (F) Distribution of cell cycle phases visualised using Seurat cell cycle scoring.

Unsupervised clustering was used to classify cell populations (Butler et al., 2018). Projection onto UMAP plots revealed 10 separate clusters (Figure 1C). These clusters were characterised by expression of classic marker genes, which assigned cluster identity and tissues. These included for the lateral plate mesoderm (Prrx1 and Krt18) (Bell et al., 2004), for neural cells (Wnt4 and Olig2) (Cauthen et al., 2001; Zhou et al., 2001), for epithelial somites (Meox1 and Tcf15) (Stockdale et al., 2000; Reijntjes et al., 2007; Berti et al., 2015), for maturing somite (Nkx3-2 and Twist1) (Nielsen et al., 2001; Tavares et al., 2001), for ectoderm (Fabp3 and Wnt6) (Schubert et al., 2002), for blood (Hbm and Hba1), for endothelial cells (Lmo2 and Sox18) (Minko et al., 2003; Jaffredo et al., 2005; Javerzat et al., 2009; Anderson et al., 2019), for tailbud/pre-segmented mesoderm (Cdx4 and Msgn1) (Joshi et al., 2019), for notochord (Tf and Shh) (Lobjois et al., 2004; Yanai et al., 2005) and for neural crest cells (FoxD3 and Sox10) (McKeown et al., 2005) (Figures 1C, D). Seurat cell cycle scoring determined cell cycle activity. This showed overall the cell clustering was not due to cell cycle phase, although neural cells were predominantly in S phase and G2M phase (Figure 1F).

RNA-tomography profiling of the developing chick embryonic trunk

Spatial information is lost in single cell sequencing data following tissue dissociation (Griffiths et al., 2018). To address this, we next used RNA-tomography to quantify the transcriptomes of a series of individual cryogenic sections along the HH14 embryonic trunk (Figures 2A, B). This enabled a systematic investigation of spatial RNA profiles along the axis and allowed us to resolve the anterior to posterior dimension. Libraries were generated from 20 micron consecutive cryosections of a HH14 chick embryo. Each section was collected into lysis buffer and mRNA captured and amplified using a modification of the G&T-seq protocol (Macaulay et al., 2015). This method allows separating genomic DNA (G) and full-length mRNA (T) from the same sample. Here, we focussed on sequencing transcriptomes. The resulting cDNA libraries had high complexity and enabled us to confidently determine spatial gene expression along the axis (Figure 2C).

Figure 2. RNA-tomography reveals distinct gene expression profiles along the embryonic axis. (A) Stage HH14 chick embryo trunk was sectioned along the anterior-to-posterior axis, from extraembryonic tissue at the posterior end, through the tailbud and pre-somitic mesoderm towards maturing somites. (B) Individual sections were collected in wells followed by RNA isolation and cDNA preparation using section specific barcodes. After that, samples were pooled for linear amplification and sequence library preparation. (C) Spatial expression traces are shown as line plots, left to right represent the posterior to anterior positions. The representative genes shown correspond to the tissue types identified from the scRNA-seq clustering.

Using the same markers as in the previous scRNA-seq analysis (Figures 1C–E), we identified 10 different clusters and established profiles of the lateral plate mesoderm, neural tissue, early somite, maturing somite, ectoderm, blood, endothelial, tailbud, notochord and neural crest cells (Figure 2C). Line graphs indicate spatial patterns of localised anterior-to-posterior restricted gene expression. These were evident for all tissue types with exception of the ectoderm, suggesting that this tissue has few distinguishing markers along the A-P axis. In the most posterior samples, we identified blood (Hbm, Hba1) and endothelial marker genes (Lmo2 and Ramp2), consistent with this region comprising extra-embryonic tissue (Minko et al., 2003; Jaffredo et al., 2005; Javerzat et al., 2009; Anderson et al., 2019). Subsequent sections showed the onset of tailbud genes [T (=brachyury) and Fgf8] (Dubrulle et al., 2001). Some notochord markers were more highly expressed in posterior sections (Shh), while transferrin (Tf) was expressed along the axis, with lower levels around the tailbud region (Riddle et al., 1993; Yanai et al., 2005). Across neural tissue (Ckb, Pax6) (Bhattacharyya et al., 2004; Roy et al., 2013), early somites (Meox1, Tcf15) (Stockdale et al., 2000; Reijntjes et al., 2007; Berti et al., 2015) and maturing somites (Cnmd, Shisa2) (Filipe et al., 2006; Shukunami et al., 2008), gene expression profiles showed a gradual increase towards the anterior. Neural crest cells are beginning to migrate and become distinguishable in the more anterior sections with discreet profiles detected for FoxD3 and Sox10 (McKeown et al., 2005) (Figure 2C).

RNA-tomography resolves mRNA localisation patterns

To identify gene expression patterns systematically, we clustered the spatial gene expression data based on a self-organising heatmap. This sorted the cumulative gene expression traces along a linear axis of 180 profiles and identified 3 major groups of localised mRNA (Figure 3A). The first group of transcripts localised to the most posterior, the second group displayed an increase in the tailbud region and across the pre-somitic mesoderm, and the third group was most highly expressed in the anterior sections comprising epithelial and maturing somites (Figure 3B). Transcripts enriched posteriorly (e.g., Hbm, Epas1, and Lmo2) (Minko et al., 2003; Ota et al., 2007; Anderson et al., 2019) were related to Gene Ontology (GO) terms such as hematopoiesis, erythrocyte differentiation and myeloid homeostasis, consistent with the presence of extraembryonic blood islands (Figure 3C). The second profile showed genes enriched for GO processes such as anterior-posterior pattern specification, embryo morphogenesis and tissue morphogenesis. This overlaps spatially with tailbud and pre-somitic regions. Markers with enriched expression included Wnt5a, Msgn1, Tbx6, and Fgf8 (Dubrulle et al., 2001; Bell et al., 2004; Sweetman et al., 2008) (Figure 3C). The third profile, which overlaps with the formation of somites but also neural tube development, included genes enriched for pattern specification, neurogenesis and animal organ morphogenesis, such as Meox1, Aldh1a2, and Shisa2 (Filipe et al., 2006; Reijntjes et al., 2007) (Figure 3C). Interestingly, the profiles for somites and neural tube are very similar and genes were clustered together, suggesting these tissues mature at a similar rate along the anterior-posterior dimension examined here. However, GO analysis did separate genes associated with generation and differentiation of neurons (Figure 3C). The genes listed for the GO terms, neurogenesis and pattern specification, in k-means cluster 3 are shown in Supplementary Table S1. They include neural markers, such as Wnt4, Nkx6.2, and paraxial mesoderm markers, such as Meox1, Tcf15, and Nkx3.2.

Figure 3. K-means clustering identifies biological components along the posterior-to-anterior axis. (A) Hierarchical cluster analysis of gene expression per section (total 180). Distinct gene expression clusters correspond to different regions along the axis, characterised by extraembryonic tissue, tailbud, pre-somitic mesoderm and epithelial somites—indicated by boxed areas. RNA sequencing reads per gene were normalised against the total read count per section. (B) Spatial expression traces for representative genes in each corresponding cluster, left to right represent the posterior to anterior positions. Cluster 1 represents extraembryonic tissue, cluster 2 tailbud and psm, cluster 3 epithelial somites. (C) Gene ontology on genes enriched in the extraembryonic tissue, tailbud and pre-somitic mesoderm, and epithelial somites and neural tube. (D) Spatial expression traces for signalling pathways associated with anterior-posterior patterning such as WNT (Wnt5A), FGF (Fgf8) and retinoic acid (Aldh1A2), which are expressed in opposing gradients. Spatial expression traces for genes associated with paraxial mesoderm differentiation (Meox1, Tcf15, and Mesp1), with Mesp1 highlighting a discreet region in the psm. Spatial expression traces for Hox genes involved in anterior-posterior patterning (Hoxd11, Hoxc9, and Hoxb1). Expression of Hoxd11 and Hoxc9 displays a clear boundary within the axial region sectioned, while Hoxb1 expression remains high. Line plots for the other Hox genes are in Supplementary Figure S1. (E) Dot plot showing average expression of genes and percentage expressed in each cell cluster associated with the WNT, FGF, NOTCH, and BMP signalling pathways.

RNA-tomography also detected the known opposing gradients across the psm and somites of Wnt5a, FGF8, which are increased posteriorly, versus Aldh1a2, increasing anteriorly and encoding an enzyme involved in retinoid acid synthesis (Olivera-Martinez and Storey, 2007). This is shown in line graph representation of gene expression along the axis (Figure 3D). The transcripts for Meox1 and Tcf15 become upregulated in anterior psm and epithelial somites, whereas Mesp1 transcripts are restricted to an anterior region in the psm, comprising the next but one prospective somite (Figure 3D). Transcripts for Hoxd11, Hoxc9 and Hoxb1 show the expected expression boundaries along the anterior-posterior axis (Figure 3D). Line plots for the members of all four Hox clusters, A–D, show the expected expression cut-off along the axis (Supplementary Figure S1). The expression profiles detected by RNA-tomography for anterior Hox gene expression boundaries and epithelial somite marker genes were in agreement with our previous bulk RNA-seq analysis of paraxial mesoderm tissues (Mok et al., 2021).

Furthermore, the spatial resolution provided by the RNA-tomography overlapped with scRNA-seq expression of components of important signalling pathways (Figure 3E). For example, analysis of WNT and FGF pathways, shows high levels of Wnt5a and Fgf8 in tailbud cells, whilst the receptors, Fgfr1 and Fgfr3, were highly expressed in neural cells (Figure 3E). We detected high levels of Wnt6 in the ectoderm and of Wnt4 in neural cells, consistent with their known expression (Cauthen et al., 2001; Schubert et al., 2002). In addition, the analysis highlighted expression of protein phosphatase, PPP2CA and scaffold protein, SLC9A3R1 in ectodermal cells (Figure 3E). Analysis of NOTCH and BMP pathways confirmed expression of the Notch1 receptor, the glycosyltransferase, Lfng, and the Hes5 transcription factor in neural cells, while neural crest cells expressed Notch1 and Notch2 receptors. The Bmp4 ligand is highly expressed in the lateral plate, while its antagonists follistatin (Fst) and chordin (Chrd) are expressed in early somite and notochord respectively. Expression of the receptor, Bmpr1b increases in maturing somites (Figure 3E).

Correlation of single cell RNA-sequencing with RNA-tomography

Axis patterning is characterized by the progressive differentiation of cell types with anterior-to-posterior identity. To validate genes identified in specific clusters obtained from the single cell RNA-sequencing, we correlated spatial patterns and confirmed gene expression by in situ hybridisation. This provides a visual point of reference for the tomography sections. Expression of Hbm, a haemoglobin gene, is representative for the most posterior group of transcripts, which is enriched for genes involved in hematopoietic differentiation. Consistent with the line graph, we confirm that Hbm transcripts are restricted to blood islands in posterior extraembryonic tissue (Figure 4A). The second spatial cluster values correlate with the tailbud and pre-somitic mesoderm. We validated this data using Wnt5a, which is expressed in the posterior regions at these stages (Baranski et al., 2000) where it regulates cell movement behaviour during axis elongation (Sweetman et al., 2008). In situ hybridisation confirmed that Wnt5a expression was restricted to the tailbud region (Figure 4B), where neuromesodermal progenitor (NMP) cells are located. The third spatial cluster, where gene expression increased towards the anterior, was validated using Tbx22—a key gene for somite boundary formation. Expression of Tbx22 was restricted to caudal somite domains resulting in an in situ signal with periodicity (Figure 4C). All three markers are well characterised (Bell et al., 2004) and help to benchmark the line plots.

Figure 4. Transcriptome map of the embryonic trunk with high spatial resolution. (A–C) UMAP plot and violin plot of representative genes for different tissues, with comparison to Spatial expression trace along the posterior-to-anterior axis. The colour coding in the violin plots corresponds to the clusters identified in the scRNA-seq. Corresponding in situ hybridisation are shown for (A) Hbm—blood islands, (B) Wnt5A—tailbud, (C) Tbx22—caudal somite halves. This provides a visual point of reference when comparing with the spatial expression along the axis shown by the line plots. See also published gene expression patterns in the data base: http://geisha.arizona.edu/geisha/.

Identification of new genes in paraxial mesoderm and somites

Next, we investigated previously unexplored genes uncovered by our approach. We focused on candidates likely to display restricted expression in the pre-somitic mesoderm and somites. Using the subset feature in Seurat, we profiled three clusters from the scRNA-seq dataset—tailbud, early somite and maturing somite (Figure 5A). We re-ran the clustering, findneighbours and pca tests on this new Seurat object. Classic markers for each cluster were re-plotted on UMAPs, such as Msgn1 for the pre-somitic mesoderm, Meox1 and Tcf15 for early somites and Tbx22 for maturing somites (Figure 5A). Sub-clustering revealed restricted expression of follistatin (Fst), in a group of cells potentially representing myogenic cells located in the dorsal part of epithelial somites (Nimmagadda et al., 2005). We identified three genes for further analysis. UMAP plot of Olfml3 and Foxd1 gene expression suggested they are likely to be expressed in maturing somites, whereas Lrig3 was predicted to be expressed in some pre-somitic mesoderm cells, early somites and less in maturing somites (Figure 5A).

Figure 5. Differentially expressed genes in pre-somitic mesoderm and somites. (A) Sub-clustering identifies specific genes within the psm, early somites and maturing somites. Psm is characterised by MSGN1 expression while Tcf15, Meox1 and Tbx22 markers represent somites. Restricted Fst expression may correlate with epaxial myoblasts. (B–D) UMAP, violin plots and spatial expression traces show expression for Foxd1, Olfml3 and Lrig3, not previously identified in somites. The colour coding in the violin plots corresponds to the clusters identified in the scRNA-seq. Violin plots show higher levels of expression in maturing somites for Foxd1 and Olfml3 and a broader distribution for Lrig3. Whole mount in situ hybridisation confirmed the spatially restricted expression of Foxd1 and Olfml3 in epithelial and maturing somites and of Lrig3 in the psm and somites.

Interrogation of the profiles for these genes in the RNA-tomography data showed that expression of Foxd1 and Olfml3 increased towards the anterior regions, with some Foxd1 expression peaks posteriorly. Foxd1 and Olfml3 were identified in spatial cluster 3 (Figure 3). For Lrig3, the spatial data showed increasing expression from the most posterior to the anterior embryonic regions, suggesting that it is expressed from the tailbud to maturing somites. Spatial validation for all three genes by in situ hybridisation in essence confirmed these observations: Foxd1 and Olfml3 are restricted to somites whilst Lrig3 is expressed in the tailbud, pre-somitic mesoderm and in somites. In epithelial somites, all three genes are restricted medially. In maturing somites Foxd1 is broadly expressed, Olfml3 remains medially restricted and Lrig3 is downregulated (Figures 5B–D). We did not detect any posterior expression of Foxd1 by in situ hybridisation, despite the signal in the line plot.

Discussion

The chick embryo is a classic model for developmental biology studies due to the versatility of in vivo experimental approaches (Stern, 2005; Sauka-Spengler and Barembaum, 2008; Gandhi and Bronner, 2018). For example, it has served to better understand the processes of body axis formation, segmentation/somitogenesis and differentiation (Benazeraf and Pourquie, 2013). Here we use single cell transcriptomics and RNA-tomography spatial transcriptomics to map cells that arise in the emerging trunk as it extends. This adds to the growing body of literature, which includes scRNA-seq data of the chick embryo from primitive streak to neurula stages (Vermillion et al., 2018; Williams et al., 2022). Previous work from other labs characterised the molecular signature of neuromesodermal progenitors (NMP) in detail by micro-dissecting anterior PS in stage HH5 and in 6-somite embryos (HH9^-) as well as the tail bud of 35-somite embryos (∼HH18) (Guillot et al., 2021). We did not detect NMPs, which are characterised by co-expression of T-brachyury and Sox2 (Wilson et al., 2009; Henrique et al., 2015; Wymeersch et al., 2021), as a discreet population. This is most likely due to the number of cells sequenced. A recent report examined the anterior most part of the main body axis, including occipital and cervical somites at several stages of development, 4-somites, 7-somites, 10-somites and 13-somites (Rito et al., 2023). This scRNA-seq data identified similar cell populations when the body extends and the cervical-thoracic region forms (HH14, 22-somite embryo) and is in agreement with our data presented here.

Furthermore, we show that a combined analysis of scRNA-seq data with RNA-tomography spatial transcriptomics can reveal novel genes involved in specific aspects of axis extension. As an example, we focused on paraxial mesoderm and discovered the previously unknown relevance of Foxd1, Olfml3 and Lrig3 in developing somites. All three genes showed restricted expression in the medial somite domain suggesting a possible role in early myoblasts, however, this remains to be confirmed with functional approaches, such as gain- and loss-of function experiments. It is noteworthy that Foxd1, a member of the fork-head family of transcription factors, is associated with pluripotency and seems to be required for successful reprogramming (Koga et al., 2014). In addition, Foxd1 protects senescence in human mesenchymal stem cells (hMSC) and is regulated by YAP (Fu et al., 2019). Interestingly, RNA-tomography line graphs indicated Foxd1 expression peaks posteriorly, where progenitors reside. However, this was not confirmed by in situ hybridisation. The reason for this is unclear and could be due to the sensitivity of those different approaches. Not much is known about Olfml3 function in development. It is a secreted glycoprotein of the Olfactomedin-family, which organises the extracellular matrix and has pro-angiogenic properties. Olfml3 deficient mice exhibit abnormalities in the vasculature causing lethality (Imhof et al., 2020). Olfml3 has also been implicated in pre-natal muscle development in pig (Jin and Li, 2019) and in Xenopus it is involved in dorso-ventral patterning by enhancing chordin degradation (Inomata et al., 2008). Finally, Leucine-rich repeats and immunoglobulin-like domains 3 (Lrig3) plays a role in neural crest development (Zhao et al., 2008) in Xenopus. This is consistent with studies in mice, which showed Lrig3 is involved in inner ear morphogenesis by restricting the expression of Ntn1 (Abraira et al., 2008). However, the roles of these genes in developing somites have not yet been investigated.

RNA-tomography is a spatial transcriptomics approach first used in zebrafish embryos (Junker et al., 2014; Kruse et al., 2016). We modified and automated the approach using the G&T-seq protocol (Macaulay et al., 2015) and applied it to the posterior half of a HH14 whole embryo, along the main body axis dimension. As reported previously in the zebrafish heart (Wu et al., 2016; Burkhard and Bakkers, 2018), we obtained high spatial resolution and sensitivity as shown by hierarchical cluster analysis. Known marker genes were expressed in the anticipated spatio-temporal patterns and identified the appropriate regions along the body axis. The extraembryonic tissue was characterized by Hbm, Hbz, Lmo2 and Epas1, the tailbud region by Hoxa13, Wnt5A, the presegemented mesoderm by Msgn1, Tbx6 and epithelial somites by Meox1 and genes involved in retinoic acid (RA) signalling.

This study combines scRNA-seq with spatial information from RNA-tomography improving our understanding and validating the gene expression patterns within the avian elongating body axis and tailbud. While the scRNA-seq provides information on gene expression at the single-cell level, combining it with RNA-tomography bulk sequencing of cryosections allows preservation of the spatial context of gene activity along the anterior-posterior axis. Other technologies are now available for obtaining high-resolution spatial transcriptional profiles of tissues, such as MERFISH, 10X Genomics Visium and Xenium, although these also do not necessarily capture cellular resolution or multiple axial dimensions. Furthermore, all these approaches are more costly to implement, specifically with multiple samples as used here to reveal gene dynamics along the anterior-to-posterior axis. Overall, our results illustrate how combining different approaches can be advantageous when addressing fundamental questions in embryonic development.

Materials and methods

Chicken embryos

Fertilised chicken eggs (Henry Stewart & Co.) were incubated at 37°C with humidity. Embryos were staged according to (Hamburger and Hamilton, 1992). All experiments were performed on chicken embryos younger than two-thirds of gestation and therefore were not regulated by the Animal Scientific Procedures Act 1986.

Preparation of single cells from chicken embryos

The trunk of HH14 embryos were dissected into Ringer’s solution in silicon lined petri dishes and pinned down using the extra-embryonic membranes. Embryonic tissue was transferred into low binding tubes and Ringer’s solution was replaced with Dispase (1.5 mg/mL) in DMEM 10 mM HEPES pH7.5 at 37°C for 7 min prior to treatment with Trypsin (0.05%) at 37°C for 7 min. The reaction was stopped with Ringer’s solution with 0.25% BSA. Cells were spun down and resuspended in Hank’s solution prior to passing through 40 μm cell strainer to obtain single cell suspension.

scRNA-seq library preparation

A suspension of approximately 10,000 single cells was loaded onto the 10X Genomics Single Cell 3′ Chip. cDNA synthesis and library construction were performed according to the manufacturer’s protocol for the Chromium Single Cell 3′ v2 protocol (PN-12033, 10X Genomics). Samples were sequenced on Illumina HiSeq 4,000 100 bp paired-end runs.

scRNA-seq data analysis

Cell Ranger 3.0.2 (10X Genomics) was used to de-multiplex Illumina BCL output, create fastq files and generate single cell feature counts for the library using the Gallus gallus transcriptome (Ensembl release 94) with 82.4% reads mapped to the genome. Subsequent processing was performed using the Seurat v3.1.0 (Butler et al., 2018) package within R (v3.6.1). Cell quality was assessed using simple QC metrics: total number of expressed genes, mitochondrial RNA content and ribosomal RNA content. Identification of chicken mitochondrial RNA content and ribosomal RNA content. Outlier cells were identified if they were above or below three median absolute deviations (MADs) from the median for any metric in the dataset. Data was normalised across all cells using the “LogNormalize” function with a scale factor of 1e4. A set of genes highly variable across the cells was identified using the ‘FindVariableGenes’ function (using “vst” and 2000 features) before being centred and scaled using the “ScaleData” function with default parameters. PCA analysis was performed on scaled data using variant genes and significant principal components were identified by plotting the standard deviation of the top 50 components. The first 2 principal components showed high enrichment for mitochondrial genes and were subsequently regressed and only principal components 3:25 were used to create a Shared Nearest Neighbour (SNN) graph using the “FindNeighbours” function with k.param set to 10. This was used to identify clusters of cells showing similar expression profiles using the FindClusters function with a resolution set to 0.6. The Uniform Manifold Approximation and Projection (UMAP) dimensional reduction technique was used to visualise data from principal components 3:26 in two-dimensional space (“RunUMAP” function). Graphing of the output enabled visualisation of cell cluster identity and marker gene expression. Biomarkers of each cluster were identified using Wilcoxon rank sum tests using Seurat’s “FindAllMarkers” function. It was stipulated that genes must show a logFC of at least 0.01 to be considered for testing. Only positive markers were reported. The expression profile of top markers ranked by average logFC were visualised as heatmaps and dotplots of the scaled data. Cluster identity was determined using visual inspection focussing on the expression of known marker genes. For cells identified in tailbud and somite clusters, the “Subset” function was used to create a new Seurat object which narrowed down to 1,359 cells. Using the “FindNeighbors” feature with dimensions set to 3:20, “FindClusters” resolution of 0.4 and “RunUMAP” set with dimensions 3:16.

RNA sequencing from serial sections—RNA-tomography

Embryos were embedded in Jung tissue freezing medium (Leica), orientated and rapidly frozen on dry ice, and stored at −80°C prior to cryosectioning. Embedded embryos were cryosectioned at 20 μm thickness, collected into 96-well plates (on ice) prior to the addition of 10 μL of RLT plus lysis buffer (Qiagen, Hilden, Germany). All instruments and surfaces were cleaned with 80% v/v ethanol, RNAse-free water and lastly RNAse-out solution after each sample to reduce cross-contamination and RNA degradation. Samples were stored at −80°C until cDNA preparation using the G&T-seq method as previously described (Macaulay et al., 2015) with minor modifications to accommodate the larger volume of lysis buffer. The polyA mRNA capture step of the G&T-seq protocol was used, which enables automated and parallel capture and amplification of RNA from serial sections. The genomic DNA component was not analysed in this study. cDNA was normalised to 0.2 ng/μL before Nextera (Illumina, San Diego, CA, United States) library preparation in a total reaction volume of 4 μL. Libraries were pooled by volume and sequenced on a single lane on the Illumina HiSeq 2,500 (150-bp paired-end reads).

Low-input RNA sequencing analysis

For RNA-seq analysis, we used Refseq version GRCg6a for genome assembly and gene annotation. Reads were trimmed and adapters were removed using tim-galore version 0.4.2. Heatmap based on hierarchical clustering was generated in R-Studio version 1.2.1335 and plotted as a heatmap using the R package DeSeq2 (Love et al., 2014).

Whole mount in situ hybridisation

Whole mount in situ hybridisation using DIG-UTP labelled antisense RNA probes was carried out using standard methods. Probes were generated from amplificons of chicken cDNA using the following primers: Wnt5a (GCA​GCA​CTG​TGG​ACA​ACA​AC/CAC​CGT​CTT​GAA​CTG​GTC​GT), Olfml3 (GGG​AGT​TCA​CGC​TCT​TCT​CG/GAT​GAT​CTG​GTA​GCC​GTC​GT) Hbm (CAT​CAC​ACA​TTG​CCA​CCA​G C/GCA​GCA​ATG​GTG​TCT​TTA​TTG​A), Tbx22 (GGA​TGT​TCC​CAT​CGG​TCA​GG/AGA​CTT​AGC​GCT​CTT CAGGC), Lrig3 (GTC​CTG​ACG​CCT​GGG​AAT​TT/AAT​CTG​TGG​GAC​AGG​ATG​CC), Foxd1 (CCCGC ATC​TCT​AAC​TGT​TAA​GGG/ATT​AAC​CCT​CAC​TAA​AGG​TCC​CTC​AAG​CCT​TCT​CTG​TTC). Briefly, following fixation in 4% PFA embryos were treated with Proteinase K, hybridised with the probe over night at 65°C. After post-hybridisation washed and blocking with BMB (Roche), embryos were treated with anti-DIG antibody coupled to alkaline phosphatase (Merck) and signal developed using NBT/BCIP (Melfords Laboratories).

Data availability statement

The data presented in the study are deposited in the European Nucleotide Archive repository, accession number PRJEB74465 https://www.ebi.ac.uk/ena/browser/view/PRJEB74465.

Ethics statement

The animal study was approved by the AWERB—Animal Welfare and Ethical Review Body. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

GM: Conceptualization, Formal Analysis, Investigation, Validation, Visualization, Writing–original draft, Writing–review and editing. ST: Investigation, Writing–review and editing. ES: Investigation, Validation, Writing–review and editing. LM: Investigation, Writing–review and editing. AL: Investigation, Writing–review and editing. JL: Investigation, Writing–review and editing. VU: Data curation, Formal Analysis, Writing–review and editing. WH: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Supervision, Writing–review and editing. IM: Conceptualization, Funding acquisition, Investigation, Supervision, Writing–review and editing. AM: Conceptualization, Funding acquisition, Project administration, Supervision, Writing–original draft, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. The authors acknowledge support from the Biotechnology and Biological Sciences Research Council (BBSRC), part of UK Research and Innovation, BBSRC Core Strategic Programme Grant (Genomes to Food Security) BB/CSP1720/1 and the National Capability BBS/E/T/000PR9816 (NC1—Supporting EI’s ISPs and the UK Community with Genomics and Single Cell Analysis). WH was additionally supported by the BBSRC Grant BBS/E/T/000PR9818. IM was supported by a BBSRC New Investigator Grant BB/P022073/1. ST was supported by the JIC/NRP International Undergraduate Summer School programme. AM acknowledges BBSRC funding (BB/N007034/1) and a NRPDTP studentship for ES.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2024.1382960/full#supplementary-material

References