Single-cell reconstruction with spatial context of migrating neural crest cells and their microenvironments during vertebrate head and neck formation

During vertebrate development, multipotent neural crest cells emigrate from all regions of the dorsal neural tube, but travel in discrete streams and interact with mesoderm, ectoderm and endoderm tissues throughout their migration to peripheral targets (Etchevers et al., 2019). Proper patterning of the head, neck and heart crucially rely on stereotypical neural crest cell invasion and reciprocal interactions between the neural crest and surrounding tissue microenvironments (Trainor and Krumlauf, 2000; Frisdal and Trainor, 2014). To understand how a healthy individual develops, it is important that we gain insight into how this coordinated movement of cells is controlled. Moreover, the neural crest is a well-known in vivo model to identify and test hypothetical mechanisms described in other cell invasion phenomena, such as cancer metastasis (Kulesa et al., 2013; Maguire et al., 2015; Wislet et al., 2018), as migrating neural crest cells are accessible to in vivo time-lapse imaging, molecular perturbation and single-cell profiling. However, there are serious technical challenges to coordinating the dynamic gene expression changes of highly migratory cells with the spatial context of their collective cell subpopulation in the embryo and the surrounding tissues through which cells interact and travel.

Gene profiling analyses of the neural crest and/or tissues within the embryonic neural crest microenvironments have started to shed light on cell fate restriction and lineage commitment. Bulk-RNAseq of the cranial and trunk neural crest has identified transcription factors important to neural crest cell axial identity (Simoes-Costa et al., 2014; Simoes-Costa and Bronner, 2016; Lumb et al., 2017). Ectoderm analysis has focused on head development at early stages of the neural crest (Plouhinec et al., 2017). Single-cell RNA-seq (scRNA-seq) analyses of either all migrating cranial neural crest cells (Williams et al., 2019; Soldatov et al., 2019) or the first branchial arch (BA1) (Xu et al., 2019) have elucidated cell lineage choices and identified enhancers within a global neural crest gene regulatory network. These data typically collected at a single developmental stage lead to natural questions as to how travel through different microenvironments alters gene expression and how the relationship of gene expression changes with neural crest cell position within a stream.

Cranial neural crest cells maintain a spatially ordered migration such that cells at the leading edge of the invasive front tend to stay in front with only local cell neighbor exchanges (Ridenour et al., 2014; Richardson et al., 2016). Blocking of leaders by a foil barrier (Kulesa et al., 2005) or photoablation (Richardson et al., 2016) revealed that follower cells continue to pathfind to peripheral head targets, suggesting all neural crest cells have the plasticity to respond to local microenvironmental signals. When chick follower neural crest cells were transplanted into the lead position, subsequent imaging and RT-qPCR confirmed the followers pathfind to the second branchial arch (BA2) and rapidly changed the expression of genes to match the lead invasive signature (McLennan et al., 2012). This supported data that neural crest cell behaviors and gene expression are driven by local microenvironmental signals and not hard-wired in premigratory cells. Single-cell RNA-seq profiling of the chick BA2 stream with isolation of leader and follower cell subpopulations at three progressive developmental stages identified molecular heterogeneities within the stream that captured the spatio-temporal transition to directed cell migration and a novel transcriptional signature that is consistent within the lead cells throughout invasion (Morrison et al., 2017a,b). What has remained unclear is a comprehensive spatio-temporal reference map of the transcriptional landscape in BA1-4 during patterning of the vertebrate cranial to cardiac region, and whether there is conservation of the invasion signature within all BA1-4 streams and its relationship to other cell invasion phenomena.

Here, we take advantage of novel single cell technology that overcomes the requirement to isolate specific cell populations of interest, resulting in the analysis of the first four chick branchial arches (BA1-4) at two developmental stages (Hamburger and Hamilton, 1951; HH13 and HH15) and determination of the transcriptional states of migrating neural crest cells and surrounding tissues. We identify cell type composition of the neural crest microenvironments and the spatio-temporal gene expression profile differences associated with the distal (front 20%) and proximal (back 80%) cell subpopulations within BA1-4. We find profiles associated with the molecular transition to the most invasive neural crest cells and identify a distinct transcriptional signature of 266 genes shared within lead cells of all four branchial arch streams. Perturbation of a subset of neural crest cell invasion genes showed significant changes in cell speed, direction, displacement and reduced invasion, suggesting their functional relevance. Strikingly, comparison of the neural crest cell invasion signature with 34 published signatures from a broad range of other cell invasion phenomena identified a shared subset of 45 out of 266 genes. Transplantation and profiling of aggressive human neuroblastoma cells (LAN5, SHSY5Y and NB1643) placed into the chick embryonic neural crest microenvironment revealed a link between invasive ability and upregulation of a subset of the 45 gene panel, providing a resource for downstream functional studies with direct relevance to cancer. Together, these results represent a comprehensive resource of the cellular hierarchy and molecular heterogeneity of the migrating neural crest and BA1-BA4 neural crest microenvironments during vertebrate development.

