Neural crest migration is driven by a few trailblazer cells with a unique molecular signature narrowly confined to the invasive front

In examples that range from primitive streak formation to mechanosensory organogenesis, several embryonic cell populations undergo persistent, directed migration in coordinated groups (Tarbashevich and Raz, 2010; Voiculescu et al., 2014; Piotrowski and Baker, 2014). When migratory cells fail to reach a target or populate an incorrect location, this often leads to improper cell differentiation or uncontrolled cell proliferation. Thus, studies of embryonic cell migration mechanisms are important for a better understanding of birth defects and tumor formation.

One of the key features of embryonic cell migration is the persistent, directed movement of cells in multicellular streams. During multicellular streaming, a cell autonomously controls its cytoskeleton but moves with its neighbors as a population. A long-standing question concerns the nature of the mechanisms that regulate the persistence of direction and cohesion of multicellular streams. Large-scale genomic analyses of premigratory or migrating embryonic cells, such as the neural crest (NC), have shed light on the genes expressed in migratory versus non-migratory cells during embryonic development (Gammill and Bronner-Fraser, 2002,, 2003; Molyneaux et al., 2004; Adams et al., 2008; Gallardo et al., 2010; Simoes-Costa and Bronner, 2013; Simoes-Costa et al., 2014). However, what remains unclear is how gene expression varies as migrating cells respond to different microenvironments and how this transduces into observed cell behaviors. Therefore, investigative efforts that correlate molecular interrogation with in vivo cell behavior analyses will yield important insights into embryonic cell migration events at the level of both individual cells and the population.

The NC is one of the most striking examples of long-distance embryonic cell migration that is accessible to manipulation and in vivo observation. NC cells emerge from the dorsal neural tube and are sculpted by intrinsic and extrinsic signals into discrete, multicellular streams throughout the head and trunk (Kulesa and Gammill, 2010). Analysis of cell behaviors from in vivo time-lapse imaging in chick (Kulesa and Fraser, 1998, 2000; Teddy and Kulesa, 2004; Kulesa et al., 2008; McKinney and Kulesa, 2011; Ridenour et al., 2014) and in chick and mouse intestine (Young et al., 2004, 2014; Nishiyama et al., 2012) has shown that there are regional differences in cell speed, direction, proliferation, calcium activity and cell morphology depending on cell position within an NC migratory stream. We previously showed that a computational model of NC cell migration, using the chick cranial NC cell behavioral data, predicted that successful cell persistence of direction and stream cohesion would result from the presence of unique lead and trailing subpopulations (McLennan et al., 2012). Further, tissue transplantation studies in which lead cells were placed into the trailing stream or vice versa showed that NC cell behaviors and gene expression profiles are not predetermined but depend on stream position. What remains unclear is how the size of the lead NC subpopulation and changes in its molecular profile affect the persistence of direction and stream cohesion.

To address these questions, we first isolated NC cells from eight discrete subregions of a typical migratory stream in the chick head using laser capture microdissection (LCM) and analyzed the expression of 77 genes using real-time quantitative PCR (RT-qPCR). We examined regional differences in gene expression and used a newly emerged fluorescent hybridization chain reaction (HCR) strategy for multiplexed imaging of mRNA expression (Choi et al., 2010, 2014). We next isolated and profiled single lead NC cells (within the most distal portion of the invasive front) at two distinct phases of migration. We compared single-cell gene expression of leaders with cells isolated from the entire stream. We identified a molecular signature unique to these selected lead NC cells (which we define as ‘trailblazers’) narrowly confined to the invasive front. We then used gain- and loss-of-function of transcription factors upstream of the trailblazer molecular signature to determine whether genes expressed as part of this signature are crucial to NC cell migration. These experiments were performed in parallel with computational modeling that simulated our experimental scenarios. Our results highlight the molecular heterogeneity of cells during cranial NC cell migration and the requirement for only a small subset of trailblazer cells to ensure persistence of direction and stream maintenance.

We previously used computational modeling, cell morphometrics and broad molecular profiling to reveal the existence of at least two subpopulations of cells with distinct molecular profiles within a typical cranial NC cell migratory stream (McLennan et al., 2012). This suggested there may exist a rich and dynamic heterogeneity of the cell population whereby NC cells alter cell behaviors and gene expression as they travel through different microenvironments to precise targets.

