During the development of the central nervous system, all chordates deploy a suite of genes to separate neural and non-neural ectoderm, including homologs of Tfap2, Zic, Msx1/2, Pax3/7, and Dlx3/5 [2,5,26,27]. In vertebrates, these genes activate a suite of neural crest marker genes, including Snail, Id, Tfap2, Twist, FoxD3, and SoxE within neural plate border (NPB) cells [2,5]. These NCC specifier genes consist of transcription factors that are deployed at different stages of NCC development. In invertebrate chordates, homologs of most NCC specifiers, with the exception of snails, are absent from the neural border and expressed in other germ layers or ectodermal domains, suggesting they were co-opted to neural border cells in vertebrates [28,29]. As neurulation proceeds, the vertebrate NPB cells expressing a combination of NCC specifier genes become true neural crests.

There are many subnetworks (also known as subcircuits) within the NC GRN that regulate different aspects of the NCC phenotype during different phases of their development including pluripotency, migratory capabilities, and differentiation capacity. These subcircuits use many of the same transcription factors and signaling molecules as the core NC GRN at different times to attenuate the activation and suppression of downstream target genes. NCC marker genes regulate downstream genes that initiate their delamination from the ectoderm during the epithelial-to-mesenchymal transition (EMT), which is marked by the expression of Snail, FoxD3, Twist, Lmo4, and Zeb2 [27]. Regulation of the EMT subcircuit leads to the dynamic expression of cadherin proteins regulated by transcription factors, such as Snail and Zeb2 [30]. Some vertebrates have slight differences in cadherin protein activity during the EMT, but many express type I cadherins prior to delamination and switch to type II cadherins prior to and during migration, orchestrating changes in the cytoskeleton that allow the proper NC movement [27,30,31]. These genes, along with Ets1, c-Myc, Tfap2, Id, and SoxE, are expressed after NCCs separate from NPB and remain active during their migration throughout the body [5,27]. It is important to note there are differences between the gnathostome and cyclostome NC GRN, as pre-and post-migratory 22NCCs express various NC-specific components at different times and with different expression boundaries, but gene swap experiments show that these homologous proteins are functionally similar between lineages [2,9,27,28]. The evolutionary significance of the heterochronic expression of these components of the NC GRN between vertebrates remains unknown.

The initial positioning of NCCs along the anterior–posterior axis affects their ultimate post-migratory destination and fate. In gnathostomes, NCCs form four major subpopulations: cranial, vagal, cardiac, and trunk/sacral, all of which give rise to different structures throughout the body [28,32,33]. In the other major lineage of extant vertebrates, the jawless cyclostomes, the precise boundaries, and derivatives of non-cranial NCC subpopulations are less clear. Cranial NCCs migrate into the pharyngeal arches (PAs) in a conserved pattern across vertebrates, where they receive various signals dependent on their anterior–posterior position that determine their skeletal fate after migration [16,28,34]. At the NPB, amniote cranial NCCs (CNCCs) are marked by the co-expression of Bm3, Lhx5, and Dmbx1 that activate SoxE, Tfap2, and Est1, which is maintained throughout the migration to the PAs [10]. Recent work in mice revealed that CNCCs reactivate pluripotency marker Oct4 after delamination and reset their positional information and become transcriptionally equivalent prior to migration into the PAs [35]. Whether this is a conserved feature of the NC GRN, or unique to amniotes, is unclear.

Skates, a representative of jawed cartilaginous fish, lack the early CNCC markers but deploy SoxE, Tfap2, and Ets1 prior to and during migration, while zebrafish share a majority of their CNCC specification and migration GRN with amniotes, with the exception of bm3 [10]. This evidence highlights a stepwise acquisition of the CNCC GRN in gnathostomes as well as a high conservation of their expression patterns. Lamprey lack the vagal stream of NCCs, with trunk NC-derived Schwan precursor cells giving rise to their enteric nervous system, which stems from vagal NCCs in gnathostomes [28,36,37]. Generally, lamprey CNCCs exhibit a GRN more comparable to a trunk neural crest than those in gnathostomes, with migratory CNCCs being marked with SoxE and Tfap2 [9,10,38]. They also express some genes orthologous to amniote early cranial specifiers, Lhx5 and Dmbx1, later in the PAs [10]. Additionally, the migration patterns of CNCCs in lamprey are less restricted, with the cells destined for the posterior pharyngeal arches migrating initially as a sheet rather than distinct streams [34]. This may be a result of heterochronies between the NCC migratory GRN of jawed and jawless vertebrates [34].

