Molecular Analyses of V0v Spinal Interneurons and Identification of Transcriptional Regulators Downstream of Evx1 and Evx2 in these Cells

Background V0v spinal interneurons are highly conserved, glutamatergic, commissural neurons that function in locomotor circuits. We have previously shown that Evx1 and Evx2 are required to specify the neurotransmitter phenotype of these cells. However, we still know very little about the gene regulatory networks that act downstream of these transcription factors in V0v cells. Methods To identify candidate members of V0v gene regulatory networks, we FAC-sorted WT and evx1;evx2 double mutant zebrafish V0v spinal interneurons and expression-profiled them using microarrays and single cell RNA-seq. We also used in situ hybridization to compare expression of a subset of candidate genes in evx1;evx2 double mutants and wild-type siblings. Results Our data reveal two molecularly distinct subtypes of V0v spinal interneurons at 48 h and suggest that, by this stage of development, evx1;evx2 double mutant cells transfate into either inhibitory spinal interneurons, or motoneurons. Our results also identify 25 transcriptional regulator genes that require Evx1/2 for their expression in V0v interneurons, plus a further 11 transcriptional regulator genes that are repressed in V0v interneurons by Evx1/2. Two of the latter genes are hmx2 and hmx3a. Intriguingly, we show that Hmx2/3a, repress dI2 interneuronal expression of skor1a and nefma, two genes that require Evx1/2 for their expression in V0v interneurons. This suggests that Evx1/2 might regulate skor1a and nefma expression in V0v interneurons by repressing Hmx2/3a expression. Conclusions This study identifies two molecularly distinct subsets of V0v spinal interneurons, as well as multiple transcriptional regulators that are strong candidates for acting downstream of Evx1/2 to specify the essential functional characteristics of these cells. Our data further suggest that in the absence of both Evx1 and Evx2, V0v spinal interneurons initially change their neurotransmitter phenotypes from excitatory to inhibitory and then, later, start to express markers of distinct types of inhibitory spinal interneurons, or motoneurons. Taken together, our findings significantly increase our knowledge of V0v and spinal development and move us closer towards the essential goal of identifying the complete gene regulatory networks that specify this crucial cell type.


Abstract Background
V0v spinal interneurons are highly conserved, glutamatergic, commissural neurons that function in locomotor circuits.
We have previously shown that Evx1 and Evx2 are required to specify the neurotransmitter phenotype of these cells.
However, we still know very little about the gene regulatory networks that act downstream of these transcription factors in V0v cells.

Methods
To identify candidate members of V0v gene regulatory networks, we FAC-sorted WT and evx1;evx2 double mutant zebra sh V0v spinal interneurons and expression-pro led them using microarrays and single cell RNA-seq. We also used in situ hybridization to compare expression of a subset of candidate genes in evx1;evx2 double mutants and wild-type siblings.

Results
Our data reveal two molecularly distinct subtypes of V0v spinal interneurons at 48 h and suggest that, by this stage of development, evx1;evx2 double mutant cells transfate into either inhibitory spinal interneurons, or motoneurons. Our results also identify 25 transcriptional regulator genes that require Evx1/2 for their expression in V0v interneurons, plus a further 11 transcriptional regulator genes that are repressed in V0v interneurons by Evx1/2. Two of the latter genes are hmx2 and hmx3a. Intriguingly, we show that Hmx2/3a, repress dI2 interneuronal expression of skor1a and nefma, two genes that require Evx1/2 for their expression in V0v interneurons. This suggests that Evx1/2 might regulate skor1a and nefma expression in V0v interneurons by repressing Hmx2/3a expression.

Conclusions
This study identi es two molecularly distinct subsets of V0v spinal interneurons, as well as multiple transcriptional regulators that are strong candidates for acting downstream of Evx1/2 to specify the essential functional characteristics of these cells. Our data further suggest that in the absence of both Evx1 and Evx2, V0v spinal interneurons initially change their neurotransmitter phenotypes from excitatory to inhibitory and then, later, start to express markers of distinct types of inhibitory spinal interneurons, or motoneurons. Taken together, our ndings signi cantly increase our knowledge of V0v and spinal development and move us closer towards the essential goal of identifying the complete gene regulatory networks that specify this crucial cell type.