To address this question, we subdivided the pre-otic, cranial NC migratory stream into eight subregions (Fig. 1A-D). Hierarchical clustering of 77 genes revealed distinct differences in gene expression between each of the eight subregions (Fig. 1E,F). Subregions 1-3, 4-5 and 6-8, clustered together (Fig. 1F). When we organized the genes according to their level of linear expression in subregion 1 (migratory front) relative to subregion 8 (newly emerged), we found non-linear transitions of gene expression between the clustered subregions (Fig. 1F). Some of the genes that were highly expressed within the migratory front (subregion 1) were dramatically diminished in expression towards the proximal subregions (Fig. 1F). Likewise, some genes expressed at low levels within the migratory front (Fig. 1F, bottom half of graph) showed a significant increase in expression throughout the more proximal subregions. The lead three subregions (1-3; corresponding to the invasive front and distal portion of the cranial NC cell migratory stream) cluster very closely with one another and are distinct from the newly emerged NC cells (subregions 6-8) (Fig. 1G). Of the profiled genes, 44% (34/77) are differentially expressed between these lead three subregions and the more proximal trailing cell subpopulations (subregions 6-8) (Fig. 1G). These data reveal a widespread regional diversity in gene expression profiles based upon cell position within the cranial NC cell migratory stream.

To determine whether the NC cells within the most advanced portion of the migratory front display a unique molecular profile, we compared the gene profiles of cells from subregion 1 (lead 12.5% of migratory stream) with those of cells immediately proximal in subregions 2-3 (Fig. 1H). We found that the migratory front (subregion 1) shows a molecular profile distinct from that of cells within subregions 2-3, with significant upregulation or downregulation of 18% (14/77) of the genes analyzed (Fig. 1H). Genes significantly upregulated in subregion 1 compared with subregions 2 and 3 included heart and NC derivatives expressed 2 (HAND2), aquaporin 1 (AQP1), BMP and activin membrane-bound inhibitor homology (BAMBI), glypican 3 (GPC3) and matrix metalloproteinase 2 (MMP2) (Fig. 1H). When we compared the molecular profile of lead (subregion 1) NC cells with the rest of the stream (subregions 2-8), we found that similar genes were enriched in the migratory front (HAND2, GPC3, MMP2) (Fig. 1I). Thus, the most advanced migratory front has a unique molecular profile.

To determine whether the NC molecular profiles along the migratory stream could be categorized, we examined the expression of individual genes (Fig. 2). First, we found patterns where genes were more highly expressed in the migratory front and expression decreased proximally (Fig. 2A,E). Second, we observed patterns that peaked mid-stream (Fig. 2B,E). Third, we found that some genes were highly expressed in the newly emerged NC cells in a pattern that diminished towards the migratory front (Fig. 2C,E). Lastly, some genes had uniform expression within all subregions of the stream (Fig. 2D). Thus, the regional differences in NC gene expression varied and displayed at least four characteristic patterns.

To more carefully validate regional differences in gene expression within the NC cell migratory stream, we applied a newly developed technology known as fluorescent HCR (Choi et al., 2010, 2014). The HCR method detects and amplifies mRNA signals within cells and provides a fluorescent readout of expression (Choi et al., 2010, 2014). We applied HCR to simultaneously visualize the expression profiles of two genes highly expressed by cells within the most invasive front (HAND2 and BAMBI) and combined HCR with HNK-1 immunolabeling to distinctly detect NC cells (Fig. 3A,B). Visual observations confirmed that lead NC cells expressed high levels of BAMBI and HAND2 at HH stage 15 (Fig. 3A,C). Quantitative analysis of the HCR signals further confirmed the patterns of BAMBI and HAND2 expression (Fig. 3D). Thus, HCR analysis allowed us to better confirm, both visually and quantitatively, the regional differences in gene expression, including the high expression of specific genes within lead cells at the migratory front.