To characterize neural crest cell migration and differentiation systematically, we sought to derive a comprehensive map of the transcriptional landscape in the first four branchial arches (Fig. 1A,B). We dissected chick BA1 and BA2 at HH13 and BA3 and BA4 at HH15, representing developmental timepoints during neural crest cell migration and branchial arch formation. Each branchial arch microenvironment was subdivided into the front (distal) and the back (Fig. 1A,B). Stage-, branchial arch- and tissue portion-matched tissues from multiple embryos were pooled and barcoded as individual samples for scRNA-seq (10X Genomics Chromium). Bioinformatic analyses (Seurat) of label-free and unsorted branchial arch tissues produced ∼95,000 single cell transcriptional profiles (Butler et al., 2018) (Fig. 1B and Fig. S1). We used Uniform Manifold Approximation and Projection (UMAP) and k-means clustering to display the transcriptomes as seven distinct clusters, containing cell types such as ectoderm, mesoderm and migrating neural crest cells, with each of the clusters marked by increased expression of unique genes (Becht et al., 2019; McInnes et al., 2018) (Fig. 1C-F; Table S1). To ensure our clusters were not the product of read depth or cell cycle variations (Kiselev et al., 2019; Whitfield et al., 2002), we performed several quality control analyses (Fig. S1) and find continuity in both sequencing depth and, surprisingly, cell cycle markers among our seven clusters. Feature plots of the mean normalized expression values for lists of genes known to be associated with different cell cycle phases are shown (Fig. S1). The evenness of the expression of these sets of genes indicates a lack of cell cycle-driven expression.

To fully characterize these seven clusters, including the spatial location of the clusters within the context of the developing embryo conferred by the eight spatially unique samples (Fig. 1A,B), we determined the contribution of each spatially distinct sample to each cluster (Fig. 1G; Table S2). The paraxial mesoderm cluster is dominated by cells from the back of BA3 and BA4 and most of the cardiac mesoderm cluster is made up of cells from BA3 and BA4 (Fig. 1G). Contributions to endothelial, ectodermal, SOX10⁺ neural crest, glial and mesodermal cells are observed from all branchial arches analyzed. To better understand the types of cells and tissues that exist at each anatomical location within BA1-4, we calculated the percentage composition of each of the eight scRNA-seq samples (Fig. 1H). The front of each branchial arch contained a higher percentage of ectodermal cells than the back. Conversely, the back of each branchial arch had a higher percentage of glia and cranial mesoderm than the front. The percentage of endothelial cells was consistent across all samples. As noted in Fig. 1G, paraxial mesoderm was almost exclusive to the back samples from BA3 and BA4, and cardiac mesoderm was preferentially found in cells from BA3 and BA4 (Fig. 1H).

To better understand similarities and differences along the axis within BA1-4, we highlighted pre-otic and post-otic cells within the context of the seven cluster UMAP (Fig. 2). We visualized the contribution of each spatially restricted sample against the entirety of the UMAP (Fig. 2D-I) and found that BA1 and BA2 cells, and BA3 and BA4 cells supplied ectoderm, otic placode, cardiac mesoderm and endothelial clusters. For example, endothelial cells are equally distributed among pre- and post-otic regions (Figs 1G,H and 2A-C,E,H). Conversely, cells within SOX10+ neural crest, cardiac mesoderm and paraxial mesoderm clusters segregate based upon their pre-otic and post-otic origin (compare differences in black and magenta regions in Fig. 2A-C). Cranial mesoderm marked by CYP26C1 showed a dominant contribution from pre-otic cells, as previously observed (Bothe et al., 2011). MEIS2 expression is specific to post-otic cells, especially in SOX10⁺ neural crest and paraxial mesoderm (Fig. 2A-H).

In a previous study, we identified distinct leader and follower neural crest molecular profiles within the BA2 stream (Morrison et al., 2017a). For comparison, we distinguished proximal and distal subregions of the arches within the context of the seven-cluster UMAP. There is a significant overlap between the front and back cells in ectoderm and endothelial clusters (Fig. 2E,H), suggesting a similar transcriptional signature for front and back cells within these clusters. In contrast, front and back cells displayed some spatial segregation within the cluster of SOX10⁺ neural crest cells at each branchial arch level (Fig. 2E,H). For example, CST3 preferentially marks front cells of the cardiac mesoderm (cluster 4) and expression of the neural marker PRTG is restricted to back cells within most clusters (Fig. 2I).