Background
For the Central Nervous System (CNS) to operate correctly, neurons with appropriate functions need to be precisely created and accurately connected into circuits. However, we still do not understand how this vital aspect of neural development is achieved. The spinal cord is a powerful system for establishing fundamental principles of neuronal fate speci cation and circuit assembly, as it is relatively simple and experimentally tractable compared to the brain. It is also an essential part of the CNS as the spinal cord controls locomotion and receives and processes sensory information from the trunk and limbs. In addition, spinal cord dysfunction caused by abnormal development, injury or disease can profoundly impair quality of life. Therefore, it is essential that we better understand neuronal speci cation in the spinal cord, so that we can develop more effective therapies to treat these debilitating conditions.
To elucidate how spinal cord circuitry develops, we rst need to establish how neuronal functional properties are speci ed, as these properties determine the circuits that speci c neurons participate in and their functions within those circuits. Most of the neurons in the spinal cord are interneurons, so called because their cell bodies and axons reside entirely within the CNS. Interneurons have essential roles within most spinal circuits. One of the most important functional properties that helps to de ne distinct interneurons, and their speci c functions in neural circuitry, is which neurotransmitter they use to communicate with other cells. Spinal interneurons use three major neurotransmitters: glutamate, which is excitatory; glycine, which is inhibitory; and GABA, which, with a few exceptions at particular stages of development, is usually inhibitory. If these neurotransmitter phenotypes are wrongly speci ed then the affected interneurons will inappropriately inhibit rather than excite, or vice versa, their synaptic partners within neuronal circuits, and the functional outputs and behaviors of those circuits will be dramatically disturbed, usually with pathological consequences (1).
All of the data so far, suggest that neurotransmitter phenotypes, and other aspects of cell fate speci cation are determined by the transcription factors that an interneuron expresses when it becomes post-mitotic and starts to differentiate (e.g. [2][3][4][5][6][7]. In most cases analyzed so far, several transcription factors act together, or in succession, as part of a Gene Regulatory Network (GRN) that speci es a particular functional property (e.g. [3][4][5][6]. We already know several transcription factors that are part of GRNs that specify spinal interneuron inhibitory (4)(5)(6)(7)(8)(9)(10) or excitatory (4,8,(11)(12)(13)(14) phenotypes. However, there are still fundamental gaps in our knowledge of how spinal interneuron neurotransmitter phenotypes are speci ed and maintained. For example, it is very unlikely that we have identi ed all the transcription factors that play crucial roles in the GRNs in any speci c class of spinal interneurons, or all of their epistatic relationships.
In this study, we concentrate on V0v spinal interneurons. These are highly conserved glutamatergic, commissural neurons that exist in all vertebrates examined so far. V0v spinal interneurons are located in the middle of the dorsal/ventral axis of the spinal cord, and they are required for correct left-right alternation during fast locomotion (15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26). Given their important role in locomotion, it is vital that we understand how V0v interneurons are speci ed, and, in particular, how they adopt their excitatory neurotransmitter fate. Within both the mouse and zebra sh spinal cord, two highly related transcription factors, Evx1 and Evx2 (Evx1/2), are exclusively expressed in V0v interneurons, and, crucially, these Evx transcription factors are required for the glutamatergic neurotransmitter phenotypes of these cells (14,17). However, our previous analyses of zebra sh evx1;evx2 double mutants did not detect changes in other V0v interneuron functional characteristics, such as axon trajectory, and, unlike in previous mouse studies, V0v interneurons in evx1;evx2 double mutant zebra sh did not trans-fate into V1 interneurons (14,17). These data suggest that Evx1/2 are speci cally required for the glutamatergic phenotype of V0v spinal interneurons in both zebra sh and mouse, but that these transcription factors may have additional functions in mammalian V0v spinal interneurons.
Alternatively, as we only analyzed early developmental stages in our zebra sh studies, it is also possible that the zebra sh evx1;evx2 mutant phenotype is more similar to that of mouse at later stages. To explore this possibility, in this study, we provide the rst temporal analysis of V0v mutant phenotypes in any species.
In previous studies, we have started to de ne the GRNs downstream of Evx1/2 in V0v neurotransmitter fate speci cation. We have shown that Lmx1ba and Lmx1bb (Lmx1ba/b) act downstream of Evx1/2 in the speci cation and/or maintenance of V0v glutamatergic fates, and Evx1/2, but not Lmx1ba/b, are also required to repress inhibitory phenotypes in V0v interneurons (13,14). These data suggest that there may be two different important GRNs downstream of evx1/2 in V0v spinal interneurons, one that speci es glutamatergic fates and one that represses inhibitory / glycinergic fates. Tg(hmx CNEIII:cfos:GAL4-VP16,UAS:EGFP) SU41 line which we used in the experiments in this paper. Note though that this line also does not recapitulate endogenous hmx3a expression in either the ear or lateral line primordium, suggesting that the enhancer region(s) driving expression in these tissues is not present in CNE III (data not shown).
in situ Hybridisation and Immunohistochemistry We xed embryos in 4% paraformaldehyde/phosphate-buffered saline (PBS) and performed single in situ hybridizations and immunohistochemistry plus in situ hybridization double-labelling experiments as previously described (37,38). Sources of in situ hybridization probes are provided in Supp. Table 1. PCR-based in situ probes were created with cDNA from 27 h WT zebra sh embryos. We extracted total RNA by homogenizing 50-100 mg of embryos in 1 ml of TRIzol reagent (Ambion, 15596026). We con rmed RNA integrity (2:1 ratio of 28S:18S rRNA bands) and quality (A260/A280 ratio of ~ 2.0) using agarose gel electrophoresis and spectrophotometry respectively.
We synthesized cDNA using Bio-Rad iScript Reverse Transcription Supermix kit (Bio-Rad, 1708891). We ampli ed each sequence using Phusion High-Fidelity DNA Polymerase (M0530L, NEB) and 5 µl cDNA in 50 µl total reaction volumes, using the primers and annealing temperatures shown in Supp. Table 1. To avoid cross-reactivity, whenever possible, riboprobes were designed against 3'UTR or coding sequence lacking all conserved protein domains in Pfam (39). Primers were designed using Primer3 web version 4.1.0 at https://primer3.ut.ee (40,41) and the following design parameters: optimum primer size: 22 bp (minimum: 20 bp, maximum: 25 bp), optimum annealing temperature: 58.0 o C (minimum: 57.0 o C, maximum: 60.0 o C), and optimum GC content: 50% (minimum: 40%, maximum: 60%). The preferred product size range was 800-1100 bp. This was not always possible, if there was little or no novel coding and/or 3' UTR sequence available (see Supp. Table 1). The PCR conditions were: 98.0 o C for 30 seconds, 35 cycles of: 98.0 o C for 10 seconds; annealing temperature in Supp. Table 1 for 30 seconds and 72.0 o C for 30 seconds, followed by a nal extension for 5 minutes at 72.0 o C. The PCR product was assessed on a 1% TAE gel, before purifying through phenol:chloroform:isoamyl alcohol extraction and precipitation with 0.2 M NaCl and ice-cold ethanol. If non-speci c banding was generated in addition to the desired PCR product, the speci c product was puri ed from the agarose gel using the Monarch DNA Gel Extraction Kit (NEB, T0120S). Each reverse primer contains the T3 RNA Polymerase minimal promoter sequence (shown in bold and underlined in Supp. Table 1). in situ probe synthesis was performed using 1 µg puri ed PCR product, T3 RNA Polymerase (Roche, 11031171001) and DIG RNA Labeling Mix (Roche, 11277073910). Embryos older than 24 h were usually incubated in 0.003% 1-phenyl-2-thiourea (PTU) to prevent pigment formation.
For some experiments we added 5% of Dextran Sulfate to the hybridization buffer (* in Table 1). Dextran sulfate can increase speci c staining in in situ hybridization experiments as it facilitates molecular crowding (42,43).
In cases where we did not detect expression of a particular gene in the spinal cord, we checked for low levels of expression by exposing embryos to prolonged staining. In some cases, this produced higher background (diffuse, nonspeci c staining), especially in the hindbrain, where ventricles can sometimes trap anti-sense riboprobes. Sometimes we were unable to perform double uorescent staining experiments due to very weak labelling with our RNA probes. In these cases, we combined single uorescent labelling with NBT/BCIP chromogenic staining. However, whereas stronger in situ signals from weak RNA probes can be obtained using NBT/BCIP, visualization of coexpressing cells becomes more di cult with this method. To preserve the integrity of both the NBT/BCIP chromogenic in situ signal and the weaker IHC signal, embryos were stored and mounted in VECTASHIELD Antifade Mounting Medium (Vector Laboratories, H-1000-10).

Fluorescent-Activated Cell Sorting (FACS)
The microarray expression pro ling experiments are described in detail in (13,44). P-values were corrected for multiple testing (45)(46)(47). These data have been deposited in the NCBI Gene Expression Omnibus with accession number GSE145916. For all these experiments, embryos were deyolked, and trunks were dissected and dissociated as described in (44), with the following modi cations: Trunk tissue was dissected anteriorly at the boundary between the hindbrain and spinal cord, and posteriorly, immediately above the end of the yolk extension. Embryos were processed in batches of 50 and stored on ice for a maximum of two hours prior to dissociation, to preserve cell and mRNA viability. D8537) xative. A tube of EGFP-negative cells was also collected to assess xation e ciency. Sorted cells were incubated on ice for 1 hour before assessing xation e ciency of the EGFP-negative control tube using Trypan Blue (ThermoFisher Scienti c, 15250061) and a hemocytometer. Samples with intact, fully xed cells, containing little or no cell debris were stored at + 4 o C for up to six days prior to rehydrating and performing single-cell capture with the 10x Genomics Chromium system (see below).

Single-Cell RNA-Seq
To rehydrate our xed EGFP-positive 48 h evx1 i232 ;evx2 sa140 ;Tg(evx1:EGFP) SU2 cells (stored at + 4 o C for up to six days post-FACS -see above), we rst centrifuged each sample at 300 rcf for 10 minutes at + 4 o C using a swing-bucket centrifuge. We recommend against using a xed rotor centrifuge as this can severely reduce recovery yields. Note that we never observed cell pellets during this step. Therefore, we marked the outer side of each 5 ml round-bottomed sample tube prior to centrifuging and avoided decanting from this side of the tube to prevent disruption and loss of the cell pellet. Next, we carefully removed most of the supernatant with a sterile p1000 tip, until approximately 100 µl remained in the tube. Samples were always kept on ice. Each cell pellet was then gently resuspended by adding 2 ml of freshly made, pre-chilled Rehydration Buffer (1x Dulbecco's Phosphate-Buffered Saline, no calcium, no magnesium (Merck, D8537), 1.0% UltraPure BSA (ThermoFisher Scienti c, AM2616), 0.5 u/µl Roche Protector RNase Inhibitor (Merck, 3335402001)) and gently pipetting 10 times. It is important to avoid making foam. We repeated the centrifugation and resuspension in Rehydration Buffer steps as previously. After the second Rehydration step, we again centrifuged at 300 rcf for 10 minutes at + 4 o C before carefully removing all but 30-40 µl of supernatant. Using a sterile p200 tip, we carefully resuspended the cell pellet and immediately measured the cell concentration in triplicate using a Bio-Rad TC20 automated cell counter (Bio-Rad, 1450102). We also checked a small aliquot using a compound microscope to ensure we had single cell suspensions. As described by 10x Genomics in their Sample Preparation Demonstrated Protocol "Methanol Fixation of Cells for Single Cell RNA Sequencing" (https://www.10xgenomics.com), we too recovered approximately 50% of the sorted cells after rehydration. Therefore, we rehydrated cells from four separate FACS experiments prior to performing single-cell capture.
We isolated single cells using a 10x Genomics Chromium system, aiming for capture of 10,000 cells per well (Chromium Next GEM Chip G Single Cell Kit, 1000127). We loaded 4 wells in total. This, and all subsequent library preparation steps were performed at the SUNY Upstate Medical University Molecular Analysis Core. We prepared libraries using a 10x Genomics Chromium Next GEM Single Cell 3' GEM, Library and Gel Bead Kit (v3.1, 10x Genomics, 1000128) and sequenced them on an Illumina NextSeq500 to a depth of at least 50,000 reads per cell (Illumina NextSeq 500/500 High Output Kit, v2. 5, 150 cycles, 20024907). We then performed demultiplexing and counts analysis as per the manufacturer's instructions using Cell Ranger v4.0.0 software (https://www.10xgenomics.com) and the Lawson Lab zebra sh transcriptome annotation model V4.3.2 (50). We analyzed the data using Partek Flow Genomic Analysis Software (51). Multiplets were removed by ltering out cells with > 12,000 counts and > 2,500 detected genes. Sick and/or "leaky" cells were removed by ltering out cells with < 500 detected genes and > 6% mitochondrial transcripts. We normalized the data using a counts per million (CPM) algorithm and applied a logarithmic transformation to improve data visualization. The outcome of normalization was assessed by principal components analysis (PCA), graph-based clustering and Uniform Manifold Approximation and Projection (UMAP) plotting, using the NN-Descent method of nearest neighbour type calculation and Euclidean distance metrics. We manually inspected UMAP plots to assess clustering quality based on expression of known V0v spinal interneuron markers. We excluded immature spinal cells that had begun to express the transgene but otherwise lacked expression of V0v post-mitotic genes to focus our analysis only on post-mitotic, differentiated V0v cells. We then ne-tuned the clustering by manually deducing and extrapolating cell fate assignments by comparing expression pro les of 48 h single-cell clusters with the molecular phenotypes of V0v spinal interneurons in 24 h and 30 h wild-type, evx1 i232/i232 and evx2 sa140/sa140 single mutant and evx1 i232/i232 ;evx2 sa140/sa140 double mutant embryos, as described by Juárez-Morales and colleagues (14) and data in this study (see Results). To perform gene-speci c analyses (GSA) of differential expression, we used the statistically robust Hurdle Model with default parameters in Partek Flow (51).
Under these conditions, the Hurdle Model in Partek Flow is equivalent to the widely used Model-based Analysis of Single-cell Transcriptomics (MAST) framework, which also incorporates Hurdle modelling (52). Hurdle models deal e ciently with the sources of nuisance variation commonly associated with single-cell datasets, such as sparsely detected cells (which in uence the cellular detection rate, an indicator of technical and/or biological variability between samples) and bimodal gene expression values (where many genes have zero expression values in the matrix, which can bias the interpretation of how much genes above the detection threshold are really expressed) (53,54).
When we examined our single-cell data, we observed that there were several small subsets of cells within mutant group 3 with distinct transcriptional pro les. Each of these subsets had so few cells, that we were concerned that this would underpower the Hurdle model by compromising the effectiveness of the variance modelling. To overcome this limitation, and further aid determination of differential expression, we also performed analysis of variance (ANOVA) using default parameters in Partek Flow (51). Unlike the Hurdle model, ANOVA models the expression of each gene independently of all the others, and Nault and colleagues have shown that it is the best method for calculating differential expression in scRNA-seq data when cell numbers are small (55). Therefore, we provide the data from both Hurdle modelling and ANOVA for analytical rigor.
We analyzed each of our four libraries separately. We omitted two of our libraries from further downstream analysis due to signi cant presence of notochord cells (which variably and ectopically express our Tg(evx1:EGFP) SU2 construct). This ectopic expression was much less abundant in our remaining two libraries and so these were used for the analysis shown in this paper. The sequencing depth for these two libraries approaches saturation (81.9% for one library and 89.3% for the other library), so that we have the best chance of detecting transcripts expressed at very low levels. We combined the data from these two libraries using the Counts Aggregation pipeline in Cell Ranger v4.0.0 and reanalyzed the data as described above. For the combined data, we identi ed 2860 cells that passed quality controls and V0v cell fate assignment (see Results).
Bulk RNA-Seq EGFP-positive cells were FAC-sorted from 27 h uninjected control and hmx2;hmx3a DKD Tg(hmx CNEIII:cfos:GAL4-VP16,UAS:EGFP) SU41 embryos as described above. RNA extractions were performed using a method based on that of (56) with the following modi cations. Prior to performing RNA extractions, all work surfaces and pipettors were treated with RNaseZAP (ThermoFisher Scienti c, AM9780). Throughout the process, samples were stored on ice unless otherwise stated. Frozen FAC-sorted cell lysates were removed from storage at -80 o C and thawed in a 37 o C water bath, before transferring to sterile microcentrifuge tubes. If necessary, sample volumes were increased to 250 µl with UltraPure DNase/RNase-Free distilled water (ThermoFisher Scienti c, 10977035). 750 µl TRIzol LS Reagent (ThermoFisher Scienti c, 10296028) was added to each 250 µl sample, before homogenising by gently pipetting up and down ten times with a sterile p1000 pipette tip. Samples were immediately transferred to Phasemaker tubes (which had been pre-centrifuged as per the manufacturer's instructions (ThermoFisher Scienti c, A33248)), before incubating for 5 minutes at room temperature. 200 µl chloroform was added to each sample. The tubes were then shaken vigorously for 15 seconds and incubated for a further 5 minutes at room temperature. The samples were then The data was analyzed using Partek Flow Genomic Analysis Software (51). We rst performed pre-alignment quality control assessment and recovered a minimum average read length of 74 bases and a minimum average read quality (Phred score) of 33.99, suggesting we recovered high quality, accurate sequencing data. Next, we trimmed the adapter sequence "CTGTCTCTTATACACATCT" from the 3' end using default parameters, before trimming bases from the 5' end. We selected an end minimum quality value (Phred) score of 32, and a minimum read length of 65 bases. Consequently, we trimmed an average of 1.14-1.15 bases. We aligned reads using default parameters and the STAR-2.6.1d algorithm, together with the Lawson Lab zebra sh transcriptome annotation model V4.3.2 (50). We aligned a minimum of 95.54% of all reads, with a minimum of 92.31% of reads aligning uniquely to the genome. Of these, a minimum 81.28% of reads aligned fully within an exon, a maximum 6.59% of reads aligned partly within an exon, a maximum 4.34% of reads aligned fully within an intron, and a maximum of 7.79% of reads were fully intergenic. We normalized the log expression ratios using a Trimmed Means of M-values (TMM) weighted algorithm (Robinson & Oshlack, 2010). We performed differential expression analysis using the Gene-Speci c Analysis (GSA) algorithm in Partek Flow. We used GSA because it makes no assumptions in advance about the data distribution nor the model choice necessary to deal with any nuisance factors present in the data. Rather, GSA describes transcript expression by calculating the data distribution and appropriate statistical model for each transcript in turn. As such, GSA can yield more accurate and reproducible expression data across the entire dataset, rather than just for the most pronounced expression outliers, as may be obtained with more common differential expression tools such as DEseq2 and limma (see Partek Flow Gene-Speci c Analysis white paper: https://documentation.partek.com/display/FLOWDOC/Genespeci c+Analysis). The outcome of GSA was assessed by hierarchical clustering (heatmap) plotting, clustering by features, using average linkage and Euclidean cluster distance and point distance metrics respectively.
qRT-PCR Analyses EGFP-positive and EGFP-negative cells were FAC-sorted from 27 h Tg(hmx CNEIII:cfos:GAL4-VP16,UAS:EGFP) SU41 embryos as described above. Total RNA was extracted as per the protocol used for our bulk RNA-seq experiments, with the exception that the nal eluted total RNA was divided in to 2-3 µl RNA aliquots in sterile PCR tubes and stored at -80 o C. cDNA was synthesized using the Bio-Rad iScript Reverse Transcription Supermix kit (Bio-Rad, 170-8891) and 3 µl of puri ed total RNA. We also included controls lacking Reverse-Transcriptase to assay for the presence of genomic DNA contamination. qRT-PCR was performed in triplicate for each sample using iTaq Universal SYBR Green Supermix (1725121, Bio-Rad) and a BioRad CFX96 real-time PCR machine. The following qPCR primers were used: The slc17a6b and slc32a1 primers were generated in this study. We generated the hmx2 and hmx3a primers in a previous study (12). The actb2 primers were generated by Hu and colleagues (57). To generate ampli cation data the program used was: 95. Imaging Embryos from single NBT/BCIP in situ hybridization experiments were mounted in 70% glycerol:30% distilled water between coverslip sandwiches (24 mm x 60 mm coverslips; VWR, 48393-106), with 2-4 coverslips (22 mm x 22 mm; VWR, 16004-094) on either side of the sample to avoid sample compression. Differential Interference Contrast (DIC) pictures were taken using an AxioCam MRc5 camera mounted on a Zeiss Axio Imager M1 compound microscope. Embryos from uorescent in situ hybridization + immunohistochemistry experiments were mounted in VECTASHIELD Antifade Mounting Medium (Vector Laboratories, H-1000-10) between coverslip sandwiches. Fluorescent images were taken on a Zeiss LSM 710 confocal microscope. Embryos from NBT/BCIP in situ hybridization + uorescent immunohistochemistry experiments were also mounted in VECTASHIELD Antifade Mounting Medium (Vector Laboratories, H-1000-10) between coverslip sandwiches. NBT/BCIP and uorescent images were captured using the T-PMT and 488 nm channels respectively on a Zeiss LSM 710 confocal microscope. Images were processed using Adobe Photoshop software (Adobe, Inc) and Image J software (58). NBT/BCIP confocal images (captured from in situ hybridization + immunohistochemistry experiments) are grayscale and were subsequently pseudo-colored in Photoshop by converting the image mode from Grayscale to Duotone. A custom purple ink tone (R = 48, G = 5, B = 107) was then applied and the image mode switched once more to RGB. The coloring now reproduces that of endogenous NBT/BCIP staining. In some cases, different focal planes were merged to show labelled cells at different mediallateral positions in the spinal cord. All images were processed for brightness-contrast and color balance using Adobe Photoshop software (Adobe, Inc.). Images of control and mutant embryos from the same experiment were processed identically. Figures were assembled using Adobe Photoshop (Adobe, Inc.).

Cell Counts and Statistics
In all cases, cell counts are for both sides of a ve-somite length of spinal cord adjacent to somites 6-10. Embryos were mounted laterally with the somite boundaries on each side of the embryo exactly aligned and the apex of the somite over the middle of the notochord. This ensures that the spinal cord is straight along its dorsal-ventral axis and that cells in the same dorsal-ventral position on opposite sides of the spinal cord will be directly above and below each other. Embryos from mutant crosses were counted blind to genotype. Labelled cells in embryos analyzed by Differential Interference Contrast (DIC) microscopy were counted while examining embryos on a Zeiss Axio Imager M1 compound microscope. We adjusted the focal plane as we examined the embryo to count cells at all medial-lateral positions (both sides of the spinal cord; also see (6, 13, 14, 37, 59)).
In some cases, cell count data were pooled from different experiments. Prior to pooling, all pairwise combinations of data sets were tested to determine if there were any statistically signi cant differences between them as described below. Data were only pooled if none of the pairwise comparisons were statistically signi cantly different from each other. In addition, as in situ hybridization staining can vary slightly between experiments, we only compared different mutant results when the counts from their corresponding WT sibling embryos were not statistically signi cantly different from each other.
To determine whether differences in values are statistically signi cant, data were rst analyzed for normality using the Shapiro-Wilk test. Data sets with non-normal distributions were subsequently analyzed using the Wilcoxon-Mann-Whitney test (also called the Mann Whitney U test). For data sets with normal distributions, the F-test for equal variances was performed, prior to conducting a type 2 (for equal variances) student's t-test. P-values generated by Wilcoxon-Mann-Whitney test and type 2 student's t-tests are indicated by ^ and + . Data are depicted as individual value plots and the n-values for each experimental group are also shown. For each plot, the wider red horizontal bar depicts the mean, and the red vertical bars depict the standard error of the mean (standard error of the mean (S.E.M.) values are listed in Table 1). Individual data value plots were generated using Prism version 9.4.0 (GraphPad Software, San Diego, California USA, www.graphpad.com). To assess whether mutant phenotypes occurred at Mendelian frequencies, we performed Chi-squared tests. Shapiro-Wilk and Wilcoxon-Mann-Whitney testing was performed in R version 3.5.1 (60). The F-test, student's t-test, and Chi-squared test were performed in Microsoft Excel version 16.62. Data and Reagent Availability Plasmids and zebra sh strains are available upon request. Microarray data were previously deposited in the NCBI Gene Expression Omnibus with accession number GSE145916. Single-cell and bulk RNA-Seq data have been deposited in the NCBI Gene Expression Omnibus with accession numbers GSE240239 and GSE240238 respectively. Results skor1a, skor1b, skor2, ebf3a, uncx, nefma, nefmb neff1, andinabare all expressed in the V0v spinal cord region.
We have previously analyzed lmx1ba and lmx1bb expression and function in V0v interneurons (13). To investigate the other genes, we performed an in situ hybridization time-course of skor1a, skor1b, skor2, ebf3a, uncx, uncx4.1, nefma, neff1, and inab expression in wild-type (WT) embryos to further con rm that they are expressed in the V0v region of the spinal cord. We analyzed 17, 20, 24, 36 and 48 h, as, within the spinal cord, evx1 and evx2 are expressed exclusively by V0v interneurons at all these time points, and all the data so far suggest that transcription factor genes important for specifying functional characteristics of spinal interneurons are expressed during these key stages of development (e.g. 5,12,13,14,37,70,71). Given that duplicated genes retained from whole genome duplication events (known as ohnologs) are often expressed in similar domains and may be functionally redundant, we also analyzed the spinal cord expression of inaa, ebf3b, and nefmb, ohnologs of inab, ebf3a, and nefma respectively. inaa did not show statistically signi cant differential expression in V0v interneurons in our microarray analysis (P > 0.17; Supp. Figure 1A'), and probes for nefmb and ebf3b were not present on the microarray because these genes were not accurately annotated in the Zv8 version of the zebra sh genome used for microarray construction. neff1 does not have an ohnolog in zebra sh (http://ohnologs.curie.fr/). We found that inaa and ebf3b are not expressed in the spinal cord at any of the stages examined (Supp. Figure 1A, C and data not shown). In contrast, all the other genes are expressed in the region where V0v cells are located (middle of the dorsal-ventral spinal cord axis, see evx1 expression Fig. 2A-E), during at least some of these crucial developmental stages. However, unlike evx1 and evx2, all of these genes are also sometimes expressed in other dorsal-ventral regions of the spinal cord, demonstrating that they are also expressed by at least one other spinal cord cell type (Fig. 2). Even skor1b, which is mainly expressed in the V0v domain, is also expressed by a few cells dorsal to V0v interneurons (Fig.  3A).
With respect to the three skor genes, only skor1b is expressed in the V0v spinal cord region at 17 h (cf. Figure 2A & 2K). skor1a (Fig. 2F) and skor2 (Fig. 2P) are both expressed in the spinal cord at this stage but in a more dorsal location. skor1b continues to be expressed in a similar dorsal-ventral spinal cord region to evx1 at 20 h, 24 h, 36 h and 48 h (cf. Figure 2L-O to B-E). By 20 h, skor1a is also expressed in the V0v spinal cord region, although it is still expressed in additional regions (cf. Figure 2G & 2B). In contrast, skor2 is still only expressed in the dorsal spinal cord (Fig. 2Q). By 24 h, all three skor genes are expressed in a similar spinal cord region to evx1 and this expression persists through 48 h (cf. Figure 2H-J, M-O, R-T to C-E). At 24 h, skor1a is still expressed in the dorsal spinal cord and it is also transiently expressed in a subset of ventral spinal cord cells (Fig. 2H). Expression in both these domains is lost by 36 h (Fig. 2I). skor2 is also still expressed in the dorsal spinal cord at 24 h ( Fig. 2R) but this dorsal expression begins to diminish by 36 h. (Fig. 2S-T). Taken together, these data suggest that these skor genes have distinct temporal patterns of expression in the V0v domain, with skor1b expression preceding skor1a and skor2 expression by 3 h and 7 h respectively. This raises the possibility that these genes have epistatic relationships with each other and/or function in different aspects of V0v cell development. ebf3a and uncx are expressed in a similar dorsal-ventral region of the spinal cord to evx1 at all the stages that we examined (cf. Figure 2U-Y & Z-AD to A-E). However, ebf3a is consistently expressed in more cells than evx1, suggesting that it is expressed by other cell types, in addition to V0v interneurons (cf. Figure 2U-Y to 2A-E). In contrast, uncx is only expressed by more cells than evx1, including cells in the ventral-most spinal cord, from 36 h onward (cf. Figure 2AC-AD to D-E). As previously described, uncx4.1 is strongly expressed in the somites at 17 h (Fig. 2AE, (67)). Nittoli and colleagues documented uncx4.1 expression in somites and brain at several different developmental stages but did not report spinal cord expression at any of the stages that they examined (67). Initially, we also did not detect spinal cord expression in our in situ hybridization experiments. However, when we used the molecular crowding agent Dextran Sulphate (see Methods), we were able to detect weak spinal cord expression in the V0v domain of the spinal cord at 20 h and 24 h (cf. Figure  . By 20 h, neff1 is expressed in a few dorsal spinal cord cells (Fig. 2AU), which is very similar to both nefma ( Fig. 2AK) and nefmb, although nefmb is also expressed in some ventral spinal cord cells (Fig. 2AP). By 24h, nefma, nefmb and neff1 are expressed in the V0v spinal cord domain, although neff1 is still only expressed by a few cells (cf. Figure 2AL, AQ, AV to C). All three of these genes continue to be expressed in this domain at 36 h and 48 h, although at these stages they are clearly expressed by more cells than evx1 (cf. Figure 2AM-AN, AR-AS, AW-AX to D-E), suggesting that they are expressed in additional spinal cell types. In contrast to these other NIF genes, inab is expressed in a similar dorsal-ventral region of the spinal cord to evx1 at all the developmental stages that we analyzed (cf. Figure 2AY-AAC to A-E). However, from 36 h onwards, inab is very broadly expressed in the spinal cord suggesting that while its earlier expression may be more speci c, at 36 and 48 h it is expressed by the majority of post-mitotic spinal cells (cf. Figure 2AAB-AAC to D-E). The temporal expression pattern of inab in the V0v domain is, therefore, similar to that of skor1b, ebf3a and uncx. In contrast, the delayed onset of neff1, nefma and nefmb expression in the V0v domain, some 7 h later, is similar to that of skor2.
We also performed uorescence in situ hybridization, for a subset of the genes that had the strongest expression in the single in situ hybridization experiments, in Tg(evx1:EGFP) SU1 embryos, in which EGFP spinal cord expression is exclusively in V0v cells (14). These double-staining data con rm that skor1b, skor2, uncx, and nefma are expressed by V0v interneurons (double-labelled cells) as well as non-V0v spinal cells (red but not green cells; Fig. 3 and also see (14) for different complementary data for skor2).  Table 1). We were unable to reliably count spinal cells expressing uncx4.1, because the expression was so weak and punctate, even after using the molecular crowding reagent Dextran Sulfate (see Methods), and prolonged staining. However, we did not observe any differences in the spinal cord expression of this gene between WT and evx1;evx2 double mutant embryos (Fig. 4P-Q).  Table 1). In contrast, we did not nd statistically signi cant differences in the number of cells expressing nefmb or inab ( Fig. 5D-F, J-L, Table 1).
As we had identi ed three NIF genes in our microarray analyses (nefma, neff1, and inab), for completeness, we decided to also examine the expression of the remaining two zebra sh NIF genes, ne a and ne b. Our microarray data suggested that these were both expressed in the spinal cord but not in V0v interneurons (Supp. Figure 2A). Consistent with this, we found no change in the number of cells expressing either ne a or ne b in evx1;evx2 double mutants compared to WT siblings (Supp. Figure 2B-G, Table 1).
While we would expect Evx1 and Evx2 to act cell autonomously, as they are both transcription factors, it is still important to con rm that the spinal cells that lose expression of particular genes in evx1;evx2 double mutants are, indeed, V0v interneurons. Therefore, we FAC-sorted EGFP-positive V0v spinal interneurons from the progeny of an incross of heterozygous evx1;evx2 double mutant sh that were homozygous for Tg(evx1:EGFP) SU2 and performed single cell RNA sequencing (scRNA-seq). We cannot distinguish evx1;evx2 single or double mutant embryos from WT siblings morphologically, so our V0v interneurons were collected from all of the different genotypes generated from this cross. However, we were able to distinguish mutant cells from WT cells in our data analyses, based on each cell's individual gene expression pro le (see Methods). We performed these analyses at 48 h because EGFP is expressed in signi cantly more V0v interneurons in 48 h Tg(evx1:EGFP) SU2 embryos than at earlier developmental stages (14), enabling us to capture a larger number of rare double mutant cells. (Only 1/16 embryos from an incross of evx1;evx2 heterozygous parents will be double mutants, and, in order to maintain mRNA integrity, we can only dissect trunks for a limited amount of time for each experiment. See Methods for more details).
Following FAC-sorting and scRNA-seq, we generated a dataset of 2860 V0v cells that passed stringent quality controls (see Methods).. For improved visualization and interpretation of this single-cell atlas, we used Uniform Manifold Approximation and Projection (UMAP) plotting, since this preserves the global structure of the expression data. We manually inspected UMAP plots to assess clustering quality based on expression of known V0v spinal interneuron markers. Using these methods, we identi ed ve distinct clusters of V0v interneurons (Fig. 6A). We visually compared the expression of different genes in these clusters and performed gene-speci c analyses of differential expression using both Hurdle model and ANOVA statistical comparisons (see Methods for an explanation of why we used these statistical methods). Comparing the expression pro les of evx1, evx2 and neurotransmitter phenotype genes in these ve distinct clusters with our previously published data (14,72), suggested that the two clusters with the highest expression levels of evx1, evx2 and the excitatory (glutamatergic) transmembrane transporter gene slc17a6a, and the lowest expression levels of inhibitory neurotransmitter genes slc6a5, slc6a1b and gad1b, are WT V0v cells ( Fig. 6B-G, Table 2). In contrast, the clusters with lower expression levels of evx1, evx2 and slc17a6a, and higher expression levels of slc6a5, slc6a1b and gad1b, are most likely to be mutant V0v cells (Fig. 6B-G, Table 2). Mutant group 3 is the smallest cluster of pro led cells, and the cluster with the most severe reduction in evx1, evx2 and slc17a6a expression ( Fig. 6A-D, Table 2) and the highest level of expression of the inhibitory (glycinergic) gene slc6a5 (Fig. 6E, Table 2). This suggests that it contains double mutant cells, as we would expect these cells to have the most severe phenotype. Gene-speci c analyses of differential expression, created through (A) Hurdle model and (B) ANOVA statistical comparisons between distinct cell clusters in our 48 h evx1 i232/+ ;evx2 sa140/+ heterozygote incross single-cell atlas (see Fig 6A and also Methods for experimental details and rationale for using both statistical methods). The Hurdle model is generally the most statistically robust method for these analyses, if there are enough cells in each group. In contrast, ANOVA usually performs better when the numbers of cells being compared are very small (see Methods for more information). We are providing the data for both methods for completeness and comparison. Column 1 shows the gene symbol.  Our previously published data suggest that the phenotypes of evx1 and evx2 single mutants and evx1;evx2 double mutants differ only in their severity / penetrance. For example, both single mutants lose expression of evx2, slc17a6 and skor2 in some spinal cord cells, whereas in the double mutants, more cells lose expression of these genes, and we also see a statistically signi cant reduction in the number of spinal cord cells expressing evx1 (14). (Due to the molecular nature of the evx1 and evx2 mutations, we can still detect mRNA for both these mutated genes). Given the similarity of the single and double mutant phenotypes, before we performed these scRNA-seq analyses, we were not sure whether all mutant cells would cluster together, or whether we would see distinct clusters of less severe and more severe mutant cells. It is theoretically possible that what distinguishes the mutant group 1 cluster from the mutant group 2 cluster is the severity of the mutant phenotype. However, the data appear to be inconsistent with this hypothesis. We do not see a statistically signi cant difference in evx1 expression between mutant group 1 and mutant group 2, and while evx2 is expressed at slightly higher levels in mutant group 2, this difference is only statistically signi cant when we use the Hurdle model ( Table 2). The most striking difference between the mutant group 1 and mutant group 2 clusters is the much higher, statistically signi cant expression of inhibitory marker slc6a1b in mutant group 1 cells (Table 2). It is not clear though what the functional signi cance of this difference is.
Statistically, mutant group 1 cells also have more expression of a different inhibitory marker gad1, than mutant group 2 cells, but statistically they also have less expression of the inhibitory marker slc6a5, and more expression of the excitatory marker slc17a6a (Table 2). Based on these data we think that it is more likely that, as discussed below, the mutant group 1 and 2 clusters are mutant versions of the two molecularly distinct WT clusters that we have identi ed.  (14), the latter explanation would not be surprising. Consistent with this, we detected expression of the inhibitory marker slc6a5 in a few cells in both WT clusters (Fig. 6E). slc6a1b is also detected in a very small number of both WT cell types, although interestingly, gad1b is not (Fig, 6F-G). Three-way differential gene expression shows that a few of the cells in WT groups 1 and 2 that express slc6a5 also express either evx1, evx2 (white/grey cells, Supp. Figure 3C) or the excitatory marker slc17a6a (turquoise cells, Supp. Figure 3D), suggesting that they may be cells with a partial mutant phenotype.
Similarly  Table 2), suggesting that these genes may be expressed by fewer V0v interneurons at later stages of embryogenesis, or that they may be expressed at low levels, and hence drop out of the pro les of some cells. nefmb was also only detected in a small number of cells in the two WT groups, although it was expressed by more cells in mutant groups 1 and 3 (Fig. 6O, Table 2). In contrast, skor2, ebf3a, uncx, nefma, neff1, and inab ( Fig. 6J-L, N, P-Q, Table 2) expression was detected in many V0v cells. Interestingly, we detected skor1b, skor2, ebf3a, uncx, uncx4.1, nefma, nefmb, and neff1 expression in more cells in the UMAP plots, and at statistically-signi cantly higher levels in our gene-speci c analyses of differential expression, in WT group 1 compared to WT group 2 ( Fig. 6I-P, Table 2A).
Using the Hurdle model (as both WT clusters contain a relatively large number of cells, and so the variance can be effectively modelled (please see Methods for further explanation)) this increase is most dramatic for skor2, but also substantial for neff1 and ebf3a. The other increases are more modest, although still statistically signi cant. In contrast, skor1a and inab were expressed at slightly higher, statistically signi cant, levels in WT group 2 cells than in WT group 1 cells (Fig. 6H, Q, Table 2A).
Consistent with our in situ hybridization data at 30 h (Figs. 4-5, Table 1), our 48 h scRNA-seq data suggest that expression of skor1a, skor1b, skor2, ebf3a and neff1 in V0v interneurons requires Evx1/2. We detected fewer V0v cells expressing skor1a, skor1b, skor2, ebf3a, and neff1 in the three mutant clusters, with the exception that neff1 is still expressed in many mutant group 3 cells (Fig. 6H-K, P). These results are also con rmed by our gene-speci c analyses of differential expression, although no statistical analyses could be performed for some of the skor1b comparisons as no cells expressing it were detected in mutant groups 1 or 2. In addition, the difference between mutant group 1 and WT group 1 is not statistically signi cant for skor1a, using the Hurdle model (Table 2A). In general, statistical analyses using both Hurdle and ANOVA models give similar results, although there are some subtle differences and in most of these comparisons, the Hurdle model is probably more reliable because of the increased e ciency of variance modelling over ANOVA when there are large cell numbers in the groups under comparison (see Methods).
In contrast, for uncx, nefma, nefmb, and inab, we identi ed differences between our in situ hybridization data at 30 h (Figs. 4 and 5, Table 1) and our 48 h scRNA-seq data (Fig. 6L, N-O, Q, Table 2). For example, while we detected no change in the number of cells expressing uncx, in evx1;evx2 double mutant embryos compared to WT siblings in our 30 h in situ hybridization experiment (Fig. 4M-O, Table 1), in our 48 h scRNA-seq data, we found that uncx is expressed in a higher proportion of cells, and at statistically-signi cant higher levels, in mutant groups 1 and 2 than the two WT groups (Fig. 6L, Table 2). In contrast, there are hardly any cells expressing uncx in mutant group 3 (Fig. 6L, Table 2). Similarly, at 30 h ( Fig. 5A-F, Table 1), we detected a slight, but statistically signi cant, reduction in the number of spinal cells expressing nefma in evx1;evx2 double mutants, but no change in the number of spinal cells expressing nefmb, whereas in our 48 h scRNA-seq experiment, we detected a higher proportion of nefmaand nefmb-expressing cells, and a statistically signi cant higher level of expression of these genes in all three mutant clusters (Fig. 6N-O and Table 2, the comparison between mutant group 2 and WT group 2 is not statistically signi cant for nefmb). In addition, at 30 h, we did not nd a difference in the number of spinal cord cells expressing inab in evx1;evx2 double mutants ( Fig. 5J-L, Table 1), whereas in our 48 h data, inab is expressed in a lower percentage of cells in the mutant clusters than in the WT clusters (Fig. 6Q, as this result is harder to see from the UMAP plot alone, we con rmed it by quantifying the number of inab-expressing cells per group in a dotplot showing the number of inab reads detected per cell per group (Table 3)) and there is also a statistically signi cant reduction in its expression in the mutant groups compared to the WT groups (Table 2). Our data for uncx4.1 is less conclusive than for these other genes as we were not able to count the number of spinal cord cells expressing this gene at 30 h. However, we didn't see any obvious change in its expression at this stage, which is consistent with our 48 h data, where the only statistically signi cant change that we see in mutant cells, is for mutant group 2 cells compared to WT group 2 cells. Our scRNA-seq data suggests that very few V0v interneurons express uncx4.1, in any of the different clusters at 48 h. Although, based on the UMAP plots and differential gene expression analyses, there is statistically slightly more expression in WT group 1 than WT group 2. Taken together, these data suggest that Evx1/2 may repress uncx4.1 expression slightly in WT group 2.  Third row shows the equivalent data for en1b.  (1) 28.36% (57) scRNA-seq analysis identi es two distinct V0v WT clusters, each of which is most similar to a different mutant cluster.
Our analyses of differential gene expression in the ve distinct clusters of V0v cells, suggest that mutant group 1 is more similar to WT group 1 than WT group 2 and, conversely, mutant group 2 is more similar to WT group 2 than WT group 1 (Figs. 6 & 7, Tables 2 & 4). For example, there is a statistically signi cant increase in expression of nefma and nefmb in WT and mutant group 1 cells compared to WT and mutant group 2 cells (Fig. 6A, N-O; Table 2). Many other genes, including neff1, anos1a, chrna2b, fndc4b, plpp4, cnih3, and drd2b are also expressed at higher levels in WT and mutant group 1 cells compared to WT and mutant group 2 cells (Fig. 6P, Fig. 7A-G, Tables 2 & 4, Supp. Tables 2 &   3). Similarly, several genes, including esrrb, scxa, and svild, are expressed at statistically signi cant higher levels in WT and mutant group 2 cells, compared to WT and mutant group 1 cells (Fig. 7A, H-J, Tables 2 & 4, Supp. Tables 2 &   3). Taken together, these data suggest that there are two molecularly distinct subsets of WT V0v cells and a mutant version of each. Gene-speci c analyses of differential expression, created through Hurdle model (columns 3-4) and ANOVA (columns 5-6) statistical comparisons between distinct cell clusters in our 48 h evx1 i232/+ ;evx2 sa140/+ heterozygote incross single-cell atlas (see Fig 7A and also Methods for experimental details and rationale for using both statistical methods). For these comparisons, the Hurdle model is probably the most statistically robust method as there are su cient cell numbers in each group to effectively model the variance (see Methods for more information). We also provide the ANOVA data for completeness and comparison.  Cells in mutant group 3 appear to represent the most severely affected evx1/2 mutant cells, based on several key criteria. These cells have the lowest levels of evx1 and evx2 expression, almost none of them express the glutamatergic gene, slc17a6a, most of them express the inhibitory glycinergic gene, slc6a5, and some cells also express inhibitory GABAergic markers, including slc6a1b and gad1b ( Fig. 6A-G, Table 2). To our surprise, differential expression analyses between cells in this cluster and all the other clusters, identi ed several transcription factor genes that are usually expressed by either distinct populations of inhibitory spinal cord interneurons, or motoneurons, upregulated in mutant group 3 cells (Fig. 8, Table 5, Supp. Tables 2 & 3). For example, there was a statistically signi cant increase in the expression of gata2a, gata3, tal1 and sst1.1 (expressed by KA and V2b cells (37,71,73,74)), en1b (expressed by V1 cells (6, 70)), dmrt3a (expressed by dI6 cells (56, 75) (76)), lbx1a (expressed by dI4 and dI6 cells, (4, 7, 77-80)) ( Fig. 8A-G, Table 5) and isl1a, isl2a, isl2b, mnx1, mnx2a, and mnx2b (expressed by motoneurons (81-84) (Fig. 8H-M, Table 5)). Interestingly, UMAP analysis suggests that these markers of different cell types are expressed by small distinct groups of mutant group 3 cells. For example, the mutant group 3 cells that express the KA and V2b genes gata2a, gata3, tal1, or sst1.1 overlap with each other (Fig. 8A-D), but these cells are spatially distinct from those that express markers of motoneurons, V1 interneurons or more dorsal inhibitory interneurons (cf. Figure 8E-M and see 8S-X). Similarly, motoneuron genes isl1a, isl2a, isl2b, mnx1, mnx2a, and mnx2b are expressed by a group of overlapping mutant group 3 cells, and these cells do not co-express interneuron markers (cf. Figure 8H-  Gene-speci c analyses of differential expression, created through Hurdle model (columns 4-6) and ANOVA (columns 7-9) statistical comparisons between distinct cell clusters in our 48 h evx1 i232/+ ;evx2 sa140/+ heterozygote incross single-cell atlas (see Fig 6A and also Methods for experimental details and rationale for using both statistical methods). For these comparisons, the ANOVA data is probably the most robust, as the number of cells in each subcluster is relatively small (see Methods for more information). We also include the Hurdle model data for completeness. Column 1 shows the gene symbol. Column 2 indicates the spinal neuron types that normally express the gene. Column 3 indicates the neurotransmitter phenotype for these neurons. Columns 4 and 7 show least squares mean read counts for cells in mutant group 3, and columns 5 and 8 show least squares mean read counts for cells in all the other clusters combined, respectively. Columns 6 and 9 show fold-change values. = fold-change increase, ¯ = fold-change decrease in mutant group 3 (the antecedent ( rst) population) compared to all other populations combined (the consequent (second) population). Statistically signi cant (P<0.05) values are indicated in bold. *** P<0.001. ** P<0.01. N.C. = Not Calculated. Hurdle model of differential expression analysis cannot be calculated, usually because expression was too low or not present in one of the groups being compared. Additional data for this comparison is available in Supp. Data Tables 2 (Hurdle model data) and 3 (ANOVA data).  In contrast, with the exception that there may be an occasional cell at the boundary between mutant group 3 and WT group 1 that expresses slc17a6a and either vsx2 or tlx3b (cf. Figure 6D & Fig. 8O-P), we do not observe upregulation of genes expressed by spinal excitatory interneuron populations in the UMAP analysis. With the exception of tlx3b, there is also no statistically signi cant change in the expression of genes expressed by spinal excitatory interneuron populations, in mutant group 3 cells compared to the other clusters (Table 5). Speci cally, we examined sim1a (V3 interneurons (85, 86)), vsx2 (V2a interneurons, (37, 87)), tlx3b (dI3 and dI5 interneurons (4, 11, 88)), foxp2 (dI2 interneurons (89, 90)) and barhl2 (dI1 interneurons (91)) expression ( Fig. 8N-R, Table 5). Taken together, these data suggest that distinct subsets of mutant group 3 cells have started to express markers of distinct inhibitory spinal neurons (V1, V2b, KA, dI4 or dI6 neurons), or cholinergic motoneurons.
Additional transcription factor genes are either downregulated or upregulated in V0v interneurons that lack Evx1 and/or Evx2 function.
In addition to the genes that we had already identi ed, our scRNA-seq data also identi ed additional transcription factor genes that may be part of V0v GRNs, as they are expressed in at least some V0v WT cells and are downregulated in V0v mutant cells (Fig. 9, Table 6). ccdc3a, dachc, luzp1, mycb, nr5a2, pou3f1, pou3f2b, pou3f3b and scrt2 are expressed in both WT clusters and this expression is reduced in all three mutant groups (Fig. 9B-J and Table  6A, with the exceptions that the reduction of expression in the mutant group 3 cluster is not statistically signi cant using the Hurdle model for dachc and luzp1, and expression of mycb, and possibly scrt2, is upregulated in the mutant group 3 cluster compared to WT cells). There is also a statistically signi cant reduction of pou2f2a, pou2f2b, mafba, pbx1b, scrt1a and zfhx3b expression in all three mutant groups, according to the Hurdle model. These genes are all expressed by many cells in WT and mutant group 2, as well as a signi cant number of cells in the WT group 1 cluster and some cells in the other two mutant groups ( Fig. 9K-P, Table 6A). Interestingly, pou2f2a, pou2f2b and zfhx3b appear to be co-expressed in a subset of mutant group 3 cells at the top of the cluster, which also express makers of V1 or dI6 cells (cf. Figure 9K-L& P to Fig. 8V, X). nr2f5 has a similar expression pattern to these six genes, except that far fewer cells express this gene in the WT or mutant group 1 clusters (Fig. 9Q). In contrast, ebf1a and pitx2 are mainly expressed by cells in the WT group 1 cluster and there are very few cells that express either of these genes in either the WT group 2 cluster or any of the mutant clusters ( Fig. 9R-S, Table 6). Gene-speci c analyses of differential expression, created through (A) Hurdle model and (B) ANOVA statistical comparisons between distinct cell clusters in our 48 h evx1 i232/+ ;evx2 sa140/+ heterozygote incross single-cell atlas (see Fig 9A and also Methods for experimental details and rationale for using both statistical methods). For these comparisons, the Hurdle model is probably the most statistically robust method as there are su cient cell numbers in each group to effectively model the variance (see Methods for more information). We also provide the ANOVA data for completeness and comparison. Column 1 shows the gene symbol.  We also identi ed some transcription factor genes that are upregulated in V0v mutant cells compared to WT cells (Fig.   10, Table 7). As with the downregulated genes, these have a few different patterns of expression. For example, hmx2, hmx3a, otpb, and znf385c are all highly expressed in mutant group 1 and 2 cells, but only a few cells express these genes in either of the WT clusters or, in the case of hmx2, hmx3a and otpb, in the mutant group 3 cluster. znf385c is expressed in more mutant group 3 cells than the other three genes and it is expressed at statistically signi cant higher levels in this cluster than in the other two mutant clusters (Fig. 10B-E, Table 7). znf385a has a similar pattern of expression to these four genes, but it is expressed in fewer mutant group 1 and 2 cells, and like znf385c, it is expressed at highest levels in the mutant group 3 cluster (Fig. 10F, Table 7). In contrast, zmat4b is predominantly expressed by cells in mutant group 1, although there is still a statistically signi cant increase in expression of this gene in the mutant group 2 cluster compared to the WT group 2 cluster and in the mutant group 3 cluster compared to WT cells (Fig. 10G, Table 7). In contrast, bhlhe22 and irx1a are predominantly expressed by cells in the mutant group 2 cluster, although there is still a statistically signi cant increase in expression of bhlhe22 in the mutant group 1 cluster compared to WT group 1, and of bhlhe22 and irx1a in the mutant group 3 cluster compared to WT cells (Fig. 10H-I, Table 7). Gene-speci c analyses of differential expression, created through (  V0v interneurons ectopically express hmx3a in evx1;evx2 mutants. Two of the genes that were upregulated in the majority of V0v mutant group 1 and mutant group 2 cells were hmx2 and hmx3a (Fig. 10B-C, Table 7). We have previously shown that Hmx2 and Hmx3a are co-expressed in dI2 and V1 spinal interneurons, and that Hmx3a is required for the excitatory fates of dI2 interneurons. (Hmx2 also has a role in this process, but it is much more subtle than that of Hmx3a (12)). In the absence of hmx3a function, many dI2 interneurons change their neurotransmitter fates from glutamatergic (excitatory) to GABAergic (inhibitory) (12). Therefore, we were surprised to discover that hmx3a is upregulated in mutant V0v cells that have changed their neurotransmitter phenotype from excitatory to inhibitory.
To further con rm this intriguing result, and examine whether it is also the case at earlier developmental stages, we analyzed hmx3a expression in 27 h evx1;evx2 double mutants. We found that there is a statistically signi cant increase in the number of hmx3a-expressing cells in the double mutants, compared to WT siblings ( Fig. 11A-C, Table   1). Double-labelling experiments with Tg(evx1:EGFP) SU2 con rmed that while V0v interneurons do not express hmx3a in WT embryos, many of them turn on hmx3a expression in evx1;evx2 double mutants (Fig. 11D-E'''), suggesting that Evx1/2 normally repress hmx3a expression in V0v interneurons. To assess whether Hmx3a might reciprocally repress evx1 and evx2 expression in V1 and dI2 spinal interneurons, we examined expression of evx1 and evx2 in a deletion mutant that lacks both hmx2 and hmx3a. However, we found no change in the number of spinal cord cells expressing either evx1 or evx2 in hmx2;hmx3a deletion mutants compared to WT siblings ( Fig. 12B-G Figure 4)). To our surprise, we found that expression of the V0v gene, skor1a, was upregulated more than 13-fold in dI2 and/or V1 interneurons from hmx2;hmx3a double morphants compared to uninjected controls (Fig. 12A, Supp. Table 4). This suggests that Hmx2 and/or Hmx3a normally repress skor1a expression in dI2 and/or V1 cells. Intrigued by this result, we examined whether expression of any of the other ve genes downregulated in evx1;evx2 mutants at 30 h, skor1b, skor2, ebf3a, nefma or neff1 (Fig. 4D-L, Fig. 5A-C, G-I, Table 1) was upregulated in hmx2;hmx3a double morphants. The bulk RNA-seq data suggested that skor1b and nefma expression might also be upregulated (Fig. 12A, Supp. Table 4). Too few reads were reliably detected to assess whether skor2 is also upregulated and the fold-change for both ebf3a and neff1 was less than two, suggesting that any differences might be due to noise in the experiment.
As some of these results were inconclusive, and also because our previous analyses identi ed some differences between the phenotypes of hmx2;hmx3a double morphants and homozygous hmx2;hmx3a deletion mutants (12), we decided to analyze the spinal cord expression of skor1a, skor1b, skor2, ebf3a, nefma, and neff1 in these mutants. We found that there was a statistically signi cant increase in the number of spinal cord cells expressing either skor1a or nefma in hmx2;hmx3a deletion mutants compared to WT siblings ( Fig. 12H-J, T-V, Table 1). However, there was no statistically signi cant difference in the number of spinal cord cells expressing skor1b, skor2, ebf3a or neff1 (Fig. 12K-S, W-Y, Table 1). As discussed above, hmx2 and hmx3a are expressed in both dI2 and V1 cells. Therefore, to test which of these two cells types was upregulating expression of skor1a and nefma, we performed double-labelling experiments with Tg(pax2a:GFP), which, in the spinal cord, is exclusively expressed in V1 cells (6). We observed no coexpression of GFP and either skor1a or nefma in either WT siblings or hmx2;hmx3a deletion mutant embryos, suggesting that the cells that ectopically express skor1a and nefma in hmx2;hmx3a deletion mutants are dI2 interneurons and not V1 interneurons (Fig. 12Z-AC'"). This suggests that Hmx2/3a normally repress skor1a and nefma expression in dI2 interneurons. Taken together with our V0v data described earlier in this paper, these data raise the intriguing possibility that Evx1/2 might regulate expression of skor1a and nefma in V0v interneurons by repressing expression of Hmx3a, which itself normally represses skor1a and nefma.  Tables 1 & 2). Therefore, these ve genes are all potential members of GRNs downstream of Evx1/2 in V0v interneurons. It is also likely, given that each of these genes has distinct patterns of spinal expression in addition to the V0v domain, that they also have roles in different subsets of other spinal neurons (Fig. 13A & C).
Consistent with this hypothesis, skor2 and skor1a, unlike skor1b, are both expressed in the dorsal spinal cord at 17 h, but they have very different expression patterns to each other at 20 h (Fig. 2F-G, P-Q, Fig. 13A), ebf3a is expressed in cells ventral and dorsal to the V0v domain at all stages examined (Fig. 2U-Y, Fig. 13A) and neff1 is expressed in a broad dorsal-ventral stripe of cells at 36 h and 48 h (Fig. 2AW-AX, Fig. 13A). In addition, while the number of spinal cells that express each of these genes is reduced in evx1;evx2 double mutants, in each case some cells remain  Table 1). This suggests that fewer V0v interneurons express skor1a and skor2 than skor1b at this stage. This could be because the former two genes are only expressed in distinct subsets of V0v interneurons, or it could re ect the later expression of these genes in V0v cells. i.e., maybe at this stage of development these genes have only turned on in the "older" V0v interneurons and not yet in the "younger" ones. On the other hand, the reduction in the number of spinal cells expressing ebf3a in 30 h evx1;evx2 mutants (25 cells) is very similar to that for skor1a and skor2, even though our in situ hybridization experiments suggest that ebf3a is expressed in V0v interneurons as early as skor1b (Fig. 2K, U, Fig. 4J-L, Fig. 13A, Table 1). Interestingly, while our 30 h data suggest that more V0v interneurons express skor1b than skor1a or skor2, at 48 h, skor2 is expressed in many more V0v interneurons than skor1a or skor1b, and, while there is not much difference between skor1a expression in the two different WT groups, skor1b and skor2 are expressed by more WT group 1 cells than WT group 2 cells (Fig. 6H-J, Table 2). This suggests that skor gene expression in V0v interneurons is dynamic, that speci c subsets of V0v interneurons express distinct combinations of skor genes at different times in development, and that skor1a and skor1b may only be expressed by V0v interneurons for a relatively short period.
In addition to being expressed in different subsets of other spinal cord cells, these genes also differ in when they are rst expressed in the V0v spinal domain. For example, skor1b is expressed in the V0v spinal region as early as 17 h (Fig. 2K, Fig. 13A), whereas skor1a is not expressed in this region until 20 h and skor2 not until 24 h (Fig. 2G, R,  Fig. 13A). ebf3a is, like skor1b, expressed in the V0v region as early as 17 h (Fig. 2U, Fig. 13A) whereas neff1 is not expressed in this spinal cord domain until 24 h (Fig. 2AV, Fig. 13A). This is important because temporal differences in expression may re ect different positions in the hierarchy of interactions that make up GRNs, with earlier expressed genes regulating the expression of later expressed genes (Model II, Fig. 13B and cf. Model I). This could be tested in future work by examining epistatic relationships between these genes. Alternatively, later-expressed genes may require a later-expressed transcription factor in addition to Evx1/2, in order to be expressed (Model III, Fig. 13B and cf. Model II). Interestingly, our data suggest that there is limited overlap between the varied spatial and temporal expression patterns of these genes. This in consistent with a model where different types of spinal interneurons are speci ed by distinct GRNs, rather than there being a cassette of genes that speci es the same functional characteristic in different neurons.
Currently, there is very little known about the functions of any of these genes in spinal cord development. neff1 encodes a NIF protein (also referred to as Type IV Intermediate Filament proteins). In mature mammalian CNS, NIF proteins are important for axon function and maintenance, but their functions during development are less clear (92)(93)(94)(95)(96). (See also more detailed discussion of NIF proteins below). Ebf and Skor proteins have been implicated in development of different neurons in the brain (e.g. 62, 64, 66, 97-100), but very little is known about their spinal cord functions in any vertebrate. Mouse Skor1 is expressed in dI4 and dI5 spinal interneurons, where it binds Lbx1 (101), and human Skor1 has been implicated in Restless legs syndrome (also known as Willis-Ekbom disease) (102). Ebf3 is expressed in spinal interneurons in mouse and chick (99,103) and correct Ebf3 spinal expression requires Lmx1bin mouse (65). This is interesting, given that Lmx1ba/b are also downstream of Evx1/2 in zebra sh V0v interneurons (13). However, none of these previous data indicate what the functions of these transcription factor genes are in V0v spinal interneurons.
Our data, as discussed above, show that all these genes are downstream of Evx1/2 in V0v interneurons. To date, the only abnormal V0v interneuron phenotype that we have found in zebra sh evx1;evx2 single or double mutants is their change from being glutamatergic to glycinergic (14). This suggests that these genes may encode members of GRNs that regulate this phenotype. However, it is also possible that some of these genes function downstream of Evx1/2 in other aspects of V0v differentiation and/or function that we have not yet detected. So far, most of the identi ed transcription factor genes that specify spinal neuron functional characteristics are expressed throughout the stages of development that we examined by in situ hybridization (e.g. 3,5,13,14). However, skor2 and neff1 are not expressed in the V0v spinal cord domain until 24 h (Fig. 2R, AV, Fig. 13A), which is probably too late to be required for specifying V0v cell glutamatergic fates (14). These genes may, instead, be involved either in maintaining correct neurotransmitter phenotypes, or in specifying later aspects of V0v interneuron development. To test these different hypotheses, future studies will need to analyze the phenotypes of V0v interneurons in mutants for each of these genes. Evx1/2 may regulate the expression of uncx, nefma, nefmb, andinabin different ways at different developmental time points.
Our data also suggest that uncx, uncx4.1, nefma, nefmb, and inab are expressed by V0v interneurons. uncx and inab are expressed in the V0v spinal cord domain at all the stages that we examined by in situ hybridization as well as being expressed by cells in both scRNA-seq WT clusters at 48 h (Figs. 2Z-AD, AY-AAC, 6L, Q, Fig. 13C, Table 2). nefma and nefmb are also expressed in the V0v domain at 24 h and at older stages ( Fig. 2AL-AN, AQ-AS, Fig. 13C), although nefma and nefmb are expressed by more WT group 1 cells that WT group 2 cells in our scRNA-seq data at 48 h and nefmb is only expressed by a few WT cells in this data set ( Fig. 6N-O, Table 2). In contrast, uncx4.1 is only expressed transiently in the V0v spinal cord domain at 20 h and 24 h in our in situ hybridization data, although we do detect a very small number of V0v cells expressing this gene in our 48 h scRNA-seq data (Fig. 2AF-AG, Fig. 6M, Fig. 13C).
It is less clear whether any of these genes require Evx1/2 function for their expression in V0v interneurons. We did not see any obvious difference in uncx4.1 expression in 30 h evx1;evx2 double mutants (Fig. 4P-Q) and the only difference in the mutant clusters in the 48 h scRNA-Seq data was a slight, but statistically-signi cant increase in expression in the mutant group 2 cluster compared to the WT group 2 cluster (Fig. 6M, Table 2). In addition, while the phenotypes of skor1a, skor1b, skor2, ebf3a, and neff1 spinal expression in the absence of Evx1/2 function were generally consistent between our experiments at 30 h and 48 h (the one exception being the expansion of neff1 expression in mutant group 3 at 48 h, Fig. 6P, Table 2), this was not the case for uncx, nefma, nefmb, and inab. While we saw no statistically signi cant change in the number of cells expressing uncx, nefmb or inab in 30 h evx1;evx2 mutants (Figs. 4-5, Table 1), at 48 h we see an increased number of cells expressing uncx and nefmb in mutant group 1 and 2 clusters (although the increase is not statistically signi cant for nefmb in mutant group 2) and, in contrast, the expression of inab is statistically-signi cantly reduced in mutant clusters (Fig. 6L, O, Q, Table 2). These differences suggest that the expression of these four genes is regulated differently at these two developmental stages. For example, maybe inab only requires Evx1/2 for its expression in V0v interneurons at later stages, and uncx and nefmb require Evx1/2 to be turned off in V0v interneurons at later stages. However, this does not explain the more surprising difference in nefma expression, where we saw a statistically-signi cant reduction in the number of cells expressing nefma in 30 h evx1;evx2 mutants ( Fig. 5A-C, Table 1), but we see a statistically-signi cant increase in its expression in all three mutant clusters at 48 h (Fig. 6N, Table 2A). While it is possible for transcription factors to function as both activators and repressors of transcription (e.g. 104,[105][106][107][108][109][110][111][112], it is unusual that transcription factors like Evx1 and Evx2 would be required both to turn a speci c gene on, and then to later turn it off, in the same cells. uncx and uncx4.1 encode paired-type homeodomain transcription factors and while spinal cord expression of uncx in zebra sh and mouse, and uncx4.1 in mouse has previously been reported (113)(114)(115)(116), there is currently no data available on the function of these genes in spinal interneurons. nefma, nefmb, and inab are, like neff1, NIF genes. These three genes have distinct expression patterns from each other and from neff1. They are all expressed in the V0v spinal cord domain from 24 h onwards, but only inab is expressed in this domain at the earlier stages that we examined. All four of these genes are also expressed in other spinal cord domains at all the stages that we analyzed ( Fig. 2AJ-AAC, Fig. 13A, C). These expression patterns suggest that all of these genes function not just in V0v interneurons, but also in other distinct subsets of spinal neurons.
As mentioned above, NIF proteins are important for axon function and maintenance at later stages of differentiation in mammals, but a developmental function for these proteins has not yet been described (92)(93)(94)(95)(96). Our data suggest that expression of nefma, nefmb, inab, and neff1 genes is regulated by Evx1/2 in V0v INs, albeit in different ways. The NIF proteins encoded by these genes have DNA-binding domains (117,118), and it has been suggested that these domains may regulate gene expression (119,120). While these proteins are generally thought to be cytoplasmic, it is possible that either full-length proteins or shorter forms of the proteins enter the nucleus (121,122). For example, the N-terminal DNA-binding domain of Vimentin, a different intermediate lament protein, can enter the nucleus and regulate nuclear architecture and chromatin structure (122). Therefore, it is possible that neff1 and nefma, both of which are expressed by statistically-signi cantly fewer cells in 30 h evx1;evx2 mutants, are part of the GRN that regulates neurotransmitter properties in V0v interneurons ( Fig. 5A-C, G-I, Table 1). If this is not the case, then they likely function downstream of Evx1/2 in a not-yet-identi ed aspect of V0v interneuron development. However, evx1;evx2 double mutants have no obvious axon defects during development (we have examined stages up to 48 h; 14), suggesting that these NIF proteins are not required for axon outgrowth or path nding at the stages that we are examining. It is even less clear what the NIF genes that are upregulated in mutant V0v interneurons are doing. Either way, our data suggest that there are novel developmental functions for these genes.
Additional candidate GRN transcription factor genes downstream of Evx1 and Evx2 in V0v interneurons.
Interestingly, in mice, a small subgroup of V0v interneurons (less than 10%) express Pitx2. This compares to 17.45% of the WT cells in our 48 h data (312/933 WT group 1 cells and 12/924 WT group 2 cells). In mouse, the Pitx2 + V0v interneurons are preferentially clustered around the central canal and can be further subdivided into excitatory cholinergic (V0c) and glutamatergic (V0g) subtypes (134). Pitx2 has also been shown to be necessary and su cient to drive expression of the inhibitory neurotransmitter gene GAD1 in C. elegans type D GABAergic motor neurons and it is required for the GABAergic differentiation of superior colliculus cells in the mouse brain (135,136). In contrast, Pitx2 function in V0v interneurons is still unknown, but as these cells are excitatory, it may be distinct from its role in these other cell types.
Currently, there is not a lot of data on spinal cord functions of most of the other genes. One exception is scrt2. When scrt2 was knocked down with morpholinos in zebra sh embryos, the number of islet2positive motoneurons was increased compared to uninjected controls, whereas the number of glutamatergic Rohon-Beard sensory neurons, and pax2a-positive inhibitory interneurons were unchanged (other types of interneurons were not examined) (130). This is intriguing, given that some of the mutant group 3 cells ectopically express motoneuron markers.
We also identi ed eight transcription factor genes that are upregulated in mutant V0v interneurons at 48 h. hmx2, hmx3a, otpb, znf385c, znf385a, zmat4b, bhlhe22 and irx1a were all detected in a small number of WT V0v interneurons, and many more mutant cells (Fig. 10, Table 7). As far as we are aware, this is the rst report of spinal cord expression of znf385a, znf385c or zmat4b, and the rst report of irx1a expression in the zebra sh spinal cord (Irx1 is expressed in the mouse spinal cord next to the hind limb (137)). However, interestingly, otpb and bhlhe22 expression has previously been detected in the zebra sh spinal cord, including in what appears to be the V0v domain (124). This suggests that these two genes may need Evx1/2 function to be turned off in V0v interneurons. In the zebra sh brain, Otpb is required and su cient for specifying aspects of a dopaminergic phenotype (138). Bhlhe22, which was previously known as Bhlhb5, has essential roles in mouse retinal development (139), axon elongation of corticospinal motor neurons in mouse (140), and survival of inhibitory neurons in the dorsal horn in mouse (141). In addition, it has been implicated in the formation of dI6, V1 and V2a spinal interneurons in chicken (142). However, analysis of a zebra sh bhlhe22 mutant detected no obvious differences in the spinal expression of en1b, evx1 and vsx2 in mutants compared to WT siblings, suggesting that V0v, V1 and V2a interneurons still formed in normal numbers (143). Therefore, the roles of these genes in V0v interneurons remain to be discovered.
We were intrigued to discover that hmx2 and hmx3a are expressed in mutant V0v interneurons ( Fig. 10B-C, Table 7, Fig. 11A-C, E-E"', Table 1). In a previous study, we showed that hmx2 and hmx3a are co-expressed in dI2 and V1 spinal interneurons, and that Hmx3a is required for the excitatory fates of dI2 interneurons. (Hmx2 also has a role in this process, but it is much more subtle than that of Hmx3a) (12). In the absence of Hmx3a function, many dI2 interneurons change their neurotransmitter fates from glutamatergic (excitatory) to GABAergic (inhibitory) (12), suggesting that Hmx3a is required to specify glutamatergic fates in these cells. Therefore, we were surprised that hmx2 and hmx3a are expressed in evx1;evx2 mutant V0v cells, given that these mutant cells have changed their neurotransmitter phenotype from excitatory to inhibitory. We con rmed this result using a combination of in situ hybridization and immunohistochemistry. We saw a large increase in the number of spinal cells expressing hmx3a in evx1;evx2 double mutants at 30 h ( Fig. 11A-C, Table 1). In addition, we used double-labelling experiments to show that V0v spinal interneurons express hmx3a in evx1;evx2 double mutants but not WT siblings (Fig. 11D-E'"). It is likely that these results re ect ectopic expression of hmx3a in mutant V0v interneurons, rather than hmx3a expression in cells that have transfated into dI2 or V1 interneurons (see discussion of mutant group 3 cells below), as we see the expanded hmx3a expression as early as 30 h (whereas we do not see expanded expression of other cell fate markers like en1b at these earlier stages (14)), and in our scRNA-seq data, we predominantly detect expression of hmx3a in mutant groups 1 and 2 rather than mutant group 3 (Fig. 10C). This temporal difference in hmx3a expression in evx1;evx2 double mutants (expanded at 30 h but not present in mutant group 3

cells at 48 h) is probably a consequence of the double mutant cells transfating into inhibitory spinal cells types by 48 h (see discussion below).
Taken together, these data suggest that Evx1/2 repress hmx3a expression in V0v interneurons. Therefore, we wondered if any of the genes that require Evx1/2 for their expression in V0v cells might be repressed by Hmx2/3a.
Interestingly, our data suggest that two of the six genes that we identi ed as requiring Evx1/2 function at 30 h, skor1a and nefma, are upregulated in dI2 interneurons in hmx2;hmx3a deletion mutants ( Fig. 4A-C, Fig. 5A-C, Fig. 12H-J, T-V, Z-AC'", Table 1), suggesting that Hmx2/3a may usually repress these genes in both dI2 and V0v interneurons ( Fig. 13D-E). However, in contrast, there is no statistically signi cant change in the number of spinal cord cells expressing skor1b, skor2, ebf3a or neff1 in hmx2;hmx3a deletion mutants (Fig. 12K-S, W-Y, Table 1). In combination, these data suggest that there are at least two distinct GRNs downstream of Evx1/2 in V0v neurons, one that includes repression of Hmx2/3a and one that is independent of Hmx2/3a (Fig. 13D).
Two molecularly distinct subsets of WT V0v interneurons exist at 48 h.
Our scRNA-seq analysis of V0v interneurons in embryos from an incross of evx1;evx2 heterozygous mutant parents, identi ed ve distinct clusters of cells. Based on our differential gene expression analyses, two are likely to be distinct WT clusters and the other three distinct mutant clusters (Fig. 6A). As discussed in the results, there are 15% more cells than we would expect in the WT clusters (1857 observed versus 1609 expected WT cells), compared to the mutant clusters (1003 observed versus 1251 expected mutant cells). This suggests that either the mutant cells were more fragile and, therefore, had a higher probability of not making it into our data set, or some of the mutant cells are contained in what we have de ned as the WT clusters. Both these explanations are possible. Due to their altered expression pro les, the mutant cells might be more likely to lose their integrity / become sick, in which case they would have been excluded from our analyses. It is also possible that some of the mutant cells have WT-like phenotypes and ended up in the WT clusters. Our previous analyses demonstrate that the V0v phenotypes of evx1 and evx2 single mutants are not completely penetrant. In single mutants, not all mutant cells lose expression of the genes that we analyzed, whereas more cells lose expression of these genes in double mutants (14). If this lack of penetrance persists to 48 h, we would expect some mutant cells to have a "WT" phenotype. Consistent with this, a small number of cells in the WT clusters in our UMAP plots appear to have partial mutant phenotypes where they express evx1 and/or evx2 but also express inhibitory genes, or they express markers of both glutamatergic and inhibitory fates (Supp. Figure 3).
Our discovery of two distinct subtypes of V0v interneurons at 48 h, is consistent with the existence of distinct molecular and/or functional subtypes of V0v interneurons in mouse (17,134) and adult zebra sh (144), and further highlights the conservation of neuronal speci cation between zebra sh and mammals. While we did not observe any obvious subtypes of V0v interneurons in our analyses of these cells at earlier developmental stages (14), the Higashijima group has previously identi ed three subsets of V0v interneurons with distinct morphologies that form at different times during the rst four days of development (145). V0-eA (commissural ascending) interneurons form rst (145), and these correspond morphologically to the neurons that we previously analyzed (14). At later stages V0-eB (commissural bifurcating), and then V0-eD (commissural descending) cells develop (145). These researchers reliably detected neurons with a V0-eB morphology at 60 h and neurons with a V0-eD morphology at 84 h (145). However, BrdU-labelling experiments showed that most V0-eB neurons are post-mitotic by 30-36 h. In contrast, most V0-eD neurons are not post-mitotic until 42-48 h (145). Therefore, it is possible that our WT groups 1 and 2 correspond to V0-eA and V0-eB cells. Consistent with this, the gene expression pro le of the WT group 1 cluster at 48 h more closely resembles what we saw at 27-30 h than WT group 2. There is statistically-signi cantly more expression of skor1b, skor2, ebf3a, nefma, nefmb, uncx and neff1 in WT group 1 than WT group 2 cells (Fig. 6I-L, N-P, Table 2A). As mentioned earlier, a small subset of V0v interneurons in mouse are cholinergic (134). However, it is unlikely that either of our WT groups correspond to these cells as both WT group 1 and WT group 2 contain too many cells and, in addition, we do not detect any expression of chatb in any of the cells in our 48 h data set and we only detect chata expression in a small subset of WT group 1 cells (6.75% (63/933 cells)) and a very small number of WT group 2 and mutant group 3 cells (0.97% (9/924 cells) and 5.97% (12/201 cells) respectively, data not shown). evx1;evx2 double mutant cells may transfate into distinct inhibitory interneurons, or motoneurons.
The most surprising result from our scRNA-seq analyses is the phenotype of the cells in the mutant group 3 cluster.
The lack of evx1, evx2 and slc17a6a expression in these cells and the increase in expression of markers of inhibitory cells, including slc6a5, slc6a1b and gad1b suggest that this cluster contains the most severe mutant cells, which presumably are the double mutant cells (Fig. 6A-G, Table 2). As discussed in the results section, the number of cells in this cluster is consistent with this hypothesis. However, the phenotype of these cells at 48 h is surprising, compared to what we have seen by in situ hybridization at earlier stages of development ( Fig. 5G-I, Fig. 6L, P, Table 1, Table 2). In addition, distinct subsets of mutant group 3 cells express markers of either different types of inhibitory spinal cord interneurons, or motoneurons (Fig. 8, Table 5). Taken together, these data suggest that the phenotype of mutant group 3 cells is distinct from the other two mutant groups.
One possible explanation of these different phenotypes is that the mutant group 3 cells are changing their gene expression pro les not because of their genotype, but as a side-effect of the experimental procedures. That some aspect of the experiment, for example the cell dissociation, has caused a subset of cells to aberrantly turn on genes that they wouldn't normally express. However, if this was the case, we would expect to see a wider variety of genes mis-expressed. We also might expect the cells to be sick, whereas all these cells appear healthy as they express very low levels of mitochondrial transcripts, typical of healthy cells (< 6%, data not shown), and pass all other stringent scRNA-seq quality controls (please see Methods for further information). In addition, the cell numbers are not consistent with this hypothesis. If this was the case, and cells of all genotypes were equally likely to be affected this way, then we would expect the ratio of WT cells (group 1 and group 2) to mutant cells in mutant group clusters 1 and 2 to be 9:7 (this assumes that the two mutant groups include both single and double mutant cells). Therefore, we would predict there to be 1496 WT cells and 1163 mutant cells. However, what we observe is 1857 cells in WT groups 1 and 2 and 802 cells in mutant groups 1 and 2 (P value for Chi squared test < 0.0001).
We consider that it is more likely that the mutant group 3 cluster consists of double mutant cells, and that by 48 h, these cells have started to transfate into different types of spinal neurons. Given that none of the potential ectopic fates that we observe are glutamatergic, and that they are instead inhibitory or cholinergic, this would also suggest that there is a feedback mechanism between neurotransmitter phenotype and cell type identity / cell fate. This hypothesis could also explain the large number of cells in mutant group 3 that express neff1, even though expression of this gene is downregulated in mutant groups 1 and 2, as neff1 is broadly expressed in the spinal cord at 48 h, and it is possible that it is expressed by motoneurons and/or inhibitory interneurons at this stage (Fig. 2AX, Fig. 6P, Table 2A).
Several studies have suggested that spinal cord fates are speci ed via a repression of repression mechanism, where fate-speci c genes inhibit all other possible fates rather than directly specifying the fate in question (e.g. 10, 73, 146, 147-150). Our data is similarly consistent with a mechanism where loss of Evx1/2 allows other non-V0v-cell fate-specifying genes to be expressed (Fig. 13F). Our analyses suggest that these different fates are non-overlapping, which would be consistent with a model where, as cells start to express one of these other fate-specifying genes, that gene then represses expression of other fate-specifying genes. Such a mechanism could stochastically produce different subsets of cells with distinct fates. It is currently not clear why these different ectopic fates are restricted to motoneurons and inhibitory interneurons. However, the rst phenotype that we observe in evx1;evx2 mutants is the loss of glutamatergic markers and gain of inhibitory markers. This suggests the intriguing hypothesis that there is something about this early change in neurotransmitter phenotypes that in uences the later change in other aspects of cell fate (Fig. 13F).
Interestingly, in mouse Evx1 mutants, expression of Evx2 is also lost in the spinal cord and about two thirds of V0v interneurons transfate into V1 interneurons, based on their expression of En1 (which, in the spinal cord, is only expressed by V1 cells) and their changed migration pattern and axon trajectories (17,70). While this previous mouse study did not examine neurotransmitter phenotypes, V1 interneurons are inhibitory (70,151,152). In our earlier analyses of zebra sh evx1;evx2 double mutants, we saw no evidence of V0v interneurons ectopically expressing en1b or Pax2 at 24 h or 30 h (Pax2 is expressed by V1 and several other inhibitory spinal interneurons) or changing their axon trajectories at 27 h or 48 h (14). However, our data in this paper suggest that evx1;evx2 double mutant V0v interneurons express markers of different inhibitory spinal neuron or motoneuron fates by 48 h. 28.36% (57/201 cells) of mutant group 3 cells express en1b, suggesting that almost a third of the double mutant cells may be transfating to V1 cells, compared to two thirds in the mouse study (Table 3). We do not think that we missed an increase in en1b expression in our earlier analyses as we saw a reduction, albeit not statistically signi cant, in the number of spinal cells expressing en1b at 24 h (14). Previously, we suggested that the mouse and zebra sh phenotypes might be different because dbx expression persists in at least some V0v cells in zebra sh (14). This might also explain why zebra sh evx1;evx2 double mutant V0v cells still develop V0v morphologies and do not initially express genes associated with other spinal cord fates. V0v and V0 D interneurons both develop from dbx-expressing progenitor cells, and they have similar axon trajectories (17,25,26). Therefore, it is possible that Evx1/2 are initially just required to specify the glutamatergic phenotype of V0v interneurons (V0 D neurons do not express Evx1/2 and are inhibitory), as Dbx can initially repress non-V0 cell fates. However, as Dbx function wanes, Evx1/2 may also be required to maintain V0v identities by repressing other interneuron and motoneuron identities (Fig. 13F). It remains unclear why more V0v cells acquire V1 fates in mouse than in zebra sh, although it is possible that this is also due to temporal differences in when Evx1/2 is required to maintain V0v identities.

Conclusions
In conclusion, this paper identi es two molecularly distinct subtypes of WT V0v spinal interneurons at 48 h. We also identify 25 transcriptional regulators that are downstream of Evx1/2 in V0v spinal interneurons at 30 h and/or 48 h that are, therefore, strong candidates for being members of the GRNs that specify the functional characteristics of these cells, plus 11 transcriptional regulators that are upregulated in V0v spinal interneurons at 48 h when Evx1/2 activity is reduced (nefma is in both of these groups as it is downregulated in evx1;evx2 mutants at 30 h and upregulated at 48 h). Interestingly, two of the transcriptional regulators that are upregulated in evx1/2 mutants are hmx2 and hmx3a, and we show that Hmx2/3a, in turn, repress expression of skor1a and nefma in dI2 interneurons. This suggests that Evx1/2 might regulate skor1a and nefma expression in V0v interneurons by repressing Hmx2/3a expression. Finally, our data suggest that in the absence of both Evx1 and Evx2, V0v spinal interneurons initially change their neurotransmitter phenotypes from excitatory to inhibitory and then at a later point of development, transfate into distinct types of inhibitory spinal interneurons, or motoneurons. Taken together, these ndings signi cantly increase our knowledge of V0v spinal development and move us towards a greater understanding of the GRNs that regulate this essential process.

Consent for publication
Not applicable Availability of data and materials Microarray data were previously deposited in the NCBI Gene Expression Omnibus with accession number GSE145916. Single-cell and bulk RNA-Seq data have been deposited in the NCBI Gene Expression Omnibus with accession numbers GSE240239 and GSE240238 respectively.

Competing interests
The authors declare that they have no competing interests.  Heatmap analysis of gene-expression pro ling of 27 h V0v spinal cord interneurons. A three-class ANOVA analysis of differential expression was performed on different FAC-sorted populations of cells. Class 1: All trunk cells. Class 2: All post-mitotic spinal neurons. Class 3: V0v interneurons. Each column is a different biological replicate. Rows show relative expression levels for a single gene as normalized data transformed to a mean of 0, with standard deviation of +1 (highly expressed, red) or -1 (weakly/not expressed, blue) sigma units. Adjusted P-values corrected for multiple testing are shown on the left-hand side. Expression pro les for positive control genes evx1 and evx2, whose spinal cord expression is exclusive to V0v interneurons, are shown. The high level of expression of these genes in our V0v samples, compared to the other samples, con rms that we have successfully isolated V0v interneurons. Additional positive control genes slc17a6a and slc17a6b, con rm that V0v interneurons are excitatory (glutamatergic), whereas negative control genes slc6a9, slc6a5, gad1b and gad2 show that V0v interneurons do not express either glycinergic or GABAergic inhibitory neurotransmitter pathway genes and that there is no contamination of our V0v samples with inhibitory neurons. The expression pro les for slc17a6a, slc17a6b, slc6a9, slc6a5, gad1b and gad2are reproduced from (14) as per the Creative Commons Attribution (CC BY) license at Neural Development.   (14). We are showing additional skor2 data here to demonstrate reproducibility of our co-expression experiments, and for ease of comparison with the skor1b, uncx, and nefma data. White asterisks indicate double-labelled V0v interneurons. Cells that are green and not red could be V0v interneurons that do not express the gene in question, or V0v interneurons with low expression, not revealed in these experiments, of the gene detected in red. We often detect fewer cells expressing a particular gene in double-labelling experiments where the mRNA is detected with a red uorophore, than in single in situ hybridization experiments where the mRNA is detected with NBT/BCIP (viewed as an opaque blue stain under visible light), suggesting that the weakest-expressing cells may not be detected in the former, probably due to the prolonged processing of samples necessitated by uorescent double-labelling experiments, which can affect the stability of target mRNA molecules, and the lower sensitivity of the red label. Therefore, we cannot conclude for certain that single-labelled EGFP-positive cells, do not express the gene detected in red. Scale bar: (  (both values are also listed in Table 1). All counts are an average of at least three embryos. Statistically signi cant comparisons are indicated with brackets and asterisks. *** P <0.001. * P <0.05. White circles indicate wild-type data and black circles indicate evx1;evx2 double mutant data. All data were analyzed for normality using the Shapiro-Wilk test. Data in L is not normally distributed and so a Wilcoxon-Mann-Whitney test was performed. Data sets in C, F, I and O are normally distributed and so the F-test for equal variances was performed, followed by a type 2 Student's t-test (for equal variances). P-values are provided in Table 1 Table 1). All counts are an average of at least four embryos. Statistically signi cant comparisons are indicated with brackets and asterisks. *** P <0.001. * P <0.05. White circles indicate wild-type data and black circles indicate evx1;evx2 double mutant data. All data were analyzed for normality using the Shapiro-Wilk test. Data in C is not normally distributed and so a Wilcoxon-Mann-Whitney test was performed. Data sets in F, I and L are normally distributed and so the F-test for equal variances was performed, followed by a type 2 Student's t-test (for equal variances). P-values are provided in Table 1  Tg(evx1:EGFP) SU2 . Clusters are color-coded by cell identity: V0v WT group 1 (light green), V0v WT group 2 (dark green), V0v mutant group 1 (turquoise), V0v mutant group 2 (light blue), and V0v mutant group 3 (dark blue). Cell fate assignments were deduced and extrapolated by comparing expression pro les of 48 h single-cell clusters with the molecular phenotypes of V0v spinal interneurons in wild-type and evx1 and evx2 single and double mutant embryos (14). (B-Q) UMAP plots of differential gene expression between cell clusters. Black shows high levels of expression, light grey shows low levels of expression. All expression data have been normalized (see Methods). (B-D) Many of the cells in both WT clusters express (B) evx1 and/or (C) evx2, as well as the glutamatergic marker (D) slc17a6a. (B-G) evx1, evx2 and slc17a6a are all detected in fewer cells in mutant groups 1 and 2 and hardly any cells in mutant group 3. Many cells in the mutant clusters upregulate inhibitory markers, including (E) slc6a5, (F) slc6a1b, and (G) gad1b. (H) skor1a and (I) skor1b are not detected in many cells in this data set. (H) skor1a is expressed in a few WT group 1 and 2 cells, as well as a couple of mutant group 3 cells and a mutant group 1 cell. (I) skor1b is predominantly detected in a few WT group 1 cells. (J) In contrast, skor2 is expressed at high levels in most V0v WT group 1 cells, and it is also detected in multiple WT group 2 cells and a small number of mutant group 2 cells. (K) ebf3a has a similar expression pro le to skor2, except that its expression is also detected in a few mutant group 1 cells and slightly more mutant  from 48 h embryos produced from an incross of evx1 i232/+ ;evx2 sa140/+ heterozygous parents homozygous for Tg(evx1:EGFP) SU2 . Clusters are color-coded by cell identity: V0v WT group 1 (light green), V0v WT group 2 (dark green), V0v mutant group 1 (turquoise), V0v mutant group 2 (light blue), and V0v mutant group 3 (dark blue). For ease of cell type comparison, panel 7A has been reproduced from Fig 6A. (B-J) UMAP plots of differential gene expression between cell clusters. Black shows high levels of expression, light grey shows low levels of expression. All expression data have been normalized (see Methods). (B) anos1a, (C) chrna2b, (D) fndc4b, (E) plpp4, (F) cnih3, and (G) drd2b are all expressed in more cells in WT and mutant group 1 clusters than WT and mutant group 2 clusters. In contrast, (H) esrrb, (I) scxa, and (J) svild are all expressed in more cells in WT and mutant group 2 clusters than WT and mutant group 1 clusters.      Table 1). All counts are an average of ve embryos. Statistically signi cant comparison is indicated with brackets and asterisks. * P <0.05. White circles indicate wild-type data and black circles indicate evx1;evx2 double mutant data. All data were rst analyzed for normality using the Shapiro-Wilk test. Both data sets in C are normally distributed and so the F-test for equal variances was performed, followed by a type 2 Student's t-test  (Fig. 2). Model I shows