To test aspects of our hypothesis that are difficult to probe in vivo, we used our extended computational framework (Fig. 4A-H). Specifically, we asked whether the number of lead cells that can detect spatial gradients in chemoattractant is a critical factor for the success of NC cell migration. Should this be the case, then the gene expression patterns detailed in vivo (Figs 1 and 2) would be consistent with this constraint. To test the effect of varying the number of leaders, we first restricted the model simulations to only include non-plastic cell behaviors, such that individual cells that begin migration as leaders (or followers) could not switch from being a leader to a follower (or vice versa). To change the fraction of all migrating cells that are leaders, we varied the time, t_LF, at which new cells entering the domain are prescribed followers instead of leaders (Fig. 4I,J): all cells that entered the migratory domain up to time t_LF were specified as leaders, and cells that entered after that were specified as followers.

Our simulations reveal that both the median distance migrated and the stream density increase with decreasing leader fraction [Fig. 4I; supplementary material Movie 1 compare low (0.051) versus high (0.59) mean leader fraction]. Whereas the furthest distance migrated does not change noticeably, it is the movement of cells away from the entrance to the migratory domain (corresponding to near the dorsal neural tube) that proves crucial for the successful migration pattern (Fig. 4I,J). This prevents jamming near the domain entrance and enables a higher number of cells to distribute more evenly along the migratory pathway (Fig. 4I,J). This trend cannot be extrapolated to zero leaders, in which case the followers move in a random, undirected manner and stay close to the entrance of the domain (supplementary material Movie 2). Thus, our computational model predicts that a small number of lead cells can efficiently guide the migration of the entire NC cell migratory stream.

To test our computational model prediction of a small leader fraction within the migratory front, we performed single-cell analysis and looked for evidence of a stable and consistent set of highly expressed genes within leaders. We isolated an average of eight lead NC cells per embryo by manually removing the tissue containing the first few fluorescently labeled NC cells, dissociating the tissue and then cell sorting via FACS at HH stages 13 and 15 (Fig. 5). We determined FACS isolation to be the most efficient method (as compared with LCM both in vitro and in vivo) for analyzing single NC cells while maintaining the native molecular profile (Morrison et al., 2015). When we profiled 96 genes (the 77 genes described above plus additional genes of interest), we discovered a stable and consistent set of genes expressed at significantly higher levels by a small fraction of the lead NC cells that we define as trailblazers (Fig. 5).

Trailblazers have a high degree of gene expression homogeneity during both phases of migration examined, as can be seen in the profiles of individual cells (Fig. 5A). Principal component analysis (PCA) plots of the genes during each phase of migration analyzed show that over 70% of the genes are stably expressed with little variation (Fig. 5B, blue squares; supplementary material Table S2). Violin plots confirm similar levels of expression for a range of genes, including BAMBI, NOTCH1 and CXCR1 (Fig. 5C). When we focused our attention on the most highly expressed genes, with RT-qPCR Ct values of less than 22, we found that 98% (61/62) of the most highly expressed genes within the migratory front at HH stage 15 were also highly expressed at HH stage 13 (supplementary material Table S2). Many of these genes are consistently expressed by a large percentage of the profiled trailblazers at both phases of migration, including BAMBI, CXCR1, NOTCH1, plakophilin 2 (PKP2) and transcription factor AP-2 alpha (TFAP2A) (Fig. 5B,C; supplementary material Table S2).

Single-cell analysis also shows genes with expression patterns that were dramatically different in the trailblazers when comparing the phases of migration analyzed (Fig. 5B,C). For example, HAND2 is expressed in a small number of lead NC cells (7%) during the first phase of migration (Fig. 5C; supplementary material Table S2). This dramatically increases during the second phase to 59% of the lead cells (Fig. 5C; supplementary material Table S2). By contrast, the percentage of NC cells highly expressing SOX10 dramatically decreases (from 69% to 21%) over time (Fig. 5C; supplementary material Table S2).

To address whether BAMBI and HAND2 were expressed by newly emerging NC cells, we analyzed HCR expression patterns starting at HH stage 9 (supplementary material Fig. S1C). Migratory and premigratory NC cells were visualized by either FoxD3 or HNK-1 expression. As the first pre-otic NC cells delaminate from the dorsal neural tube, we found that cells express BAMBI (supplementary material Fig. S1C, box; HH stage 9). A small number of lead cells expressed high levels of BAMBI at all phases of migration examined (supplementary material Fig. S1C, box; Fig. S3). By contrast, we found that HAND2 was not expressed by newly emerged NC cells (supplementary material Fig. S1C, box; Fig. S3). Therefore, we restricted our focus to genes that are expressed highly by at least 50% of the cells at both phases of migration (supplementary material Table S2, genes in bold). Together, these data confirmed that a stable and consistent set of genes is expressed within the trailblazers, after cells encounter the NC microenvironment.