To determine whether our scRNA-seq profiling faithfully captured the typical HOX gene expression patterns of the hindbrain and neural crest (Trainor and Krumlauf, 2000), we analyzed the expression of several well-known markers (Fig. 3A-D). We find that HOXA2, HOXA3 and HOXB3 expression correctly correlated with the branchial arches (Fig. 3A-D). This faithful capture of HOX genes and other branchial arch specific genes (DLX5 and MAFB; Fig. 3E,F) led us to examine markers of migrating neural crest cells. Migrating neural crest cells have traditionally been identified by the expression of SOX10 or ITGB3 (Tucker et al., 1984; Vincent et al., 1983; Southard-Smith et al., 1998; Pietri et al., 2003). Our scRNA-seq analysis identified SOX10 and ITGB3 (Fig. 3G). We also identified novel markers of migrating neural crest cells, including UBA7, ITPK1 and COL18A1 (Fig. 3H; Table S1). Together, these data verified our spatially distinct tissue isolations, confirmed expression of known branchial arch and migrating neural crest markers, and provided novel candidates for exploring the temporal dynamic signatures in migrating neural crest cells.

We next explored whether the novel neural crest cell invasion signature associated with the most invasive neural crest cells migrating into BA2 (Morrison et al., 2017a) is conserved in all BA1-4 streams. As SOX10 and other commonly used neural crest markers are reduced in the most invasive neural crest cells at the migratory front (Fig. S2A-D; Morrison et al., 2017b), we asked whether these cells could be more accurately identified within the UMAP of BA1-4 tissue by employing the known BA2 cell invasion signature (Morrison et al., 2017a) as a reference. The cardiac mesoderm and endothelial clusters (Fig. 1C; clusters 4 and 6, respectively) were enriched for an average of the 964 genes enriched in the BA2 invasion signature (Fig. 3I). As cluster 6 distinctly expressed well-characterized markers of endothelial cells, such as SOX18, LMO2 and CDH5 (Fig. 1C-F; Fig. 3I; Table S1), we instead focused our efforts to identify invasive neural crest cells within the cardiac mesoderm cluster. As the BA2 invasion signature was originally determined by comparing lead cells within BA2 with the other neural crest cell streams, enrichment of BA2 invasion genes in the cardiac mesoderm cluster should be compared with expression within the SOX10⁺ neural crest cluster (Fig. 1C; cluster 1).

We find that, unlike the cardiac mesoderm cluster, the SOX10⁺ neural crest cluster is not enriched for invasion genes, suggesting it may correspond to a proximal rather than distal subpopulation of migrating neural crest cells (Fig. 3I). To further verify similarity of the cardiac mesoderm and BA2 invasion signature, we examined the average expression of 406 genes reduced in the invasion signature (Fig. 3J) and merged (Fig. 3K). There are clearly fewer genes reduced in the invasion signature in the cardiac mesoderm cluster compared with the SOX10⁺ neural crest cluster 1 (Fig. 1C; cluster 4 and Fig. 3J). These results suggest that the most invasive neural crest cells at the migratory stream front are not part of the SOX10⁺ neural crest cluster (cluster 1), which appear to represent the subpopulation of migrating neural crest cells in proximal positions within a stream.

The subpopulation of cells expressing the novel invasion signature within the BA2 stream was a small percentage (∼2%) of cells within the entire migratory stream that were confined to the invasive front (Morrison et al., 2017a). Here, our analysis of all BA1-4 tissue found a larger percentage (8%) of the cells that clustered as cardiac mesoderm (Table S2). The cardiac mesoderm cluster (cluster 4) also contained cells from the back (proximal) region of BA3 and BA4 (Fig. 1G,H; Table S2). These observations suggested that cluster 4 may contain both invasive neural crest in addition to migrating mesodermal cells. To explore the molecular heterogeneity within the cardiac mesoderm cluster, we independently subclustered only cells from the cardiac mesoderm cluster, and color-coded the resulting single-cell transcriptomes by sample types that confer both spatial and temporal information (Fig. 4C; Fig. S3). That is, we re-ran dimensionality reduction on cluster 4 alone. We found that single cell transcriptomes on the left side of the sub-clustering composition UMAP for cluster 4 almost exclusively comprised front and back BA4 cells (Fig. 4C; Fig. S3). Single cell transcriptomes on the right side of the plot contained front cells from all four branchial arches and very few back cells.