The vertebrate head skeleton can be divided into the neurocranium, which encases the brain and can have a mesodermal component, and the viscerocranium, which develops from NCCs in the pharyngeal arches [39]. NCCs populate the PAs during head development and give rise to unique structures in each arch. While there are many cranial NCC derivatives that are shared by all vertebrates, hagfish and lamprey lack some NCC derivatives, such as jaws and cranial sympathetic ganglia [34,37]. Due to the inaccessibility of hagfish embryos, the vast majority of what is known about NCC development in agnathans is from studies of lamprey embryos.

Intercellular signals that are received by CNCCs as they migrate into the head are generally conserved between jawed and jawless vertebrates [16]. Within the PAs, both cyclostomes and gnathostomes develop NC-derived skeletal structures. However, the cartilage composition differs in the pharyngeal skeleton, and lampreys totally lack bone [40]. Fibroblast growth factor (Fgf) expression is activated by retinoic acid (RA) signaling within the pharyngeal pouch endoderm and is crucial for cartilage development in all vertebrates. FGF receptors on CNCCs receive that signal and activate a downstream cartilage regulatory module containing Dlx, SoxE, Twist, and Ets [16,40]. A similar GRN is deployed in the cartilage of amphioxus oral cirri, implying an ancient cartilage regulatory subcircuit may have been co-opted into CNCCs during vertebrate evolution [4]. Chondrocytes differentiate and give rise to cellular and “soft” cartilages that permit the stability and flexibility of the viscerocranium. In addition to the cartilage differentiation genes mentioned above, gnathostomes require Barx and Runx for proper facial cartilage/bone development, while lamprey only deploy those genes in the branchial vasculature of the PAs [41].

All vertebrates deploy a conserved subset of transcription factors within CNCCs, filling the PAs that give rise to different structures and skeletal fates [16]. Within the head, CNCCs give rise to both dermal and endochondral bone [42] Dermal bones, such as mandibular bones forming from the first PA, ossify directly from mesenchymal cell matrices. Endochondral bone, such as the CNC-derived parietal bone of the skull, requires a cartilage intermediate prior to ossification [42]. Just before FGF signaling activates chondrogenesis, other signaling pathways including BMP, endothelin (Edn), and Notch initiate the patterning of the CNC and resulting head [16,23]. Twist, Ets, Id, Alx, and SoxE orthologs are expressed in all vertebrate skeletogenic CNCCs during migration and the population of the PAs [41]. Skeletal elements differentiate through the polarized and combinatorial expression of Alx, Hand, Msx, and Prrx around a core of nested Dlx expression along the dorsal–ventral (DV) axis that corresponds to unique structures in the head [16,41]. These genes cooperate with Hox genes expressed along the anterior–posterior (AP) axis of the PAs to give rise to transcriptionally distinct CNCC populations in each arch that give rise to skeletal elements of various shapes and properties [16,43]. Dlx and Hox genes are thought to pattern the DV and AP axes of the vertebrate head in a highly conserved, code-like pattern [16,44].

The term “ohnolog” was coined to distinguish gene duplicates that were products of vertebrate-specific WGDs from paralogs that are exclusive to a single lineage and orthologs shared by multiple lineages [12]. All genes in the NC GRN appear to have been duplicated during the vertebrate WGDs, with the resulting ohnologs being differentially retained across lineages (Figure 2). However, virtually all NC GRN ohnologs function at some point during neural crest development. It would be highly unlikely for each duplicated member of these gene families to be independently co-opted into the same GRN. Thus, the fact that virtually all NC GRN ohnologs are active at some point in the NC GRN, or its various subcircuits, serves as robust evidence that the co-option of these genes to the neural border, and the assembly of the NC GRN, occurred before the WGDs.

In this context, the differential expression of NC GRN ohnologs likely reflects the temporal and spatial subfunctionalization of NC GRN components after the WGDs (depicted in Figure 3). The functional consequences of this extensive subfunctionalization of NC GRN ohnologs remain unclear. However, in general terms, these duplications appear to have substantially increased the overall complexity of the modern NC GRN by creating temporally and spatially restricted subcircuits [2,3,7,10,11]. The dedication of individual NC GRN ohnologs to particular phases of the NC GRN, or particular NC subpopulations, may have increased the modularity of the GRN, allowing different portions to evolve without interfering with its other functions. Consider an NC GRN gene that was ancestrally involved in the initial activation of the NC GRN (i.e., NCC specification) and later during the differentiation of two NCC derivations, derivative 1 and derivative 2. Suppose this gene was duplicated during the WGDs, and its three retained ohnologs, A, B, and C, became temporally and spatially subfunctionalized: A being expressed early and dedicated to NCC induction, and B and C expressed later and dedicated to derivative 1 and 2, respectively. If a mutation affected the expression or function of the B ohnolog, causing an adaptive change in derivative 1, it would have minimal effect on NCC specification or the formation of derivative 2. Conceivably, this subfunctionalization of ohnologs could then lead to biochemical or regulatory neofunctionalization, as they evolved new biochemical properties or expression domains.