To determine the molecular signature unique to trailblazers, we compared the stable and consistently expressed genes with those expressed by NC cells within the entire stream at HH stage 15 (Fig. 6A). We measured the similarity of molecular profiles across 318 single NC cells by PCA, violin plots, hierarchical clustering, pairwise correlation and Pearson's correlation (Fig. 6B-D; supplementary material Figs S2 and S3). Our analysis revealed four main results. First, the vast majority of cells from the quartile subregions have a poor correlation with the molecular profile of the trailblazers (Fig. 6B; supplementary material Fig. S2). At the single-cell level, PCA shows incomplete overlap between genes expressed by trailblazers and quartile 1 (Fig. 6B). Second, NC cells within each quartile have high correlations with one another (Fig. 6B; supplementary material Fig. S2). Third, the gene expression profile of the trailblazers is distinct from the gene expression profile of all quartiles, but most similar to quartile 1 (Fig. 6B). Fourth, the average expression of NC cells in the quartiles cluster according to their position within the migratory stream as well as within the local microenvironment depending on the clustering method employed (Fig. 6D; Euclidean, Pearson). Violin plots reveal examples of individual genes that are expressed at higher (BAMBI, CXCR1, PKP2, HAND2) or lower [SOX10, integrin alpha 3 (ITGA3)] levels by trailblazers than cells located more proximal within the migratory stream (Fig. 6C). Measurements of HAND2 and BAMBI HCR expression at higher resolution within the four quartiles confirmed the regional differences in expression (supplementary material Fig. S1A,B). Thus, the unique molecular signature associated with trailblazer cells is distinct from that of cells in the entire first quartile and is not shared by other migrating NC cells within the rest of the stream.

To refine the unique molecular signature of the trailblazers, we examined the genes that were differentially expressed between trailblazers and quartile 1 at HH stage 15 (supplementary material Table S3, genes in bold). This resulted in a list of 17 genes that are highly expressed in trailblazers and differentially expressed between the trailblazers and migrating NC cells in quartile 1. This list of genes did not take into account comparison of the trailblazers with the remainder of the stream. Therefore, a gene expressed highly by cells in quartile 3 and trailblazers would remain in the trailblazer molecular signature. To remove any such genes, we compared the migratory front to single-cell profiles of cells isolated throughout the stream (quartiles 2-4) and found 94% (16/17) of the genes to be in agreement (supplementary material Table S3, genes in bold). Thus, these 16 genes are representative of the unique molecular signature of the trailblazer cells (see Fig. 9).

To determine the changes in the unique molecular signature associated with the trailblazers (see Fig. 9) upon loss of HAND2 or TFAP2A function, we profiled migrating NC cells 24 h after transfection with HAND2 or TFAP2A morpholinos (HH stage 15). We selected TFAP2A and HAND2 based on the fact that they are transcription factors upstream of the trailblazer molecular signature and their expression is enriched within cells in the lead three subregions (Fig. 1G). When HAND2 is knocked down, several trailblazer signature genes are significantly upregulated, including AQP1, CDH11, CDH7, CXCR4 and EPHB1 (supplementary material Table S4, genes in bold). Also, TFAP2A is downregulated in HAND2 morphant embryos, indicating crosstalk between the two transcription factors (supplementary material Table S4). Loss of TFAP2A function results in the upregulation of integrin beta 5 (ITGB5) and NEDD9 and downregulation of CXCR4 and EPHB1 (supplementary material Table S4, genes in bold). Loss of either HAND2 or TFAP2A results in upregulation of HAND2 expression, suggesting either negative regulation or activation of a compensatory pathway. From this, we conclude that loss of either HAND2 or TFAP2A function within migrating NC cells influences the unique molecular signature associated with the trailblazer NC cells.