To provide evidence that subclusters 1 and 3 within cluster 4 comprise NC cells, we asked how many of our 314 NC markers (cluster 1 markers from Table S1) were also markers of each of the four cluster 4 subclusters (Fig. S1J; Table S4). We found that none of the 64 or 147 markers of subclusters 0 or 2 (presumptive cardiac mesoderm) contained markers of NC cells. Conversely, 147 of the 390 markers (45%) of subcluster 1 (presumptive NC) also marked cluster 1 NC. Furthermore, 62 of the 197 (31%) markers of subcluster 3 (presumptive invasive NC) overlapped with the markers of cluster 1 NC (Fig. S1K,L).

To further confirm that cluster 4 contained the most invasive neural crest cells, we first examined the expression of the BA2 invasion signature, which contained 964 enriched genes compared with other neural crest cells (Morrison et al., 2017a). We find the average of BA2 cell invasion markers on the right side of the cardiac mesoderm UMAP (Fig. 4A; circled), while the bottom region of the cardiac mesoderm cluster is enriched for the average of all SOX10⁺ neural crest cell markers (Fig. 4B; circled). We used the overlap of the enriched BA2 invasion signature and neural crest cell markers with the cardiac mesoderm cluster to assign a subpopulation as BA1-4 invasive neural crest cells (Fig. 4C; intersection of circled subregions). To determine whether the putative BA1-4 neural crest invasion signature is present at each axial level and whether they exist in the more proximal or distal regions of each neural crest cell migratory stream, we used the spatial information preserved within the scRNA-seq dataset. We find that the BA1-4 neural crest invasion signature is almost exclusively at the front of BA1-3 and preferentially at the front of BA4 (Fig. 4D). This analysis of BA1-4 neural crest cell invasion signature approximated the same front-back ratio previously observed for the BA2 invasion signature (Morrison et al., 2017a).

Having identified BA1-4 invasive neural crest cells, we sought to establish their molecular properties. To identify genes unique to the BA1-4 invasion signature, we first compared BA1-4 invasion genes to all other BA1-4 cells analyzed by scRNA-seq and found 405 genes enriched in the BA1-4 invasive neural crest cells (Fig. 4E; pink circle). Next, we compared the BA1-4 neural crest invasion signature with the transcriptionally similar cardiac mesoderm cells and identified 313 genes enriched in BA1-4 invasive neural crest cells (Fig. 4E; blue circle). At the intersection of these two gene lists are 266 definitive BA1-4 invasive neural crest cell markers (Fig. 4E; magenta circle, Table S3) and the top 20 out of 266 genes are shown (Fig. 4F). This analysis clearly demonstrates that there is a shared invasion signature in all four branchial arch streams, indicating a conserved property of migrating cranial-to-cardiac neural crest cells.

The above results led us to examine the reduction of multipotency markers in cranial-to-cardiac neural crest cells by analyzing the markers of subpopulations of cells between the proximal, SOX10⁺ NC and most invasive NC (Fig. 4G; subpopulations X, Y and Z). We find that traditional neural crest cell markers are downregulated when comparing clusters 1, X, Y, Z and cluster 4 (Fig. 4H). Strikingly, more than 90% of subcluster X comes from the front and back of BA1 (Fig. 4I). About 80% of the cells in subcluster Y come from the back of BA3 and BA4, and ∼70% of subcluster Z is composed of cells from the back of BA4. These results describe transcriptionally unique subpopulations within individual axial levels. Furthermore, fewer than 6% of cells in X, Y and Z originated in BA2, suggesting these interesting subpopulations may have escaped detection within our earlier study that analyzed smaller cell numbers of the BA2 stream (Morrison et al., 2017a). Overall, these results present an attractive hypothesis for future testing; cells within X, Y and Z represent axial-level-specific transitions between proximal neural crest cells with high SOX10 expression and distal leading-edge cells.

To continue to validate the expression of novel neural crest cell invasion genes within the branchial arch streams and their spatial expression, we leveraged our expertise in combined multiplexed fluorescence in situ hybridization (RNAscope), immunohistochemistry, and tissue clearing and image segmentation methods (McLennan et al., 2015; Morrison et al., 2017a,b). We clearly identify the spatial heterogeneity of expression of selected markers within the BA1-4 streams by selecting CDH5 (expressed within follower neural crest cells), HNK-1 (all migrating neural crest) and KAZALD1 (leaders) (Fig. 5F). In higher resolution of a typical branchial arch stream (for example, BA2), we observed KAZALD1 distinctly visible within the lead neural crest cells, confirming our earlier observations of KAZALD1 expression in individual lead NC cells (Morrison et al., 2017a). SOX10 expression was strongly expressed within the migrating neural crest cells, but absent in the leaders (Fig. 5G), as was observed (Morrison et al., 2017b; Fig. 1) and quantified by both polyline kymograph measurements of the BA2 stream (Morrison et al., 2017b; Fig. 3) and specifically within migrating BA2 neural crest cells (Morrison et al., 2017b; Fig. 4). Surprisingly, we observed KAZALD1 expression in the tissue distal to the neural crest stream, which corresponded to lead cells within the mesoderm (Fig. 5G).

