The developmental hourglass model is applicable to the spinal cord based on single‐cell transcriptomes and non‐conserved cis‐regulatory elements

Abstract The developmental hourglass model predicts that embryonic morphology is most conserved at the mid‐embryonic stage and diverges at the early and late stages. To date, this model has been verified by examining the anatomical features or gene expression profiles at the whole embryonic level. Here, by data mining approach utilizing multiple genomic and transcriptomic datasets from different species in combination, and by experimental validation, we demonstrate that the hourglass model is also applicable to a reduced element, the spinal cord. In the middle of spinal cord development, dorsoventrally arrayed neuronal progenitor domains are established, which are conserved among vertebrates. By comparing the publicly available single‐cell transcriptome datasets of mice and zebrafish, we found that ventral subpopulations of post‐mitotic spinal neurons display divergent molecular profiles. We also detected the non‐conservation of cis‐regulatory elements located around the progenitor fate determinants, indicating that the cis‐regulatory elements contributing to the progenitor specification are evolvable. These results demonstrate that, despite the conservation of the progenitor domains, the processes before and after the progenitor domain specification diverged. This study will be helpful to understand the molecular basis of the developmental hourglass model.


| INTRODUC TI ON
The vertebrate neural tube can generate a diverse array of neurons in a precisely controlled and reproducible manner. During neurogenesis and neuronal differentiation, each neuron is assigned a distinct feature, such as neurotransmitter phenotype, axonal projection pathway, and cell body localization (Lai et al., 2016;Lu et al., 2015;. The initial phase of this process is the fate specification of progenitor cells, whereby molecularly defined progenitor domains are established with sharp boundaries along the dorsoventral axis of the neural tube (Briscoe et al., 2000). Each progenitor identity is specified through a highly complex gene regulatory network (GRN), which consists of a graded Shh signaling activity localized ventrally, Bmp | 373 and Wnt signaling molecules expressed dorsally, pan-neural transcriptional activator Sox1-3, and a number of domain-specific transcription factors (TFs) that function as repressors (Andrews et al., 2019;Balaskas et al., 2012;Delás & Briscoe, 2020;Kutejova et al., 2016;Nishi et al., 2015;Oosterveen et al., 2012;Peterson et al., 2012;Zagorski et al., 2017). As the output of the GRN, 11 progenitor domains, termed dp1-6 in the dorsal half, and p0, p1, p2, pMN, and p3 in the ventral half, are established (Lai et al., 2016;Lu et al., 2015;. Subsequently, multiple neuronal subtypes are generated from a single progenitor domain. For example, V2a, V2b, and V2c interneurons (INs) are differentiated from the p2 domain (Del Barrio et al., 2007;Li et al., 2005;Panayi et al., 2010;Peng et al., 2007), and dI1i and dI1c INs are differentiated from the dp1 domain (Ding et al., 2012;Wilson et al., 2008). Each subtype of post-mitotic neurons proceeds to the maturation process, such as migration, axonal projection, and circuit formation. The knowledge described above is mainly derived from chick and mouse studies, while a number of studies using zebrafish have demonstrated that the spatial architecture of the progenitor domains in the neural tube is largely conserved among vertebrates Gribble et al., 2007;Guner & Karlstrom, 2007;Lewis et al., 2005;Park et al., 2002;Schäfer et al., 2005).
However, there are overt differences between amniotes and teleosts with regard to post-mitotic neuronal maturation. For example, V2a INs, which are defined by the expression of Vsx2 (Chx10), contribute to the locomotor rhythm generation in zebrafish (Eklöf-Ljunggren et al., 2012), whereas V2a INs in mice play a role in the left-right alternation of limbs by providing excitatory input to the commissural INs (Crone et al., 2008(Crone et al., , 2009). This implies that, during post-mitotic differentiation in mice and zebrafish, V2a INs are assigned different properties or are placed in different positions within the spinal locomotor circuit (Kiehn, 2016). Another example is Robo3, an axon guidance receptor that is essential for commissural axons to cross the midline (Friocourt & Chédotal, 2017;Marillat et al., 2004;Sabatier et al., 2004). In the spinal cord of amniotes, Robo3 is expressed in V0, V1, and V3 INs in the ventral spinal cord, which encompass commissural INs (Friocourt et al., 2019;Tulloch et al., 2019). However, in zebrafish, although the double labeling of robo3 and neuronal subtype markers has not been conducted, robo3 expression is observed in the region encompassing motor neurons (MNs; Challa et al., 2005). If zebrafish MNs express bona fide robo3, the role of Robo3 has diverged between amniotes and zebrafish, as the MNs of amniotes never express Robo3 (Friocourt et al., 2019;Tulloch et al., 2019).
In addition to post-mitotic differentiation, the process before the establishment of progenitor domains also differs between amniotes and zebrafish. Shh is an essential gene for the specification of the ventral neural tube and is expressed in the notochord and floor plate in amniotes (Echelard et al., 1993;Riddle et al., 1993;Roelink et al., 1994). In contrast, three Shh-related genes, shha (sonic hedgehog), shhb (tiggy winkle hedgehog), and ihhb (echidna hedgehog), are expressed in the notochord and/or floor plate in zebrafish, conferring a functional redundancy of hedgehog (Hh) signaling (Currie & Ingham, 1996;Ekker et al., 1995;Krauss et al., 1993). The Hh signaling pathway culminates in Gli family TFs, which function as transcriptional activators or repressors, depending on the Hh signaling activity (Briscoe & Thérond, 2013).
Amniotes possess three Gli genes (Gli1, Gli2, and Gli3), while teleosts possess four Gli genes (gli1, gli2a, gli2b, and gli3). Lossof-function experiments of Gli genes led to clearly different phenotypes between mice and zebrafish. For example, Gli1 knockout mice are viable and appear normal, showing that Gli1 is not essential for embryogenesis in mice (Park et al., 2000). In contrast, gli1-mutant zebrafish displayed severe cranial MN deficiency and reduced Hh-target genes, such as ptch1 and nkx2.2a, and died at the larval stage (Chandrasekhar et al., 1999;Karlstrom et al., 1996Karlstrom et al., , 2003. In Gli2 knockout mice, the floor plate has not been specified, and MNs aberrantly occupy the ventral-most domain in the spinal cord, although MN itself is differentiated (Ding et al., 1998;Matise et al., 1998). In zebrafish, gli2a is completely dispensable for normal embryogenesis and growth to adulthood (Karlstrom et al., 2003;Vanderlaan et al., 2005;Wang et al., 2013), and gli2b knockdown causes a marked reduction of MNs in the spinal cord (Ke et al., 2008). Knockdown of gli3 in zebrafish leads to a reduction of MNs (Vanderlaan et al., 2005); however, in Gli3 mutant mice, patterning defects of the floor plate and MNs have not been observed (Persson et al., 2002).
To summarize the conservation and divergence of spinal cord development between amniotes and teleosts, the dorsoventral arrangement of the progenitor domains is well conserved, while the processes before and after the progenitor domain specification have diverged. This pattern of developmental divergence is highly consistent with the developmental hourglass model, which argues that embryonic morphology at the early and late developmental stages is divergent and that at the mid-embryonic stage is conserved (Duboule, 1994;Hu et al., 2017;Irie & Kuratani, 2014).
However, the process of spinal cord development has never been investigated from the perspective of this model. Herein, we provided evidence that the developmental hourglass model is applicable to spinal cord development. We examined the publicly available single-cell transcriptome data from mice and zebrafish (Delile et al., 2019;Farnsworth et al., 2020), and provided other examples of diverse differentiation of post-mitotic neurons. We also examined the transcriptional regulatory elements in the neural tube patterning genes based on chromatin immunoprecipitation sequencing (ChIP-seq) and assay for transposase-accessible chromatin using sequencing (ATAC-seq) data, and sequence conservation, suggesting that cis-regulatory elements contributing to the progenitor domain specification had undergone turnover (nucleotide changes) during vertebrate evolution. Based on our findings and the robust development of neuronal progenitor specification (Balaskas et al., 2012;Delás & Briscoe, 2020;Exelby et al., 2021;Xiong et al., 2013;Zagorski et al., 2017), we propose that the progenitor domain configuration in the neural tube is less evolvable owing to its canalization (Waddington, 1942(Waddington, , 1957.

| Experimental animals
Fertilized chicken eggs were obtained from Takeuchi farm (Nara, Japan) and incubated at 38℃ in a humidified incubator. Embryos were staged according to the Hamburger-Hamilton (HH) stage series (Hamburger & Hamilton, 1951

| ATAC-seq data analysis
ATAC-seq data (Metzis et al., 2018) were obtained from ArrayExpress (E-MTAB-6337). Data on the neural progenitor cells with spinal cord identity were utilized from the dataset. Read quality control was done as ChIP-seq. Read mapping was performed using Bowtie2 with the following parameters: -X 2000 --sensitive-local. PCR duplicates were marked using Picard (version 2.9.2). Reads with low mapping quality (MAPQ < 30) were removed by samtools. Peak calls were performed using MACS2 with the following parameters: -f BAMPE -g mm -q 0.05 --nomodel --keep-dup auto -B.

| Comparison of genomic sequences by VISTA
For interspecies comparisons of genomic sequences, partial genomic sequences were obtained from the UCSC genome browser (https:// genome.ucsc.edu) and Ensembl (https://www.ensem bl.org/index.html), and VISTA (Frazer et al., 2004) was used to align and visualize the results.

| Vector construction
For RNA probe preparation, partial fragments of chick ROBO3, chick SIM1, chick IRX3, chick DBX2, and chick GSX1 were amplified by PCR from chick embryonic spinal cord cDNA, and partial fragments of mouse Robo3 and mouse Sim1 were amplified from mouse embryonic spinal cord cDNA. The primer sequences used are presented in Table S2.
The fragments were inserted into pCR-XL-TOPO (Thermo Fisher) or pBluescriptKS. For the reporter assay of putative cis-regulatory module (CRMs), genomic DNA fragments were amplified by PCR from a solution that was prepared by digesting the tail tip of a C57BL/6J mouse.
The amplified fragments were inserted into the pSF-pA-MinProm-eGFP vector (Oxford Genetics). The positions of the CRMs within the mouse reference genome are provided in Table S3. For the construction of pCAGGS-mCherry, the mCherry gene cassette of pmCherry-C1 (Clontech) was inserted into pCAGGS (Niwa et al., 1991).

| In ovo electroporation
A small window was opened on top of the fertilized chicken eggshell. The electrodes were placed on both sides of the neural tube of an HH12-13 chick embryo. Plasmid DNA was injected into the neural tube. During injection, electric pulses (25 V, 50 ms, five times, 950 ms interval) were applied using CUY21EDIT (BEX). The concentrations of the electroporated plasmids are listed in Table S5.

| Immunohistochemistry
Embryos were fixed in 0.1 M phosphate buffer/4% paraformaldehyde (PFA) at room temperature for 45-60 min or 4℃ overnight. The fixed embryos were then cryoprotected in 20% sucrose at 4℃ overnight, embedded in an OCT compound, and cryosectioned. In cases of overnight fixed samples, sections were boiled for 20 min in sodium citrate buffer (10 mM sodium citrate, 0.05% Tween 20, pH 6.0) and cooled to room temperature. After washing with phosphatebuffered saline (PBS) containing 0.1% Triton X-100 (PBST), the sections were incubated with primary antibodies at 4℃ overnight, washed with PBST three times for 5 min, and incubated with secondary antibodies for 1-2 h at room temperature. After being washed with PBST three times for 5 min, the slides were coverslipped in VECTASHIELD (Vector Laboratories). The antibodies used are presented in Table S4. In the cases of anti-Pax6 and anti-Lhx1, Can Get Signal (TOYOBO, NKB-401) was used as the diluent instead of PBST. Images were captured using an Olympus BX51 fluorescent microscope equipped with an Olympus DP71 digital camera or an Olympus FluoView FV1000 confocal microscope.

| In situ hybridization
For double staining of immunohistochemistry and in situ hybridization, immunohistochemistry was performed as described above with primary and secondary antibody incubation for 1 h at room temperature, and the signal was developed using the VECTASTAIN Elite ABC kit (Vector Laboratories). Subsequently, the process of in situ hybridization commenced, as described below. The plasmid for the RNA probe template was linearized, and digoxigenin (DIG)-labeled RNA probes were generated using a DIG RNA labeling mix (Sigma-Aldrich, 11277073910) and T3 RNA polymerase (Promega, P2083).
Tissue slides were washed with PBS for 5 min, treated with proteinase K (2 µg/ml) for 10 min, washed with PBS for 5 min, fixed in 4% PFA for 10 min, and washed three times with PBS for 5 min.
The slides were then incubated in an acetylation solution (100 mM triethanolamine, pH 8.0) twice for 2-3 min, transferred to a new acetylation solution, acetylated for 15 min by adding dropwise acetic anhydride (0.3% final concentration), and washed with PBS three times for 5 min. The slides were prehybridized with a hybridization buffer (50% formamide, 5× SSC, 5× Denhardt's solution, 500 µg/ml herring sperm DNA, and 250 µg/ml yeast RNA) for 2-3 h at room temperature, and hybridized with a hybridization buffer containing DIG-labeled RNA probe. The slides were coverslipped and incubated in a humidified chamber at 70℃ overnight. The following day, the slides were transferred to 5× SSC at 70℃ for 5 min, washed twice in 0.2× SSC for 30 min at 70℃, and washed in 0.2× SSC for 30 min at room temperature. The slides were washed with buffer B1 (100 mM Tris-HCl pH 7.5, 150 mM NaCl, 0.1% Tween 20) three times for 30 min, incubated with buffer B1 containing 10% heatinactivated normal goat serum for 2-3 h at room temperature, and then incubated with buffer B1 containing alkaline phosphatase conjugated anti-DIG antibody (Roche 11093274910, 1:5000) and 1% heat-inactivated normal goat serum at 4℃ overnight. The following day, the slides were washed five times with buffer B1 for 30 min, and incubated in buffer B3 (100 mM Tris-HCl pH 9.5, 100 mM NaCl, 50 mM MgCl2, 0.1% Tween20) three times for 5 min. The slides were incubated with BM purple (Roche, 11442074001) until the signal was visualized (3 h to overnight). Next, the slides were washed three times with buffer B1 for 30 min, washed three times with PBS for 5 min, and coverslipped in Fluoromount (Diagnostic Biosystems, K024).

| Two distinct subtypes of V3 INs in amniotes
To clarify whether post-mitotic neuronal differentiation has diverged among vertebrates, we compared publicly available scRNAseq data from mice and zebrafish. We focused on V3 INs because we detected diversity in the gene expression among different species, as discussed below. First, we examined scRNA-seq data obtained from the spinal cords of mouse embryos (Delile et al., 2019). The data consisted of the whole spinal cord cells from embryonic day (E) 9.5 to E13.5. We extracted data with V3 IN profiles (Nkx2-2 or Sim1 positive), applied graph-based clustering and dimensionality reduction by tSNE to these cells, and then identified the differentially expressed genes in each cluster (Table S1). The expression levels of several markers were visualized on the tSNE plots ( Figure 1a).
V3 INs are glutamatergic neurons, and thus express Slc17a6 (vGlut2) (Zhang et al., 2008). The expression of Sox2, Neurog3, and Tubb3 (class III β-tubulin) in this order indicates the transition from progenitor to differentiated neurons (Carcagno et al., 2014). Thus, the medial-to-lateral direction in the spinal cord corresponds to the top-to-bottom direction in this plot. We found at least two distinct populations in V3 INs: one positive for Robo3, Olig3, and Cntn2 (Tag1) (cluster 5 in Figure 1a), and the other positive for Lhx1 (cluster 1 in Figure 1a, Table S1).
To distinguish these two subtypes in vivo, markers of each subtype (Robo3 and Lhx1) were examined in the spinal cord of the mouse embryo at E11.5 (Figure 1b- laterally to the commissural axons; however, these cells did not express Robo3 (Figure 1f,g, and k). Rather, Robo3 was selectively expressed by V3 INs located medially to the commissural axons ( Figure 1f,g). In contrast, Lhx1 was expressed in the laterally located population of the V3 INs (Figure 1h-j). We performed the same examination using chick embryos and obtained identical results ( Figure S1). These observations confirmed that there are at least two subtypes of V3 INs, which are transcriptionally, spatially, and probably also functionally distinct populations in the developing spinal cord in amniotes.

| The gene expression profile of V3 IN in zebrafish was distinct from that in amniotes
Next, we analyzed the scRNA-seq data derived from zebrafish  Figures 1 and 2). These results indicate that the V3 INs of zebrafish are not equivalent to those of amniotes as far as gene expression is concerned.

| Distinct gene expression profiles in other ventral INs between mice and zebrafish
To determine whether the divergence of the gene expression pro-  Barrio et al., 2007;Hayashi et al., 2018;Kimura et al., 2008;Li et al., 2005;Panayi et al., 2010;Peng et al., 2007), we annotated five V2 IN subtypes on the mouse tSNE plot as follows: (1)  CSF-cNs were identified as a distinct cluster in zebrafish tSNE plots ( Figure 3b, gata3 + , pkd2l1 + , sox1a/b + ), while in the mouse tSNE plots, there were few CSF-cNs, which were scattered and did not make any recognizable clusters, suggesting that CSF-cNs were not yet fully differentiated in mice at E13.5 (Figure 3a, plot for Pkd2l1).
Indeed, amniotes' CSF-cNs are born at a significantly late phase of spinal cord development (E14 or later in mice), namely when gliogenesis commences (Petracca et al., 2016). In contrast, the CSF-cNs of zebrafish are born early, that is, simultaneously with other INs and MNs (Park et al., 2004;Shin et al., 2007). This heterochronic development of CSF-cNs is another example of the divergence of neuronal development in the vertebrate spinal cord (Petracca et al., 2016).
We also examined the gene expression profiles of V0 (Evx1/2 + ) and V1 (En1 + ) INs briefly, which were suggestive of divergent gene expression profiles in V0 and V1 INs between mice and zebrafish (details are provided in Figure S5). Taken together, these results suggest that the divergence of gene expression profiles is not a specific feature of V3 INs; rather, it is a generally observed feature of the ventral INs of the spinal cord.

| CRM in the Robo3 locus is not conserved in teleosts
As is shown in Figure 2  Notably, when sequence conservation was checked using Multiz alignments in the UCSC genome browser (Blanchette et al., 2004) and VISTA (Frazer et al., 2004) among vertebrates, Robo3-CRM was conserved in tetrapods, but not in teleosts (highlighted in yellow, comparing tracks labeled "Tetrapods" and "Teleosts" in Figure 4a).

| Non-conservation of Gli binding sites in Gli1 and Ptch1 loci in vertebrates
Next, we examined whether the upstream process before the progenitor fate specification also diverged among vertebrates.
Previous studies have reported that the impact of loss of Gli function on neuronal progenitors is different between mice and zebrafish (Chandrasekhar et al., 1999;Karlstrom et al., 2003;Tyurina et al., 2005;Vanderlaan et al., 2005;Wang et al., 2013). This raises the possibility that some Gli binding sites (GBSs), that is, the sites of action of Hh signaling, vary among vertebrates. To test this, Gli1 and Ptch1 loci were interrogated, as these genes are direct targets of Hh signaling. We searched for GBS in the mouse genome and checked sequence conservation using Multiz alignments in the UCSC genome browser and VISTA. Using Gli1 and Gli3 ChIP-seq data, many ChIP-seq peaks were observed in mouse Gli1 and Ptch1 loci near TSS and within introns (Figure 5a,b), confirming previous reports (Dai et al., 1999;Nishi et al., 2015). We scanned the genomic regions around the ChIP-seq peaks and found five and six GBSs matching the Gli binding motif in Gli1 and Ptch1 loci, respectively (Figure 5c,d).
ATAC-seq data confirmed that most of these GBSs were in an open chromatin state in neural cells (Figure 5a,b, tracks labeled "ATAC").
We found that most of these GBSs were not conserved among vertebrates. In the case of Gli1, one GBS was found only in mice, and four other GBSs were conserved only in mammals (Figure 5a,d). In the Ptch1 locus, two GBSs were conserved in vertebrates, but four others were not well conserved (Figure 5b,d and S7). These findings imply that the GBS is a highly flexible element.

| Non-conservation of the CRMs of progenitor fate specifying TFs in vertebrates
As the regulatory elements in Robo3, Gli1, and Ptch1 were not conserved among vertebrates, we expected that non-conservation of the regulatory elements could be found in other genes. Thus, we applied the same analysis to TFs functioning as progenitor fate determinants expressed in the neural tube. Indeed, we found that We identified several CRMs that were conserved only in tetrapods; however, except for Olig2-CRM1, the in vivo functions of these CRMs have not been examined so far. Thus, we carried out a reporter assay by means of the chick in ovo electroporation system.

| Divergent properties of post-mitotic neurons in vertebrates
First, in the present study, we compared the single-cell transcriptomes of mice and zebrafish to elucidate the divergence of postmitotic neuronal differentiation. We found that, in amniotes, V3 INs can be divided into two distinct subtypes. One is located medially and expresses Robo3, Olig3, and Cntn2, and the other is located laterally and expresses Lhx1 (Figure 1). We tried to make these V3 INs correspond with those of zebrafish, but were unsuccessful (Figures 1   and 2). Likewise, we also examined V2 INs and again found that the gene expression profiles of V2a INs differed between mice and zebrafish ( Figure 3). This indicates that the spinal cord ventral INs of amniotes are not equivalent to those of zebrafish, at least regarding gene expression profiles, which supports the hypothesis that the properties of post-mitotic neurons in the spinal cord have diverged among vertebrates. This is in contrast to the case of progenitor cells, in which each progenitor domain of zebrafish readily corresponds to that of amniotes ( Figure S3).
The role of V3 INs in mice is to secure a stable locomotor rhythm (Zhang et al., 2008). Several studies have also suggested that V3 INs contribute to left-right synchronous motor output, such as gallop F I G U R E 7 Enhancer functions of CRMs of progenitor fate-specifying TFs. The CRM reporter vectors indicated were electroporated into chick neural tubes together with CAGGS::mCherry as a control vector. (a-j) The GFP and mCherry expression was examined at HH19-20 in the forelimb-level neural tube. (k-o) The endogenous expression of the indicated genes was examined by in situ hybridization (l, m, n) or immunohistochemistry (k, o) together with the GFP expression. These five CRMs displayed enhancer functions, although the GFP expression domains incompletely overlapped with the endogenous expression domain. Pax6-CRM induced GFP expression almost ubiquitously, but its expression was not observed in the roof plate (a, f, k). Gsx1-CRM induced GFP expression in the dorsal neural tube, but not in the more ventral region than the endogenous expression domain (b, g, l). The GFP expression domains induced by Dbx2-CRM or Irx3-CRM2 overlapped with the endogenous expression domain almost completely (c, d, h, i, m, and n). Olig2-CRM2 induced GFP expression in the intermediate to ventral region, which partially overlapped with, but was more dorsal than, the endogenous expression domain (e, j, and o). The number of embryos examined in the electroporation experiments is presented in Table S5. Scale bar: 100 μm and bound (Danner et al., 2019;Kiehn, 2016;Rabe et al., 2009).
These gaits can be expressed by quadrupeds, but not by fish. This suggests that the different locomotor behaviors between amniotes and teleosts are associated with the divergence of V3 IN properties.
The same may be true for V2a INs, whose roles have diverged between amniotes and teleosts (Azim et al., 2014;Crone et al., 2008Crone et al., , 2009Eklöf-Ljunggren et al., 2012;Kiehn, 2016). Given that the number of muscles and corresponding MNs has significantly increased during tetrapod evolution, especially in limb muscles (Diogo & Abdala, 2010), additional layers of neuronal circuits are required to precisely control complex locomotor behavior (Kiehn, 2016). This may be accomplished, at least in part, by tinkering with the already existing neuronal population (Jacob, 1977), eventually leading to the divergence of post-mitotic neuronal properties among vertebrates.
Distinct expression profiles of Robo3 between amniotes and zebrafish may be a strategy for neuronal tinkering (Figures 1, 2, and 4).

| Limitations in cross-species comparisons of single-cell transcriptomes
Recently, several studies have reported cross-species comparisons of single-cell transcriptomes (Tosches et al., 2018;Zhu et al., 2018).
To achieve a comprehensive and quantitative comparison of singlecell transcriptomes, these studies integrated two datasets from distinct species after one-to-one ortholog identification and filtering out nonhomologous genes (Tosches et al., 2018;Zhu et al., 2018). In contrast, instead of integrating two datasets, in the present study, we analyzed the mouse and zebrafish scRNA-seq datasets separately for the following two reasons.
(1) Before the integration of the two transcriptome datasets, one-to-one ortholog identification is necessary. However, in the case of mice and zebrafish, this is difficult because of the whole genome duplication in teleost species.
(2) In the current study, we focused on non-conservation; thus, it would have been disadvantageous for us to filter out nonhomologous genes. Fortunately, an improved method of single-cell transcriptome integration has recently been reported, which does not rely on one-to-one ortholog identification, but rather utilizes the weighted gene-gene homology graph, and can even detect paralog substitutions (Tarashansky et al., 2021). Such a new method may help overcome these difficulties in future studies.
In the present study, we used publicly available mouse and zebrafish scRNA-seq data; however, these data were not fully comparable. The mouse dataset contained data for E9.5, 10.5, 11.5, 12.5, and 13.5 (Delile et al., 2019), while the zebrafish dataset contained data for 1, 2, and 5 dpf (Farnsworth et al., 2020; 5 dpf data were not included in our analysis, as that developmental stage is too advanced).
Given the rapid embryogenesis of zebrafish, data that included more time points and shorter intervals would be favorable. In the case of zebrafish, data with spinal cord identity were extracted from whole embryonic data, based solely on the gene expression profiles, not on the anatomical data. Thus, theoretically, we cannot rule out the possibility that cells outside of the spinal cord are misannotated as the spinal cord, and vice versa. To further corroborate our hypothesis in future studies, it would be advisable to acquire new zebrafish data, obtained in shorter intervals during the neurogenic phase (for example, 6 h intervals during 12-48 dpf) in combination with the isolation of the spinal cord, by taking advantage of reporter transgenic lines and/or by manual dissection.
For example, Zhu et al. (2018) reported that, even in closely related species (human and macaque), the neuronal transcriptomes in identical brain regions diverged between the two species.

| Divergent processes leading to the conserved progenitor domains
In the latter part of this study, we investigated the transcriptional regulatory elements located around the progenitor fate specifying genes in order to reveal the extent to which the upstream process before the progenitor specification has diverged in vertebrates. In the Gli1 locus, clustered GBSs near the TSS were conserved only in mammals, but not in teleosts ( Figure 5). Even in birds and reptiles, these GBSs are conserved only partially. Nevertheless, in vertebrates, Gli1 is commonly expressed in response to Hh signaling (Aglyamova & Agarwala, 2007;Karlstrom et al., 2003;Luo et al., 2006), suggesting the presence of distinct GBSs with equivalent functions. Indeed, zebrafish possess three GBSs in the Gli1 locus, which are conserved only in some teleost species (Wang et al., 2013), indicating that regions where Hh signaling is eventually transduced has shifted during vertebrate evolution. This is supported by a similar situation in the Ptch1 locus ( Figure 5 and S7; Wang et al., 2013). These findings may partly explain the discrepancy in Gli loss-of-function phenotypes between mice and zebrafish.
It has been pointed out that many features of Shh are similarly observed in bicoid, a morphogen that plays a crucial role in patterning the anterior-posterior axis in Drosophila blastoderm (Briscoe & Small, 2015). It is noteworthy that bicoid binding sites have also undergone a rapid turnover in Diptera (McGregor et al., 2001). The flexibility of morphogen response elements might contribute to the integration of morphogen dependency into the patterning system of the embryo (Dearden & Akam, 1999;Miyamoto & Wada, 2013;Ren et al., 2020;Stauber et al., 1999).
We identified functional CRMs bound by multiple TFs in the Pax6, Gsx1, Dbx2, Irx3, and Olig2 loci in the mouse genome ( Figures 6   and 7, Fig. S8-S11). These CRMs are conserved only in tetrapods, but are lost in teleosts (Figure 6, Fig. S8-S11). These findings indicate that, although the progenitor domain organization is conserved among vertebrates, the cis-regulatory elements contributing to it are not constrained (Figure 8). This is considered a case of developmental system drift (DSD; True & Haag, 2001). Several studies have reported similar situations, in which divergent regulatory sequences in different species result in conserved gene expression (Barrière et al., 2012;Domené et al., 2013;Fisher et al., 2006;Hare et al., 2008;Ludwig et al., 1998Ludwig et al., , 2000Paris et al., 2013;Stolfi et al., 2014;Swanson et al., 2011). In these cases, including the present study, the GRN architecture is likely to be maintained as a whole despite the cis-element turnover, as demonstrated in mammalian evolution (Stergachis et al., 2014;Vierstra et al., 2014).
The process of neural tube formation differs between amniotes and zebrafish not only genetically but also morphogenetically. In amniotes, a midline groove is formed by the bending of the neural plate, the edges of which are then fused, thus forming a neural tube (Schoenwolf & Smith, 1990). On the other hand, in zebrafish, a solid neural keel with no central canal is formed by the convergent movement of the neural plate cells, and then the central canal opens secondarily (Schmitz et al., 1993). It follows that neural tube formation is a case of DSD from both genetic and morphogenetic viewpoints.

| Molecular basis of the hourglass model in the spinal cord
Taken together, our results demonstrate that the divergence of the developmental process of the spinal cord is in accordance with the developmental hourglass model (Figure 8). What then imposes the hourglass-like pattern on the evolution of spinal cord development? The progenitor regionalization in the neural tube is a highly robust developmental system and is thus insensitive to stochastic noises of graded morphogen activities and genetic mutations (Balaskas et al., 2012;Delás & Briscoe, 2020;Exelby et al., 2021;Xiong et al., 2013;Zagorski et al., 2017). Simulation studies have also proposed that developmental system robustness is an emergent property of the complex GRN (Bergman & Siegal, 2003;von Dassow et al., 2000;Siegal & Bergman, 2002). Although the authors of the aforementioned studies used in silico simulations and did not deal F I G U R E 8 The hourglass-like pattern of the developmental divergence of the spinal cord. A summary of this study is presented. In this drawing, development proceeds from the bottom to the top. Only five progenitor domains were set up for simplification. Conserved and nonconserved features are colored in blue and orange, respectively. The bottom drawing represents GRNs regulating the progenitor domain establishment. CRMs (indicated by boxes) located in the progenitor fate specifying genes (GeneA-D) diverged (turnover) between amniotes and teleosts (orange boxes). Accordingly, these GRNs have been rewired. Nevertheless, these distinct GRNs result in the same progenitor domain organization. The top drawing represents the divergence of the differentiation process of post-mitotic neurons. After individual cells leave the progenitor domains as post-mitotic neurons (indicated by circles), some neurons undergo distinct maturation processes between amniotes and teleosts (orange circles). Thus, there exist neuronal subpopulations whose function is different between amniotes and teleosts with neural tube development, their results indirectly support the robustness of the neuronal progenitor specification system, which is organized by highly complex GRN (Kutejova et al., 2016). From an evolutionary perspective, such a robust (canalized) developmental system can tolerate genetic mutations without phenotypic effects (Rutherford & Lindquist, 1998;Waddington, 1942Waddington, , 1957. In other words, the progenitor specification system in the neural tube buffers genetic variations. We speculate that this situation consequently led to the turnover of cis-regulatory elements, while the progenitor arrangement was maintained. Once cells exit and migrate away from the progenitor domains as post-mitotic neurons, individual neurons depend on the mechanisms regulating post-mitotic maturation, which is distinct from the progenitor specification GRN. As mentioned previously, the maturation process is considered to diverge in association with locomotor divergence. Eventually, through vertebrate evolution, the developmental divergence of the spinal cord might lead to an hourglass shape (Figure 8). The current study focused on the spinal cord. Thus, it is of particular interest to examine whether this scenario is also the case in other developmental systems or even at the whole embryonic level.
This suggests that the developmental system robustness of neuroectodermal regionalization has already been acquired in the last common ancestor of bilaterians, and that upstream signals of the neuroectodermal regionalization had been modified after the divergence of phyla (e.g., the ventralization factor is dl in Drosophila or Shh in vertebrates; Cornell & von Ohlen, 2000;von Ohlen & Doe, 2000). This is supported by the notion of the common origin of the central nervous system (CNS) in bilaterians (Arendt, 2018;Arendt et al., 2016;Denes et al., 2007), but is inconsistent with the convergent evolution of the CNS (Martín-Durán et al., 2018). Supporting the common evolutionary origin of the bilaterian CNS, we suggest the following evolutionary scenario. The trunk nerve cord, which develops from the regionalized neuroectoderm, has already been acquired in the common ancestor of bilaterians, and its development was so canalized that upstream regulators could be modified, while conserving the regionalized neuroectoderm (e.g., Hh signaling recruitment in deuterostomes; Miyamoto & Wada, 2013;Ren et al., 2020). Therefore, distinct upstream regulators of neuroectoderm regionalization among bilaterians can be considered to be DSD that had taken place during more than 500 million years of bilaterian evolution.

ACK N OWLED G M ENTS
We thank Y. Watanabe and T.