We hypothesize that genes expressed by trailblazers are crucial to cell direction persistence and stream migration. To test this hypothesis, we knocked down TFAP2A and HAND2 function in NC cells using morpholinos. Loss of HAND2 function results in a significant reduction in the area invaded by treated versus control NC cells (Fig. 7A,F,G). Upon close examination, the fluorescence associated with individual NC cells within HAND2 morphants was more punctate than with control morpholinos (Fig. 7A). This phenotype did not correlate with increased cell death (supplementary material Fig. S4). We also observed some cells at the distal portion of the branchial arches in HAND2 morphant embryos (Fig. 7A-E), suggesting that HAND2 is not required for NC cells to properly distribute along the migratory pathway and colonize the branchial arches.

Loss of TFAP2A function results in dramatic alterations to the NC cell migration pattern (Fig. 7A). First, the distribution of NC cells along the migratory pathway is perturbed when TFAP2A is knocked down at both phases of migration due to a statistically significant drop in cell number in the lead subregion of the migratory stream (Fig. 7B,C). During the first phase of migration (HH stage 13), the distance migrated and area invaded by TFAP2A morpholino-transfected NC cells are not statistically different to those of control embryos (Fig. 7D,F). However, during the second phase of migration, NC cells transfected with the TFAP2A morpholino stopped and failed to migrate the entire length of the migratory pathway (Fig. 7E,G). From these results, we conclude TFAP2A is necessary for NC cell migration into the branchial arches.

To model the overexpression of a trailblazer gene within the trailing subpopulation, we modified the model parameters to convert 50% of trailing cells at random into leaders. This mimics, in silico, the transfection of ∼50% of the NC cells within the trailing subpopulation with overexpression of the transcription factors of the trailblazer signature. When we forced trailing cells in silico to display a lead cell phenotype, we found that cells remain near the entrance of the migratory domain (Fig. 8A,B). This migratory pattern is similar to the model simulations scenario that introduced higher lead cell fractions (Fig. 4I,J). Thus, forcing a lead cell behavior within the trailing subpopulation in silico predicts disruption to NC cell migration that would be observed as cell bunching near the dorsal neural tube exit.

To experimentally test whether gain-of-function of the lead NC cell behavior would alter the migration pattern, we selectively overexpressed TFAP2A or HAND2 in the trailing subpopulation (Fig. 8). In controls, the majority of transfected cells reside in the trailing portion, as predicted (Fig. 8C). When HAND2 is overexpressed in trailing NC cells we found fewer migrating cells (Fig. 8C,E). However, these fewer migrating trailing cells distributed along the migratory pathway (Fig. 8D,E). When TFAP2A is overexpressed in trailing NC cells, we observed some cells bunching near the neural tube, but for the most part cells continued to distribute along the migratory pathway (Fig. 8C-E). These results agreed with our model simulations (Fig. 8A,B) but also supported our previous tissue transplantation experiments: lead cells placed into the trailing stream either stalled or re-initiated their migration to distribute along the migratory pathway, but did not overtake the migratory front (McLennan et al., 2012). Thus, the gain-of-function, model simulations and tissue transplantation experiments all support the hypothesis that distinct lead and trailing cell behaviors and gene expression must be maintained for proper NC cell migration.

In this study, we addressed questions in directed cell migration during vertebrate development using the cranial NC model. Using novel methods to isolate and profile single and small numbers of cells, we analyzed and compared NC cell gene expression profiles during distinct phases of migration. This led us to discover regional differences in gene expression within the NC stream and a consistent (Figs 1 and 3) and stable molecular signature unique to the cells within the most distal portion of the migratory front (Figs 5 and 6), which we termed trailblazers. Gain- and loss-of-function experiments then revealed insights into the roles of genes associated with the trailblazer molecular signature that ensure the proper pattern of NC cell migration. In parallel, we used a hybrid computational model to simulate and predict experimental outcomes.

Only a few lead NC cells with guidance information are required for persistence of direction and stream maintenance. Our previous computational model predicted at least two separate cell subpopulations (leaders/trailers), which was confirmed by RT-qPCR analysis that divided the stream into a 30/70 percentage split (McLennan et al., 2012). What was unclear, and difficult to manipulate experimentally, was whether a change in the number of leaders would affect stream dynamics. This is where our extended computational model proved very useful. In fact, we were able to show that the furthest distance migrated was insensitive to the number of lead cells (Fig. 4). Even a few lead cells could migrate as far as the entire multicellular stream (Fig. 4). Moreover, higher numbers of leaders were less efficient at guiding the entire stream.