Transcription factors of the SoxE family participate in multiple steps of NCC development. While there is one SoxE in invertebrate chordates, SoxE is duplicated in vertebrates with three each in non-teleost gnathostomes (Sox8, Sox9, and Sox10) and lamprey (SoxE1, SoxE2, and SoxE3) and five in teleosts [4,45,69]. SoxE3 is in lamprey and gnathostome is in Sox9, and both have a conserved role in regulating cartilage morphogenesis. SoxE2 is comparable to Sox10, as both have similar melanogenic and glial-inducing roles [45,69]. All vertebrate SoxE paralogs are expressed during induction, migration, and differentiation of NCCs but with heterochronic differences across lineages [22,45,69]. The most parsimonious explanation for this is that SoxE was co-opted into the NC GRN before the first WGD at the base of the vertebrates. 1R ohnologs were then duplicated in gnathostomes during the 2R WGD, while one or both 1R ohnologs were independently duplicated in the cyclostome lineage. Based on the similar expression of SoxE3/Sox9 in CNCCs and SoxE2/Sox10 in the trunk, the subfunctionalization of R1 SoxE duplicates likely occurred before the gnathostome/agnathan split.

Interestingly, AmphiSoxE is expressed in the fibrillar cartilage of the oral cirri in amphioxus, pointing to an ancestral role for SoxE in cartilage development [4]. AmphiSoxE is capable of activating downstream genes in NC GRN modules, such as neurogenic Phox2 and melanogenic Mitf, within transgenic zebrafish, expressing AmphiSoxE and other NC specifiers in transgenic chicks [22,45]. This indicates an ancestral function in differentiation and NCC induction in the ancestral SoxE. However, AmphiSoxE regulatory elements are not capable of driving reporter expression within NCCs in transgenic zebrafish, indicating that a cis change permitted SoxE function in neural crest GRN [4]. Additionally, transgenic mice expressing the Drosophila ortholog of SoxE (Sox100B) in place of Sox10 developed many NC derivatives, including the peripheral nervous system (PNS) and melanocytes, showing that Sox10 has retained a highly conserved role of inductive capabilities [70]. However, when Sox8 is expressed in place of Sox10, transgenic mice fail to develop melanocytes, indicating that gnathostome Sox8/9/10 proteins have diverged in functionality [3,71]. Taken together, these data suggest that after an ancestral SoxE was co-opted into the NC GRN by a cis-regulatory mutation, it underwent duplication and both temporal and spatial subfunctionalization of its ancestral expression pattern. This was followed by some divergence of the SoxE protein function, reflecting either subfunctionalization or neofunctionalization of SoxE ohnologs.

The Dlx transcription factor family, which underwent an ancient tandem duplication in the last common ancestor of vertebrates and urochordates, is homologous to the distal-less gene (Dll) in insects [72]. This tandem gene pair was duplicated in gnathostomes into pairs of multiple copies across different chromosomes in a pattern that agrees with 2R duplication events [72,73]. As with other gene families, Dlx orthology is not perfect between cyclostomes and gnathostomes, and teleosts have more Dlx paralogs than other gnathostomes. Each tandem duplicate has its own orthology group, annotated Dlx1/4/6 and Dlx2/3/5, in gnathostomes. While in lamprey, Dlx orthologs are annotated as DlxD/E/F and DlxA/B/C, respectively [16,73]. The last common vertebrate ancestor is hypothesized to have duplicated its tandem Dlx clusters during the 1R event with lineage-specific duplications followed by differential retention, resulting in all vertebrates having six Dlx genes [16,72,73].

The Dlx family has a highly conserved role in patterning, especially within the vertebrate head where it exhibits a nested expression pattern along the DV axis of the Pas; however, the precise borders of these expression patterns vary between jawed and jawless vertebrates [16,74]. The sea lamprey exhibits nested DV expression of DlxA-D paralogs and regionalized Hand1 and Msx transcription, reminiscent of the Dlx-Hand-Msx code in gnathostomes [75]. This nesting is less obvious in the Japanese lamprey with the exception of DlxF, which is differentially expressed in the nasal region [43]. As in gnathostomes, the combinatorial differences in Dlx, Hand, and Msx expression domains along the DV axis also correlate with different cartilage types within the lamprey head [75].

