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.

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).

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).

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).

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.

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).

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.

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).

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.

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.

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.

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.

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.

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).

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 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).