Further analysis of lead NC cells provided insights into the molecular characteristics of the migratory front. Our single-cell analysis led to the identification of a unique molecular signature associated with a few lead cells, which we termed trailblazers (Figs 5 and 6). Trailblazers are narrowly confined to the most invasive front. The molecular signature of the trailblazers (16/96 genes) is distinct from the molecular profiles of other migrating cells within the stream (Figs 6 and 9; supplementary material Figs S2 and S3). It is important to note that some of the other 80/96 genes did show changes in expression over time in the trailblazer cells, suggesting that a majority of the molecular profile of trailblazers is influenced by the microenvironment (Figs 5 and 6; supplementary material Figs S2 and S3). Thus, computational model simulations that predicted that only a few cells are required to direct stream migration in the presence of a chemoattractant on a growing domain led to the discovery of a small subpopulation of lead cells narrowly confined to the migratory front with a consistent and stable molecular signature.

Based on our identification of a unique molecular signature associated with trailblazer NC cells, we next asked whether trailblazer gene function is important for NC cell persistence of direction and stream maintenance. Since HAND2 expression was high in lead cranial NC cells, we expected its knockdown would affect the migration pattern. NC cells with HAND2 knockdown did invade less area, but continued to reach distal portions of the second branchial arch (Fig. 7). This might be explained, though, by the differences in HAND2 expression levels at the two distinct phases of migration analyzed. HAND2 expression is high at HH stage 13 but only in 7% of the trailblazer cells. By HH stage 15, HAND2 is highly expressed in 59% of trailblazers. Thus, we conclude that knockdown of HAND2 does not affect early phases of cranial NC cell migration but may play a role in the distribution of cells throughout the branchial arch and in subsequent cell differentiation (Thomas et al., 1998; Howard et al., 1999; Srivastava et al., 1997; Miller et al., 2003; Hendershot et al., 2007; D'Autreaux et al., 2007).

Knockdown of TFAP2A resulted in the failure of NC cells to migrate completely into the branchial arches (Fig. 7). TFAP2A has previously been shown to be expressed by premigratory and migratory cranial NC cells and be important in NC induction, proliferation and differentiation in mice, zebrafish and Xenopus, as well as avian facial tissue growth (Mitchell et al., 1991; Chazaud et al., 1996; Schorle et al., 1996; Shen et al., 1997; Pfisterer et al., 2002; Brewer et al., 2002; Luo et al., 2003; Knight et al., 2003,, 2004; Barrallo-Gimeno et al., 2004; Li and Cornell, 2007; Wang et al., 2011). We have previously shown that when leaders are prevented from migrating into the target site by a physical barrier, trailing NC cells sense the paused leaders, reroute around the barrier and become the new leaders (Kulesa et al., 2005). Here, because NC cells throughout the stream were transfected with the TFAP2A morpholino, no cells were able to take on the role of the trailblazers and migration was hindered (Fig. 7). Together, these data support our hypothesis that only a small number of trailblazer cells, with an expression profile that is influenced by TFAP2A, are necessary to guide the NC stream in a directed manner.

Computational model simulations that tested the gain-of-function of trailblazer genes within the trailing subpopulation predicted alterations to the NC migration pattern (Fig. 8). When similar experimental perturbations were performed in ovo, we observed alterations to the NC migration pattern (Fig. 8). That is, when HAND2 was overexpressed in the trailing NC cells, fewer cells were observed along the migratory pathway (Fig. 8). We initially interpreted this result as a possible NC cell delamination defect. However, after comparing the phenotype with computational model simulations (the first 50 μm of the computer model migratory domain corresponds to NC cell migration from the dorsal midline into the paraxial mesoderm), it is more likely that the majority of trailing NC cells that overexpress HAND2 properly delaminate but then fail to persist in directed migration and instead reside above the neural tube (before our in vivo measurements of the stream begin). Therefore, we hypothesize that HAND2 expression in vivo contributes to the trailblazer phenotype. This is consistent with the computational model simulations that show that leaders within the trailing portion of the stream fail to migrate due to a lack of chemoattractant to follow. Trailing NC cells that look for and follow the stalled leaders also become stalled as a result. Although we do not have direct in vivo evidence that lead cells detect spatial gradients of chemoattractants differently from trailers, our results clearly show that association of lead cell behaviors with genes expressed at the invasive front is valid. By contrast, when TFAP2A was overexpressed in the trailing portion of the migratory stream, no obvious defects were observed (Fig. 8). Since TFAP2A was more broadly expressed in the lead subregions of the NC cell migratory stream and is not restricted to the most distal migratory front (Fig. 1H), this result was not surprising.

Theoretical testing of the effect of the number of lead cells on multicellular stream migration required that we restrict our in silico experiments to the case of non-plastic lead and trailing cell behaviors. We did this even though our previous tissue transplantation experiments showed that chick cranial NC cells alter their gene expression profile depending on the position within a stream. Here, we disabled phenotype switching in our computational model simulations, for otherwise we would have only been able to change the leader fraction transiently before phenotype switching restored it towards the unperturbed case. Alternatively, we could have set up the computational model simulations with varying initial lead cell fractions and, given appropriate microenvironmental parameters, observed what leader fraction emerges naturally from phenotype switching. However, this would require at least a phenomenologically correct implementation of switching, which in turn has to be verified experimentally. This is outside the scope of the current manuscript, but will be addressed in future work.

In summary, we show that the embryonic NC microenvironment regulates the gene expression profile and pattern of cranial NC cell migration in a manner that is dependent on cell position and phase of migration. This regulation was identified by the presence of a unique molecular signature associated with trailblazer NC cells that are narrowly confined to the distal migratory front. These data support the hypothesis that a few lead NC cells interpret complex microenvironmental signals differently than other migrating NC cells within the multicellular stream. To test our hypothesis and the importance of the unique molecular signature, we showed that misexpression of transcription factors TFAP2A and HAND2 may result in significant alterations to the NC cell migration pattern. We speculate, based on experiment and computational model simulations, that only a few trailblazers drive directed migration, and perturbations in the identity and function of these trailblazers lead to migration pattern defects. When a gene associated with the trailblazer signature is overexpressed (experiment) or the trailblazer cell behavior is forced (simulation) within the trailing subpopulation, cells are stranded near the migratory domain entrance. We postulate this is due to a lack of guidance information. These findings were made possible by studying NC cellular and molecular dynamics within the embryonic microenvironment using a closely integrated experimental and theoretical approach.

Fertilized chicken eggs (Centurion Poultry) were incubated at 38°C in a humidified incubator until the desired stages of development (Hamburger and Hamilton, 1951). Plasmid DNA (2.5-5 μg/μl) or fluorescein-tagged morpholinos (0.5 mM) were injected into the neural tube and electroporated at HH stage 9 (McLennan and Kulesa, 2007). We focused the morpholino knockdowns to specific developmental stages of interest with large number n values, and we did not observe any obvious off-target effects. For electroporations of the trailing NC, embryos were incubated until 10-12 somites, at which point plasmid DNA was injected and electroporated with CellTracker CM-DiI (C-7001, Life Technologies).

Tissue was harvested by LCM (Zeiss), pre-amplified using a modified version of the Cells-to-Ct Kit (Ambion) and analyzed by microfluidic RT-qPCR on the BioMark HD (Fluidigm) as previously described (Morrison et al., 2012). The stream lateral to rhombomere 4 (r4) (in cryosections) was divided into eight subregions (the number of subregions selected for our ability to cut reproducibly) and harvested by LCM. Following LCM, RNA from residual cryosections produced RNA integrity numbers of 5.8-6.8 on a Bioanalyzer 2100 (Agilent). Seventy-seven transcripts were pre-amplified from cDNA using 14 cycles as per the Cells-to-Ct protocol (Life Technologies). Pre-amplified cDNAs were diluted with sterile 1×TE and the products analyzed on a BioMark HD at the Children's Hospital Boston, IDDRC Molecular Genetics Core. Median absolute deviation (MAD) was used to eliminate outliers, resulting in three to six biological replicates per sample. Different outliers were automatically removed when comparing subregions 1-8 than when comparing subregions 1-3. A trio of reference genes, selected from six candidates, was used to calculate normalized relative quantities and differential expression in qBASEplus software (Biogazelle). Hierarchical clustering by Pearson's dissimilarity was performed using Partek Genomics Suite.

For the quartile analysis, four regions of interest (the number of subregions selected for our ability to manually dissect reproducibly) from the cranial r4 migratory stream, electroporated with Gap43-yellow fluorescent protein (YFP), were chilled in PBS containing 0.1% DEPC and mechanically and chemically dissociated. For the trailblazer analysis, the tips of the arches, containing an average of eight electroporated NC cells, were manually removed from the rest of the embryo and then mechanically and chemically dissociated. Single and healthy YFP⁺ NC cells were isolated by FACs into Cells-to-Ct lysis solution (Ambion) with 1:100 DNase I and incubated at room temperature for 20 min. Bioanalyzer 2100 RNA integrity numbers determined from tissue remaining after LCM measured >6.8 for neural tube cultures and >8.9 for cryosectioned tissue. cDNA was synthesized directly from entire lysates using the High Capacity cDNA Kit (Life Technologies). An 18-cycle, Cell-to-Ct pre-amplification protocol was employed to selectively amplify 96 transcripts with TaqMan Gene Expression Assays (Life Technologies). Pre-amplification products were diluted with sterile 1×TE before being run on the BioMark HD. Data were analyzed using the Singular2 software package (Fluidigm) in R. Limit of detection (LoD) was set at 28 for the single-cell isolation method comparison and 26 for the trailblazer and quartile single-cell profiling. We calculated normalized relative quantities and differential expression in qBASEplus software and determined hierarchical clustering and intensity plots using Partek Genomics Suite.

Transcripts were visualized in whole embryos and tissue sections by HCR. Chick embryos were incubated to specific developmental stages, tissue rapidly collected in chilled PBS containing 0.1% DEPC and fixed in fresh 4% paraformaldehyde at ambient temperature for 1 h. Embryos were then dehydrated and rehydrated before HCR was performed according to the manufacturer's instructions (Molecular Instruments). Tissue sections were cut on a vibratome (Leica VT1000S). In some cases, immunohistochemistry was performed on the tissue following HCR for NC cell specificity. All samples were imaged by confocal microscopy (Zeiss, LSM 780) and the fluorescence signal was quantified by application of the polyline kymograph analysis (Jay Unruh, Stowers Institute for Medical Research) in Image J (FIJI).

To test the logical conclusions of our mechanistic hypotheses, we used a hybrid computational model (McLennan et al., 2012) with individually represented cells and a continuous chemoattractant concentration (Fig. 4A-D). In the model, cells undertake a two-dimensional off-lattice random walk on a growing rectangular domain that represents the NC cell migratory environment (Fig. 4A-D). New cells enter the domain from one end throughout the simulation (Fig. 4D). The chemoattractant concentration is modeled by a reaction-diffusion equation, with the cells acting as sinks, representing internalization of the chemoattractant. The key model components are illustrated in Fig. 4 and the model parameters are listed in supplementary material Table S1; further information is provided in supplementary material Appendix S1: supplementary model information.

We extended our original computational model (McLennan et al., 2012) to incorporate new experimental findings and performed simulations in parallel with the experiments presented in this paper (Figs 4 and 8). Key changes include: (1) a wider stream of cells that allows for greater cell numbers and more adequate representation of multicellular stream migration; (2) limiting the sensing accuracy of cells (for chemotaxis); based on Berg and Purcell (1977) we derived (order of magnitude) bounds of how small a local chemoattractant gradient can be relative to the bulk concentration before cells cannot sense it; (3) assuming that when cells are in contact, the length of the filopodium giving rise to the contact is not fixed; instead, a range of intercellular distances is allowed, with a maximum beyond which cells lose contact and cease to communicate directional information; these distances are based on empirical data (supplementary material Fig. S5) and this change improves stream cohesion and reduces stream breakup in model migration.

P.M.K. thanks the Stowers Institute for Medical Research for their kind generosity. L.J.S. thanks Louise Dyson for helpful discussions and access to the original code and gratefully acknowledges the UK Engineering and Physical Sciences Research Council (EPSRC) for funding through a studentship at the Life Science Interface programme of the University of Oxford's Doctoral Training Centre. R.M. and J.A.M. thank the Histology Core at the Stowers Institute for Medical Research for assistance with cryosectioning.