Previous morpholino knockdown of a subset of neural crest invasion genes led to a reduction in distance migrated by the BA2 stream in vivo (Morrison et al., 2017a), suggesting a functionally significant role. Indeed, when we performed a detailed investigation of the functional role of the novel neural crest cell invasion gene AQP1 (McLennan et al., 2015), we learned that its reduction dramatically affected the initiation and stabilization of neural crest cell filopodia, resulting in a loss of local matrix degradation and attachment that affected the ability of the cell to properly invade (McLennan et al., 2020). To continue to assess the functions of novel neural crest cell genes (Morrison et al., 2017a), but this time selected for expression within all of the BA1-4 migratory streams (Fig. 4E; total of 266 genes), we knocked down a small subset of individual genes (PODXL, CST3 and KIF26A) and measured in vitro changes in cell speed and displacement from tracked cell trajectories in neural tube explanted cultures (Fig. 5H,I). Loss-of-function of cystatin 3 (CST3) significantly decreased mean cell speed, increased the mean cell speed with loss of podocalyxin-like protein 1 (PODXL) (Fig. 5H) and increased the displacement of neural crest cells with the loss of kinesin family member 26A (KIF26A) (Fig. 5I). Moreover, and to build on the success of observing dramatic changes to the in vivo distance migrated by cranial neural crest cells into BA2 after in vivo loss of function (Morrison et al., 2017a,b; Fig. 6), we determined that knockdown of either PODXL or KIF26A by siRNA led to significant reduction of migration into BA2 (Fig. 5J). Thus, loss of function of a subset of novel genes enhanced within lead neural crest cells showed significant changes in cell migratory characteristics in vitro, and in vivo invasion towards the second branchial arch (BA2).

Cell invasion is a hallmark in many different biological phenomena, including embryogenesis, wound repair, the immune response and cancer metastasis. To better understand the relationship of the neural crest invasion signature of BA1-4 (Fig. 4E) with other processes that share similar physiological and pathological features, we curated a gene list based on 34 published cell invasion signatures (Fig. 6A,B; Table S5). To reduce bias, we ensured that the invasion signatures represented a broad range of cell types, techniques and model organisms. Because of the wide array of model organisms represented, we first ensured that the enriched genes curated from the 34 publications have an ortholog in chicken. We find that 28 out of the 34 signatures (82%) displayed overlap of enriched genes with the BA2 neural crest invasion signature (Morrison et al., 2017a). There are 252 genes enriched in the neural crest invasion signature that appear in at least one of the 34 published invasion signatures (Table S5). Increasing the stringency of analysis revealed a group of 45 neural crest invasion genes in two or more other invasion signatures and 15 genes enriched in three or more invasion signatures (Fig. 6C,D). Two genes (JUN and ITGB1) are enriched in four invasion signatures and one gene (ITGB1) is enriched in five invasion signatures (Fig. 6D). We performed canonical pathway enrichment and REVIGO analyses to explore specific functions and cellular states associated with the 45 gene invasion signature (Fig. 6E; Fig. S2E) similar to Supek et al. (2011). This condensed signature, although not encapsulating all genes necessary for invasive characteristics, further supports the hypothesis that a common set of genes is required to achieve efficient collective cell invasion.

To directly test whether the 45 gene cell invasion signature determined above has direct relevance to human cancer cell invasion, we took advantage of our established chick embryo transplant model (Kulesa et al., 2006; Bailey and Kulesa, 2014). We transplanted three neural crest-derived human neuroblastoma cell lines (LAN5, SHSY5Y and NB1643) into the chick trunk embryonic neural crest microenvironment (Fig. 6F). After re-incubating eggs for 48 h, we typically observed that a subpopulation of neuroblastoma cells invaded the chick embryonic microenvironment and a subset of cells remained at the transplant site (Fig. 6F). LAN5 and SHSY5Y human neuroblastoma cells were highly invasive into the embryonic neural crest microenvironment (Fig. 6F), but very few NB1643 cells exited from the dorsal neural tube transplant site (data not shown). To evaluate changes in gene expression, we manually isolated the chick trunk tissue containing the invasive human neuroblastoma cells and used FACs and RNA-seq analysis to compare gene expression changes to each of the cultured cell subpopulations (LAN5, SHSY5Y and NB1643) (Fig. 6F,G). We find that the LAN5 human neuroblastoma cells were the most invasive into the chick embryonic neural crest microenvironment and these cells upregulated the largest percentage (15 out of 45; 33%) of the 45 gene panel (Fig. 6C,G). Moderately invasive SHSY5Y cells upregulated a modest number of the 45 gene panel (six out of 45; 13%) (Fig. 6G). In contrast, the poorly invasive NB1643 cells tended to remain at the transplant site, and upregulated only four out of the 45 gene panel (less than 10%) (Fig. 6G).