As with other NC GRN genes, all Dlx paralogs are expressed in NCC at some point during NCC development, strongly supporting the idea that the ancestral Dlx gene pair was co-opted into the NC GRN before the WGDs [16]. The distinct expression patterns of Dlx ohnologs further suggest the subfunctionalization of an ancestral, pan-NCC expression pattern. The fact that different combinations of Dlx genes mark different portions of the head skeleton is consistent with the participation of different Dlx ohnologs in separate regulatory subcircuits in the NC GRN. The degree to which the biochemical function of the transcription factors encoded Dlx genes was affected by WGDs is unclear, as there is limited evidence that Dlx proteins have qualitatively different DNA binding properties. However, the divergent expression of Dlx ohnologs supports the idea that the WGDs allowed the regulatory landscape of the Dlx family to diverge, resulting in the nested expression of the “Dlx Code” [40]. The complex expression of Dlx ohnologs within the vertebrate head may be consistent with the idea that the patterning of NCC derivative became more elaborate following the vertebrate-specific expansion of this gene family.

The Hox gene family has been presented as evidence for the 2R hypothesis, one reason being that amphioxus possesses a single, syntenic cluster of tandem Hox homologs, while gnathostomes have four homologous clusters across separate chromosomes [76,77,78]. Hox clusters in teleosts are also evidence for teleost-specific genome duplication (3R) and are distributed across seven chromosomes [72]. Alternatively, some have used data from phylogenetic analyses of human Hox genes as evidence against WGDs, and propose that this gene family was expanded by small-scale events that occurred in vertebrates [79,80,81]. There are ancient paralogical groups of Hox genes that are annotated numerically within clusters (i.e., Hox1, Hox2, Hox3…). These groups are cis-tandem duplicates that occur throughout metazoans, and Hox ohnologs are specified by letters (i.e., Hoxa1, Hoxb1, Hoxc1…) and occupy separate chromosomes in vertebrates. These transcription factor genes are expressed in a conserved, colinear manner across bilaterians along the anterior–posterior axis and have been shown to induce homeotic transformations if mutated [82,83,84,85].

Within the pharyngeal arches of all vertebrates, the first PA (PA1) is Hox-negative, PA2 is marked by Hox2 paralogs, and PA3, along with the posterior arches, are marked by Hox3 paralogs [43,84,86,87]. The ectopic expression of Hox cis-paralogs within the PAs can lead to homeotic transformations, for instance when Hoxa2 is expressed in the first PA of mice, which is Hox-negative, embryos lose the mandibular structures or take on a second arch identity [88]. Additionally, certain levels of ectopic Hoxa2 expression can alter skeletal element identity in mice, highlighting the possibility that the regulation of Hox gene expression could have more of an influence on skeletal identity than protein function [88]. Knock-downs of hox2 ohnologs in zebrafish have shown that these genes are partially redundant, as removing the function of one ohnolog often results in a less dramatic phenotypic effect than targeting the whole group [85,89]. This highlights the functional differences between Hox tandem duplicates, yet the degree to which Hox ohnologs differ regarding DNA binding specificity is still unclear.

In cyclostomes, the Hox gene orthology is disparate from gnathostomes, providing another line of evidence for cyclostome-specific duplications subsequent to their divergence from jawed vertebrates and the possibility of segmental gene duplications in this family [44,65]. Differences in the Hox cluster number and organization between different species of lamprey and hagfish add difficulty to solving the phylogeny of Hox clusters [44,84]. Although Hox orthologues are hard to assign between jawed and jawless vertebrates, there is evidence that cis-regulatory elements (CREs) upstream of Hoxa1/Hoxb1 in gnathostomes and hoxα1 in lamprey are capable of driving similar expression in NC [44]. Similar CRE activity provides minor support for orthology assignment as well as the nested colinear expression pattern of tandem duplicates in the pharyngeal arches [44,84]. Hox ohnologs have been shown to have both redundant and diverged traits in NC patterning, making them a complex gene family for understanding their evolution after duplication events. Recent work in zebrafish has shown that within the PAs, certain levels of ectopic hoxa2 expression can alter skeletal element identity, highlighting the possibility that the regulation of Hox gene expression could have more of an influence on skeletal identity than protein function [88].