In comparing the profiles of the aggressive versus non-aggressive human neuroblastoma cells after transplantation into the chick embryo, we noted several interesting gene expression differences (Fig. 6G). First, six out of 15 genes upregulated in the LAN5 cells were the same genes upregulated in the invasive SHSY5Y cells (Fig. 6G). This included actin cytoskeletal components ARPC3 (Arp2/3-related), TUBB6 (tubulin-related) and the structural protein vimentin (VIM). Invasive LAN5 and SHSY5Y cells also showed upregulation of Vascular endothelial growth factor (VEGF) receptors, FLT1 (also known as VEGFR1) and NRP1 (Fig. 6G, green); VEGF is present in the chick embryonic microenvironment and is a neural crest cell chemoattractant (McLennan et al., 2010). NB1643 cells did not upregulate VEGF receptors, but included four genes in common with the LAN5 cells (Fig. 6G, red). These results suggest that the invasive ability of human cancer cells within the embryonic neural crest microenvironment may correlate with the number and type of genes upregulated from a common cell invasion signature.

We have molecularly characterized migrating neural crest cells and their surrounding microenvironments within the first four chick cranial to cardiac branchial arches (BA1-4) using novel label-free, unsorted single-cell RNA sequencing. These data provide a comprehensive spatio-temporal map of the transcriptional landscape during patterning of the vertebrate head and neck. The faithful capture of HOX genes and other branchial arch-specific markers led to identification of novel genes within axial level-specific and proximal-to-distal subpopulations of neural crest cell migratory streams. This included UBA7, ITPK1 and COL18A1 as more definitive markers of migrating neural crest cells from neural tube exit through branchial arch invasion. In-depth analysis of lead neural crest cells identified a novel cell invasion transcriptional signature common to all BA1-4 streams. Knockdown of a subset of novel neural crest cell invasion genes (PODXL, CST3 or KIF26A) and either in vitro time-lapse analysis or in vivo 3D confocal imaging revealed some significant changes in neural crest cell migratory behaviors and reduced invasion towards BA2. The neural crest cell invasion signature shared 45 genes in common with published invasion signatures curated from a wide range of other cell invasion phenomenon and has a direct relevance to human neuroblastoma cancer cell metastasis.

Our data support a model of axial level and spatio-temporal molecular heterogeneities within the neural crest cell migratory streams and microenvironments of the first four branchial arches that may be mined for insights into cell fate decisions. Previous knowledge derived from our sc-qPCR and scRNA-seq profiling of purified neural crest cells combined here with the innovative label-free, non-isolation cell profiling of the entire branchial arches provided confidence in mapping cell type composition (Figs 1 and 3). These results also provide a much richer dataset with which to examine both neural crest and neural crest-microenvironment signaling, rather than bulk RNA-seq analysis of the tissues or purified neural crest cells. By providing a knowledge base for discovery of novel cell type markers in space and time, future work will be able to identify distinct spatio-temporal locations corresponding to neural crest cell hierarchy.

The discovery that the most invasive neural crest cells from all BA1-4 streams clustered with mesodermal cells (Figs 1 and 2) and share expression of genes (Fig. 5) suggests a high degree of transcriptional similarity between cells of two different lineages. The co-clustering of invasive neural crest and mesodermal cells appears counterintuitive when viewed in the context of differentiation, but we suggest that the transcriptional similarity between these two progenitor populations lies in their migratory ability, an observation that was previously reported in zebrafish (Wagner et al., 2018). Mapping of the genes enriched in the neural crest-refined signature to specific signaling pathways will provide for future in-depth functional experiments that begin to tease out the complex network dynamics that underlie collective cell migration of the neural crest and mesoderm, and cellular hierarchy of the neural crest. The combination of high-throughput screening of cultured neural crest or neural crest-derived cancer cells with the in vivo neural crest model, both of which leverage the strengths of dynamic imaging, present a clear means with which to address the fundamental mechanisms that underlie cell communication and collective cell migration.

The above findings may have important implications for neural crest-derived cancers (melanoma and neuroblastoma) and other aggressive cell phenomena. The molecular signature of the neural crest compared well with other cell invasion phenomena (Fig. 6), suggesting the exciting possibility that fundamental properties underlie development, wound repair and cancer metastasis. Not all cell invasion signatures that we curated from 34 published profiles included a complete transcriptome analysis or lacked an exact orthologue between species (Fig. 6B). Nevertheless, the overlap of 45 out of 964 genes does provide unique molecular inroads, and the in vivo neural crest model is poised to study gene network dynamics and function. This is exemplified by our results that show human neuroblastoma cells transplanted into our established chick embryo model upregulated neural crest invasion genes according to invasive ability (Fig. 6F,G), offering a rapid means to test gene candidates in common with other human cancers. These types of studies have the potential to derive more accurate gene lists for risk assessment of metastasis and to test therapeutic targets to either promote (wound repair or immune response) or inhibit (cancer) cell invasion. In conclusion, these data provide a comprehensive resource of the dynamic transcriptional landscape in BA1-4 during patterning of the chick cranial to cardiac branchial arches, and offer a resource for comparative analysis with other vertebrate organisms, and studies of cell invasion and patterning.

All experiments were performed according to the Stowers Institute for Medical Research's Institute Biosafety Committee (IBC-2003-23-pmk) and federal ethical standards. Fertilized white leghorn chicken eggs (NCBI Taxonomy ID:9031; Centurion Poultry, Lexington, GA, USA) were incubated at 38°C in a humidified incubator until the desired Hamburger and Hamilton (HH) stage of development (Hamburger and Hamilton, 1951). Embryos were screened for health and harvested into chilled 0.1% DEPC phosphate-buffered saline (PBS). To capture neural crest during active migration, branchial arches (BA) 1 and 2 were manually isolated from HH13 embryos (n=15 embryos). Branchial arches 3 and 4 were manually isolated from HH15 embryos (n=10 embryos). Each branchial arch was further manually dissected into the front 20% (distal-most) and 80% back (proximal-most) regions. Stage-, branchial arch- and tissue region-matched tissues were pooled, dissociated as previously described (Morrison et al., 2015, 2017a) and treated as individual samples. The viability and concentration of each single cell suspension was quickly confirmed on a Nexelome Cellometer Auto T4 (Nexelome Bioscience) and the cell suspensions used as input for 10X scRNA-seq (10X Genomics) according to the manufacturer's recommendations.

Dissociated cells were loaded on a Chromium Single Cell Controller (10X Genomics), based on live cell concentration, with a target of 4000-10,000 cells per sample. Libraries were prepared using the Chromium Single Cell 3′ Library and Gel Bead Kit v2 (10X Genomics) according to the manufacturer's directions. Resulting short fragment libraries were checked for quality and quantity using a Bioanalyzer 2100 (Agilent Technologies) and Qubit Fluorometer (Invitrogen). Libraries were pooled and sequenced to a depth necessary to achieve 25-35,000 mean reads per cell, 125-520 M reads each, on an Illumina HiSeq 2500 instrument using Rapid SBS v2 chemistry with the following paired read lengths: 26 bp Read1, 8 bp I7 Index and 98 bp Read2.

Eight samples of 10X Genomics Chromium scRNA-seq data were sequenced on five Illumina HiSeq 2500 flowcells. Data were processed with bcl2fastq (2.20) and aligned and aggregated with CellRanger (2.1.0). Data were aligned to galGal4 using annotations from Ensembl 84 to be consistent with our previous scRNA-seq datasets (Morrison et al., 2017a). Downstream analysis was carried out in R (3.5.0) with the Seurat package (2.3.4) (Butler et al., 2018). Cells were kept for downstream analysis if they had more than 500 genes expressed, fewer than 20,000 UMIs and less than 30% mitochondrial expression. 2796 genes were selected as variable genes for downstream analysis, and the first 50 principal components were used for later steps (UMAP and cluster identification). Clusters were identified by k-means. For identifying cluster markers, we used FindAllMarkers with the minimum percentage of cells within a cluster expressing the gene set to 25% and a log fold change threshold of greater than or equal to 0.25. To further investigate some clusters, we performed subclustering by subsetting the Seurat object by cluster and running FindClusters at higher resolutions. UMAP has gained favor in recent single cell analyses because it retains more of the global architecture of the dataset (McInnes et al., 2018; Becht et al., 2019) and information may be gathered from the spatial relationship of single-cell clusters (Becht et al., 2019). Canonical pathway enrichment analysis was determined by Qiagen's Ingenuity Pathway Analysis (IPA). Cell cycle analysis was performed using 1134 cell cycle regulated clones published by Whitfield et al. (2002). From the 1134 cell cycle markers, we used a list of 276 genes with matching gene symbols between human and chick (Galgal4; Ensembl release 84). The Whitfield 2002 cell cycle markers are a field standard that have been used in many scRNA-seq publications (e.g. Macosko et al., 2015, Cell) and have also been shown to correlate with FUCCI cell cycle phases (Hsiao et al., 2020, Genome Res). Feature plots of the mean normalized expression values for lists of genes known to be associated with different cell cycle phases are shown (Fig. S1). We estimated the multiplet rate according to the 10X Genomics webpage: https://kb.10xgenomics.com/hc/en-us/articles/360001378811-What-is-the-maximum-number-of-cells-that-can-be-profiled.