Cranial neural crest cells have endothelin receptors that interpret signals from the surrounding tissues to activate downstream genes, including NCC specifier genes and Dlx [16,23,41]. The endothelin receptor genes (Ednra/Ednrb) are both critical components of the NC GRN, suggesting they are likely “1R” duplicates of a single ancestral Ednr present in the last common vertebrate ancestor [23,90]. The 2R WGD in gnathostomes was presumably followed by the loss of two Ednrs, leaving only one Ednra and one Ednrb, while the teleost 3R resulted in four Ednr genes in this group. Following the gnathostome–cyclostome split, Ednrs evolved somewhat divergent roles in CNCC differentiation [23,91]. Within gnathostomes, Ednra has a key role in skeletogenesis and vascular development, while Ednrb contributes to melanocyte and peripheral nervous system development, and Ednrb mutants have no craniofacial defects [92]. In contrast, lamprey Ednra and Ednrb both are required for proper PA formation, and those lacking either Ednra or Ednrb fail to form branchial arch cartilages properly [23,92]. Recently, it was discovered that skate embryos also express Ednrb in the ventral and intermediate NCCs in all pharyngeal arches in a way that resembles lamprey expression [93]. Together, these data imply that after the R2 WGD, Ednr ohnolog expression diverged in chondrichthyans and bony fish, with the developmental roles of Ednra and Ednrb becoming more specialized in the latter [23,94].

Endothelin receptors on the surface of NCCs activate genes that contribute to the skeletal phenotype within each arch. Within gnathostomes, Ednra activates Hand and Dlx5/6 genes within the ventral portion of the PAs in zebrafish, mice, and frogs, where they contribute to the lower jaw and joint formation [23,92,95]. Lamprey exhibit similar Dlx expression within the ventral Pas; however, hand expression is not regulated by Ednra [23]. Duplicated genes can subfunctionalize or neofunctionalize following duplications; Ednra and Ednrb appear to have diverged in terms of expression pattern and downstream gene targets [23,41]. It is possible that a stem or independent duplication provided the flexibility necessary for agnathans and gnathostomes to alter endothelin signaling, supporting the morphological specialization of their head skeletons. Gene swap experiments between cyclostome and gnathostome Ednra/Ednrb orthologues will highlight the divergent aspects of these paralogous signaling pathways. Ednra/Ednrb divergence following duplication is a potential link between head specialization and “2nd R”, but further research will have to be conducted to deem it as a necessary event for that divergence.

Alx genes are part of NC GRN subcircuits that confer cellular identity in NCC-derived skeletal elements, contributing to their unique shapes and different components [16,41,66]. This gene family has a conserved skeletal role across deuterostomes and is deployed in the skeletal GRNs of echinoderms and vertebrates [96]. Amphioxus possesses two copies of Alx (Bf-alx1, Bf-alx2)m whose loci are near each other and share highly similar intron/exon structures, indicating a lineage-specific, tandem duplication [96]. Non-teleost gnathostomes have three Alx homologs, Alx1/Cart1, Alx3, and Alx4, and teleosts possess additional paralogs of Alx4 (alx4a, alx4b) that are generated during the vertebrate and teleost WGDs, demonstrating that these paralogs have been differentially retained [66,96]. Lamprey possesses two Alx homologs that were likely generated during the first vertebrate WGD and orthologous to gnathostomes Alx1 and Alx4 [41]. Only the expression of lamprey Alx4 has been reported, and like gnathostome Alx ohnologs, it is expressed in skeletogenic CNCCs [41].

The expression of different Alx paralogs corresponds to different cellular cartilage shapes within the vertebrate head [16,66,67,75,96]. In lamprey, Alx4 expression is largely coincident with stiffer, more rigid cartilage phenotypes, indicating an ancestral role in the differentiation of skeletal tissue subtypes [41]. In zebrafish, combinations of Alx paralogs label NCC-derived chondrocytes with distinct cellular phenotypes. Detailed functional analyses suggest that this “Alx code” contributes to the intricate shapes of the gnathostome head skeleton [66]. In mice, Alx1 is activated earlier and is concentrated in the midline regions of the frontal nasal prominence (FNP), while Alx3 and Alx4 are largely limited to the lateral tissues of the FNP [67,68]. Experiments in mice have also shown that Alx paralogs have diverged in function, with homozygous deletions of Alx1 or Alx3 resulting in different craniofacial malformations, while Alx4 loss-of-function results in minimal defects [67]. Deciphering the evolution of the gnathostome “Alx code” would be aided by gene swap experiments within and between species or the overexpression experiments of various Alx1/3/4 paralogs. Furthermore, analyzing the open chromatin regions surrounding paralogs within various lineages will help determine whether regulatory changes are at the base of their divergence or if there are protein sequences that have diverged biochemically.