Our invasive neural crest cells from the first four branchial arches (BA1-4) were identified using our previously established BA2 cell invasion signature that was rigorously characterized by Morrison et al. (2017a,b). First, cells that produced the BA2 invasion signature were definitively neural crest, based upon fluorescent labeling and cytometric isolation. Second, the cells were sequenced to the maximum depth available (SMART-seq v3), which is far more sensitive than current droplet-based techniques. We confirmed similar transcriptional signatures in two types of highly migratory cells that exist within the same spatiotemporal space by detection of migratory genes in both HNK1-positive (unique identifier of migrating avian neural crest cells) and invasive neural crest cells, as well as neighboring HNK1-negative mesodermal cells (Fig. 5). Furthermore, the transcriptional similarity of subsets of NC and mesodermal cells has been previously described in zebrafish by Wagner et al. (2018).

Integrated analysis of gene expression and protein detection in single cells was performed as previously described (Morrison et al., 2017a,b). Briefly, RNAscope (Advanced Cell Diagnostics) fluorescent in situ hybridization and subsequent immunohistochemistry were carried out on fixed avian embryos. Images were acquired on a Zeiss LSM 800 laser scanning confocal microscope and analyzed in Imaris (Bitplane). After the 3D boundaries of each neural crest cell were determined by the membrane-specific HNK-1 immunohistochemistry signal (carbohydrate epitope localized to the surface of migrating avian NC cells), the number of RNAscope spots within the 3D volume of each neural crest cell were quantified. Our RNAscope analysis of invasion genes could not confirm gene expression unique to cells within the invasive front, suggesting the transcriptional signature is a cohort of genes rather than enhanced expression of any single gene.

In vitro neural tubes cultures and cell tracking analysis were performed as previously described (McLennan et al., 2010, 2020). 2 h prior to neural tube removal from the embryos, the dorsal region of each neural tube was transfected with fluorescein-tagged splice modifying morpholinos targeting the invasion genes using in ovo electroporation as previously described (McLennan and Kulesa, 2019). Briefly, we control for the timing and position of the morpholino injection into the dorsal neural tube and introduction into premigratory neural crest cells, which are well documented in our previous work and by others in the avian community. The 5′ to 3′ morpholino sequences are CGATTTTAATCTTCTGATACCTGCT for PODXL, ACATCCTGCTCAGAGCCTACCTTAG for CST3, AACAGAAAGGTGACAAACCTGATGA for KIF26A, ATCTGAGGCTCTGGGAATGGAAGAT for KAZALD1 and GACGGACCTCTAAGGCTCCTCACC for RBM38 (Gene Tools).

siRNAs were designed towards chick KIF26A (5′-AAGGUCCGGAGUUUGUUCUU-3′) and chick PODXL (5′-AAGACCAACCGAAUGCUGUAGUU-3′), and purchased from Horizon Discovery Biosciences (Cambridge, UK). siRNAs were injected and electroporated into neural tubes of HH8 chick embryos, reincubated for 24 h and stained for HNK-1. Harvested embryos were mounted on a glass coverslide and imaged by 3D confocal microscopy on an LSM 800 (Zeiss). Cell migratory distance was measured for siRNA and HNK-1-labeled neural crest cells at the axial level of r4. Statistical analysis of cell tracks was calculated using Student's t-test with two tails and two sample unequal variance. Statistical analysis and P-value calculations in Fig. 5 were performed using Bonferroni's method.

Human neuroblastoma cell lines, LAN5 and NB1643 (Children's Oncology Group, Texas Tech University, USA), and SHSY5Y (ATCC, Manassas, VA, USA; CRL-2266), were maintained as described by Kasemeier-Kulesa et al. (2018). Neuroblastoma cells were cultured in hanging drops and transplanted into Hamburger Hamilton (HH) stage 10 chick embryos (Kulesa et al., 2006). As we observed that only a few of the least invasive human neuroblastoma NB1643 cells invaded the embryonic neural crest microenvironment, single-cell RNA-seq was not deemed a good option for this analysis.

We thank Dr Robb E. Krumlauf for critical reading of the manuscript and the Stowers Technology Centers for technical assistance.

The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.199468.