Epha1 is a cell-surface marker for the neuromesodermal competent population

Vertebrate body axis extension occurs in a head-to-tail sequence and relies on populations of cells with self-renewing properties, collectively known as axial progenitors (Aires et al., 2018; Dias and Aires, 2020; Wymeersch et al., 2021). These progenitors, initially located in the caudal epiblast and later in the tailbud, include several cell pools classified according to their potential. The neuromesodermal competent population (NMC population), previously referred to generally as neuromesodermal progenitors (NMPs) (see Binagui-Casas et al., 2021) is known to harbor cells with developmental potential to generate both post-occipital neural and mesodermal derivatives (i.e. NMC cells) (Binagui-Casas et al., 2021; Cambray and Wilson, 2002, 2007; Guillot et al., 2021; Henrique et al., 2015; Tzouanacou et al., 2009; Wymeersch et al., 2021). Recently, transcriptome data analyses helped to identify another population, the notochord progenitors, that seems to maintain a stable molecular signature during axial elongation (including genes such as Noto, Shh and Foxa2), possibly acting as a stable progenitor niche, particularly for the NMC population (Wymeersch et al., 2019). Since the initial identification of axial progenitors, increasing efforts have been conducted to determine the precise molecular characteristics of each different cell population, with particular focus on the NMC population, and the capacity of individual cells to actually generate both post-occipital neural and mesodermal derivatives [i.e. the ‘true’ neuro-mesodermal progenitors (NMPs; see Binagui-Casas et al., 2021)]. Combined mapping and expression studies indicate that the early neural marker Sox2 and the mesodermal transcription factor brachyury (also known as Tbxt or T) are mainly co-expressed in regions known to contain NMC cells and, more recently, also in NMPs (Binagui-Casas et al., 2021; Cambray and Wilson, 2007; Guillot et al., 2021; Martin and Kimelman, 2012; Olivera-Martinez et al., 2012; Tsakiridis et al., 2014; Wymeersch et al., 2016, 2021). The introduction of improved high-throughput techniques allowing in vivo and in vitro transcriptomic analyses, and in particular single-cell RNA sequencing (scRNA-seq), has shed additional light on the role of these key transcription factors and has provided a deeper understanding of other molecular players (e.g. Wnt3a, Gdf11, Lin28 and Cdx2) and gene regulatory networks involved in the maintenance and differentiation of the NMC population, both in the embryo and in vitro (Aires et al., 2019; Dias et al., 2020; Edri et al., 2019; Gouti et al., 2017; Guibentif et al., 2021; Guillot et al., 2021; Koch et al., 2017; Wymeersch et al., 2019, 2021). Interestingly, some of these studies have shown that the transcriptome of the NMC population changes extensively over time, including the activation of an incomplete epithelial-to-mesenchymal transition (EMT) when entering the tailbud (Dias et al., 2020; Guillot et al., 2021; Wymeersch et al., 2019). Additionally, expression of Tbx6 in a subset of cells within the tailbud NMC population has also been reported, and this expression is involved in neuro-mesodermal fate decisions (Javali et al., 2017). Gene expression and lineage-tracing experiments have also indicated that Nkx1-2 (previously known as Sax1) is present in the NMC region, and in early neural and mesodermal progenitors throughout axial extension (Albors et al., 2018).

Despite all these studies, it is still not possible to isolate the NMC population in a physiologically active form without relying on previous modifications to introduce reporter genes, mostly because the molecules typically used to identify those progenitors are transcription factors. In the present work, we searched for cell-surface markers that could facilitate isolation of the NMC population using conventional cell-sorting approaches. For this, we exploited the intrinsic self-renewal properties of these progenitors to obtain a cell population highly enriched in NMC cells (and possibly NMPs) from the tailbud. In our approach, we introduced a fluorescent marker into the axial progenitors early in development, which permitted us to identify and isolate the long-term progenitors (i.e. cells from the NMC population that are present for long periods of axial elongation, potentially NMC/NMP cells; Cambray and Wilson, 2007) from the tailbud at later developmental stages. From the genes enriched in these cells, we focused on Epha1 because it encodes a cell-surface protein that is mostly localized in the progenitor-containing areas of the mouse embryo throughout axial extension and because of the existence of fluorescence-activated cell sorting (FACS)-tested antibodies that could facilitate the development of protocols for the isolation of NMC cells from embryonic tissues. Using both whole-mount in situ hybridization and immunofluorescence, we detected stronger Epha1 expression in regions of the caudal epiblast and the tailbud known to contain NMCs/NMPs and their early mesoderm derivatives. Importantly, we show that Epha1-positive cells isolated from the mouse tailbud were able to produce neural and mesodermal derivatives when cultured in vitro under NMP differentiating conditions. In addition, Epha1-positive cells from the tailbud were highly enriched in Sox2⁺/Tbxt⁺ cells, the most frequently used molecular signature of the NMC population (Binagui-Casas et al., 2021). FACS profiles obtained from tailbud cells revealed that Epha1-positive cells can be divided in two subpopulations, Epha1^High and Epha1^Low, according to their Epha1 protein content. Whereas cells with low Epha1 expression levels were also obtained from tissues containing early differentiated NMC-population derivatives, Epha1^High cells were exclusively found in axial progenitor regions. Analysis of the transcriptomes obtained from tailbud Epha1^High and Epha1^Low cells revealed that, in addition to an enrichment in NMC population-related transcripts, the profile obtained from Epha1^High cells suggested that progenitors were entering mesodermal routes, whereas Epha1^Low cells included a molecular signature congruent with the existence of notochord progenitors within this cell compartment. Together, our results indicate that Epha1 is a valuable cell-surface marker for the isolation of the NMC population from the mouse tailbud and that high Epha1 values might be a molecular signature of progenitors entering the mesodermal progenitor compartment.

Most transcriptome analyses of the NMC population from embryonic tissue have been performed on microdissected regions of the embryo containing NMCs/NMPs that also include other cell types (e.g. some that had already entered mesoderm differentiation routes) or on cells isolated on the basis of the expression of a gene enriched in these progenitors (Dias et al., 2020; Gouti et al., 2017; Guillot et al., 2021; Koch et al., 2017; Wymeersch et al., 2019). Here, we used an unbiased approach to isolate NMPs exploiting their self-renewing properties. In particular, as the NMC population in the tailbud is derived mainly from the cells located in the node streak border (NSB) at earlier developmental stages (Cambray and Wilson, 2007), we designed a genetic strategy to label axial progenitors when they are part of the caudal epiblast and isolated them from the tailbud at later developmental stages. This system combines the Cdx2P-Cre^ERT transgene (Jurberg et al., 2013) with the ROSA26-YFP-Cre reporter (Srinivas et al., 2001). By administering a single low tamoxifen dose at embryonic day (E) 7.5, we induced a short pulse of permanent YFP label into a subset of axial progenitors that allowed their fate to be followed at later developmental stages. We have previously used this system with the ROSA26-βgal-Cre reporter (Soriano, 1999), which proved the efficiency of the method, as estimated by the presence of labeled cells in the embryonic tissues posterior to the position of Cre-mediated recombination, all the way down to the tail tip (Aires et al., 2019) (Fig. 1A-C).

Descendants of progenitors labeled with YFP at E7.5 were then recovered from the tail region at E10.5 by FACS. Two sets of YFP-positive cells were recovered: the first set was isolated from the tailbud (Tail^Prog), which is expected to be highly enriched in progenitors; a second set was recovered from more-anterior tail regions, where labeled cells are already part of the tissues derived from the progenitors (Tail^Desc). Comparison of the transcriptome of these two cell pools, obtained by RNA-sequencing (RNA-seq), identified 2465 genes showing differential expression (P<0.05) between the two cell groups (Table S1). Of these, 847 genes were highly expressed in the Tail^Prog, whereas 611 genes were upregulated in the Tail^Desc sample. A selection of 12 differentially expressed genes, including NMC population-related genes and others encoding membrane proteins (Epha1, Cldn9 and Nkd2), was then used to validate the RNA-seq data by reverse transcription-quantitative polymerase chain reaction (RT-qPCR) (Fig. 1E).

Initial analysis of the differentially expressed genes revealed high enrichment of the Tail^Prog cells in factors that have been linked to neuro-mesodermal identity (Fig. 1D). For example, Cdx2, Cdx4 and Tbxt, which are known to be highly expressed in the NMC population and proven to be essential for their activity (Amin et al., 2016; Chawengsaksophak et al., 2004; Herrmann et al., 1990; Savory et al., 2011; van Rooijen et al., 2012), were among the most strongly upregulated genes in the Tail^Prog compartment (Fig. 1D). Similarly, other genes whose expression has been shown to be enriched in tailbud NMC region, including Mnx1, Nkx1-2, Fgf8, Fgf4, Fgf3, Evx1, Cyp26a1, Hoxb1, Sp5, Gdf11, Wnt5a or Wnt3a (Abu-Abed et al., 2001; Aires et al., 2019; Albors et al., 2018; Boulet and Capecchi, 2012; Cambray and Wilson, 2007; Dush and Martin, 1992; Greco et al., 1996; Guillot et al., 2021; Harrison et al., 2000; Javali et al., 2017; McPherron et al., 1999; Murphy and Hill, 1991; Naiche et al., 2011; Robinton et al., 2019; Sakai et al., 2001; Takada et al., 1994; Wymeersch et al., 2019; Yamaguchi et al., 1999), also showed significant differential expression in the Tail^Prog cell pool (Fig. 1D). In addition, we also found other genes previously not linked to the NMC population, such as Epha1, Arl4d, Efna1, Cldn9, Gad1 and Scara5, to be highly expressed in the progenitor compartment (Fig. 1D).

Conversely, we found enrichment of the Tail^Desc cell pool in markers for neural and mesodermal derivatives of the NMC population, including Ngn2, Sox1, Olig2, Olig3, Sox10, Meox2, Raldh2 (Aldh1a2), Pax6, Fst and Nr2f2 (Albano et al., 1994; Aubert et al., 2003; Candia et al., 1992; Gradwohl et al., 1996; Jonk et al., 1994; Kuhlbrodt et al., 1998; Niederreither et al., 1997; Takeichi et al., 2002; Walther and Gruss, 1991) (Fig. 1D), indicating that the Tail^Prog cells have the differentiation potential expected from the NMC population. Sox2, one of the components of the typical NMC-population signature, was expressed at slightly higher levels in Tail^Desc than in Tail^Prog cells (Fig. 1D). This is not surprising, as Sox2 is highly expressed in the neural tube, which is one of the tissues from which the Tail^Desc cells were recovered. Interestingly, both Tail^Prog and Tail^Desc datasets do not contain notochord markers such as Shh, Noto or Foxa2, which were readily found in a dataset obtained from similarly staged unsorted tailbuds (Aires et al., 2019) (Tail^Tot) (Fig. 1F; Table S2), indicating that the Tail^Prog cell pool is specifically enriched in cells from the NMC population. Together, these data indicate that the genetic lineage-tracing strategy described here labels the tailbud NMC populations and their neural and mesodermal descendants, and that the isolated Tail^Prog cell pool is highly enriched in cells belonging to the NMC population.

From the genes differentially expressed (P<0.05) between the Tail^Prog and the Tail^Desc cell populations, we concentrated on those encoding membrane proteins that could be used to isolate physiologically active cells of the NMC population without previous genomic modifications (e.g. transgenic reporters). Gene ontology (GO) categorization (Ashburner et al., 2000) of genes differentially upregulated in Tail^Prog cells with a log₂ fold change>2 and q<0.05, identified 61 genes assigned to the category ‘membrane’ (GO:0005886). From these, we further selected the genes with higher expression, reducing the list to 16 genes (Fig. S1A). Expression analyses at E10.5 by in situ hybridization revealed that the staining patterns for some of those genes included a strong signal in the tail region, although these patterns differed among the various genes (Fig. S1B). In addition, for most of them, tailbud expression represented a subset of more-complex expression patterns that included other embryonic regions.

From these genes, we focused on Epha1 based not only on its expression pattern, which is apparently restricted to the progenitor zone at different developmental stages (Fig. S1B) (Duffy et al., 2006), but also on the availability of FACS-validated antibodies, which are able to provide reliable data with cells obtained from solid embryonic tissues. By in situ hybridization, we detected Epha1 expression as early as E8.5, when it was mostly observed in the caudal lateral epiblast (CLE), a region containing a NMC population (Binagui-Casas et al., 2021; Cambray and Wilson, 2007; Wymeersch et al., 2016), as well as in cells entering in the mesodermal compartment (Fig. 2Aa,a′). Epha1 expression was also observed in the caudal epiblast of E9.5 embryos, fading anteriorly when entering the regions corresponding to the presomitic mesoderm and the caudal neural tube (Fig. 2Ab,b′); at E10.5 and E11.5, Epha1 expression was essentially restricted to the tailbud (Fig. 2Ac-d′). Analysis of Epha1 by immunostaining confirmed the expression patterns obtained by in situ hybridization (Fig. 2B), although slightly broader, possibly derived from higher stability of the protein than of the mRNA. In addition, although this technique is not purely quantitative, it showed that the Epha1 protein levels were not uniform in the positive domain. In particular, in E8.5 embryos, Epha1 expression was strong in the region of the primitive streak (PS) (Fig. 2Ba3, white arrow), which is known to contain progenitor cells undergoing an EMT to generate the mesodermal layer (Acloque et al., 2009; Hay, 1968; Wilson et al., 2009), and at E10.5 in the region abutting the posterior end of the neural tube, which is positive for both Tbxt and Sox2, as well as in mesenchyme caudal to this region thought to contain mesoderm progenitors (Dias et al., 2020; McGrew et al., 2008; Wymeersch et al., 2016, 2019) (Fig. 2Bc1, white arrows). Analysis of published mouse bulk RNA-seq datasets (Koch et al., 2017; Wymeersch et al., 2019) indicates that Epha1 is not only strongly expressed in NMC regions [NSB, CLE and chordoneural hinge (CNH)] and their early mesoderm derivatives (e.g. posterior CNH), but is also highly enriched specifically in isolated cells with the Tbxt⁺/Sox2⁺ phenotype, the currently most-used criterion to identify the NMC population (Fig. S2A,B). At the single-cell level, using an early organogenesis stage (Theiler stage 12) mouse scRNA-seq dataset (Pijuan-Sala et al., 2019), we also found Epha1 to be highly expressed in the clusters identified as the NMC population (previously referred to as NMPs) and the caudal mesoderm (Fig. 3), overlapping with the region co-expressing genes associated with the NMC population. Consistent with this observation, we also found that Epha1 is present in several CLE single cells (Gouti et al., 2017) that co-express Tbxt and Sox2 at E8.5 (Fig. S2C). Our expression data, together with the published transcriptomes, therefore suggest an association between Epha1 expression and the NMC population.

To further characterize the Epha1-positive cell population, we analyzed cells from the tailbud and adjacent anterior tail region of E10.5 embryos by FACS, using an antibody against Epha1. Both areas contained a high proportion of Epha1-positive cells (Fig. 4A-C), which fits with the immunofluorescence data. Interestingly, we noticed that the staining patterns in the two embryonic regions were different, as cells from the tailbud included a population with higher staining intensity that was never observed in the FACS plots obtained from the anterior tail region (Fig. 4B). These results indicate the existence of cells with different amounts of Epha1 in the tailbud, consistent with the non-uniform staining intensities observed by immunofluorescence. Based on their Epha1 content, we grouped Epha1-positive tailbud cells into two subpopulations, which will be referred to as Epha1^High and Epha1^Low (Fig. S3). We also obtained similar Epha1^high and Epha1^low cell compartments in the progenitor-containing area of E8.5 embryos, although the proportion of Epha1^high cells was lower in this tissue than in the tailbud from E10.5 embryos (Fig. S4).

We then assessed the distribution of Sox2⁺/Tbxt⁺ cells among the different Epha1 compartments obtained from the tailbud of E10.5 embryos. We observed that the proportion of cells co-expressing Sox2 and Tbxt differed significantly among the three cell compartments, being (on average) 57% for Epha1^High cells, 29% for Epha1^Low cells and almost absent in the Epha1-negative pool (3%), consistent with cells from the NMC population being present in the Epha1-positive compartment (Fig. 4D-F; Table 1). We then performed a complementary analysis by first isolating Sox2⁺/Tbxt⁺ cells from the E10.5 tailbud and analyzing their distribution among the different Epha1 compartments. We observed that, on average, 79% of the cells co-expressing Sox2 and Tbxt were Epha1^High, whereas only 9.8% and 4% of these cells were Epha1^Low or Epha1 negative, respectively (Fig. 4G-I; Table 2). These results are consistent with the Epha1^High subpopulation being enriched in the NMC cell population.

A similar distribution of Sox2⁺/Tbxt⁺ cells among the various Epha1 compartments and of Epha1^High and Epha1^Low from the Sox2⁺/Tbxt⁺ gated cells was also observed at E8.5 (Fig. S4). On average 53.4% of Epha1^High and 22.5% of Epha1^Low were found within the Sox2⁺/Tbxt⁺ cells (Fig. S4A-C; Table 1), and 70% and 18.4% of the Sox2⁺/Tbxt⁺ cells showed Epha1^High and Epha1^Low profiles, respectively (Fig. S4D-F; Table 2). In addition, NMC populations obtained from mouse ES cell cultures according to standard inducing conditions (Gouti et al., 2014; Turner et al., 2014) were preferentially found in the Epha1^High compartment, from which an average of 91% were positive for both Tbxt and Sox2 (Fig. 4J-L; Table 1), also further suggesting enrichment in cells from the NMC population in this Epha1-positive compartment of in vitro differentiated ES cells. Consistent with this, about 95% of Sox2⁺/Tbxt⁺ cells obtained from ES cell-derived NMC populations mapped to the Epha1^High compartment (Fig. 4M-O; Table 2). The proportion of Epha1-negative cells in the in vitro-derived NMC populations was almost residual, matching the high proportion of Sox2⁺/Tbxt⁺ cells in this compartment (Fig. 4N-O; Tables 1 and 2). Together, taking the simultaneous expression of Tbxt and Sox2 as the defining criterion for the NMC population, these results suggest that these cells might be mostly located within the Epha1^High compartment, although a significant number can also have an Epha1^Low phenotype, at least in embryonic tissues. This might reflect heterogeneity within the NMC population, when defined only as Sox2⁺/Tbxt⁺ cells, and is in line with observations made by other laboratories (e.g. Romanos et al., 2021; Wymeersch et al., 2016).

To test whether the Epha1-positive compartment have the functional properties expected for the NMC population, we FACS to separate Epha1-positive and -negative cells from E10.5 embryonic tailbuds, and evaluated their capacity to differentiate in vitro along the neural and mesodermal routes (Fig. 5, Figs S5 and S6). When Epha1-positive cells were incubated under neural-promoting conditions, we observed sustained expression of Sox2 together with the absence of mesodermal markers (Fig. 5B and Fig. S6B,E). This result matches the molecular phenotype observed during neural differentiation of in vitro-derived NMC populations (Gouti et al., 2014, 2017). Interestingly, after 3.5 days under neural differentiation conditions, we observed the presence of Sox2-positive structures resembling neural rosettes, as well as the presence of neurite-like cellular extensions, suggesting the presence of neurons in these cultures. This was confirmed by the identification of Tuj1 expression in the periphery of the rosettes associated with the neurites (Fig. 5C). Conversely, when Epha1-positive cells were incubated under mesoderm-promoting conditions, they activated Tbx6 expression, indicating that they were able to enter paraxial mesodermal fates (Fig. 5A and Fig. S6A,E), which is also consistent with the patterns observed during mesoderm differentiation of in vitro-derived NMC populations (Gouti et al., 2014, 2017). It should be noted that clumps of Sox2-positive cells were often observed mixed with Tbx6-positive cells (Fig. 5A), although they both segregated from each other, even when present in the same cluster. A similar observation has been previously reported with in vitro-differentiated NMC populations (Gouti et al., 2017). So far, we have no explanation for the presence of those Sox2 clusters under mesodermal differentiation conditions. Importantly, Epha1-negative cells mostly failed to form colonies in vitro under the same culture conditions, remaining as single cells, and the few colonies formed by these cells were unable to follow neural or mesodermal fates when incubated under differentiation conditions (Fig. S6E). Epha1-positive sorted cells from in vitro-derived NMC populations, obtained from mouse embryonic stem (ES) cells, were also able to enter neural and paraxial mesodermal differentiation routes (Fig. 5D-F and Fig. S6C,D). In this case, however, it was not possible to compare with the Epha1-negative population from the same cultures, because they were present in only residual amounts (Fig. 4J and Fig. S5). Together, these results indicate a close functional association between Epha1 expression and the NMC population.

To further characterize the Epha1-positive cell populations in the tailbud, we isolated the Epha1^High and Epha1^Low compartments from this region of E10.5 embryos, and analyzed their transcript content by RNA-seq. Both cell compartments showed enrichment in genes that have been associated with the NMC population (e.g. Tbxt, Wnt3a, Tbx6 and Nkx1-2) to levels similar to those observed in Tail^Prog cells, an effect particularly clear in Epha1^High cells (Fig. 6A). However, Epha1^High and Epha1^Low cells seemed to represent two different cell compartments, as 875 genes showed differential expression (P<0.05) between these two Epha1 cell populations (Table S3). These included several early mesoderm-related genes, such as Tbx6, Dll1, Msgn1, Cited1, Lfng, Snai1 or Wnt3a (Carver et al., 2001; Dias et al., 2020; Gouti et al., 2017; Koch et al., 2017; Serth et al., 2003; Wymeersch et al., 2019), that were upregulated in the Epha1^High compartment (Fig. 6B), suggesting that this cell population might contain cells from the NMC population that are already entering mesodermal fates. Interestingly, we also found that this cell population displays a more mesenchymal state than that observed in the axial progenitor pool (Tail^Prog) (Fig. 6C), thus suggesting that cells in the tailbud NMC population might complete the partial EMT that brought them into the tailbud (Dias et al., 2020) when entering mesodermal fates.

The transcriptome profile of Epha1^Low cells also showed substantial enrichment in some markers for the NMC population, which is consistent with the presence of Sox2⁺/Tbxt⁺ cells in this cell compartment (Fig. 6A). Interestingly, these cells also contained a high abundance of transcripts for Shh, Noto, Foxa2, Lmxa1, Krt8 and Krt18, which are not present in either the Epha1^High or the Tail^Prog cell populations (Figs 6B and 1D). These genes are commonly expressed in the notochord (Abdelkhalek et al., 2004; Ang et al., 1993; Echelard et al., 1993; Rodrigues-Pinto et al., 2016; Wymeersch et al., 2019), suggesting that the tailbud Epha1^Low cell population might also be enriched in notochord progenitors. In addition, this compartment also expressed markers for other midline structures, including those in common with the notochord (Foxa2 and Shh), and others specific for the floor plate, such as Slit2 and Spon1 (Brose et al., 1999; Klar et al., 1992) (Fig. 6B), suggesting that the Epha1^Low cells might also contain progenitors for the ventral spinal cord.

Together, our results indicate that Epha1 is a valuable cell-surface marker that can be used to label and isolate cells from the NMC population, and is able to generate neural and mesodermal derivatives. In addition, although higher Epha1 expression might be a hallmark of cells of the NMC population entering mesodermal routes, the midline progenitor population contributing to notochord and floor plate might be included within the low Epha1 expression compartment.

In this study, we have used a genetic strategy to label and isolate physiologically active cells from the NMC population of the tailbud of mouse embryos and analyzed their mRNA content using high-throughput methods. Although the Cre driver used to label the cells is not expected to be exclusive for the progenitors at early stages, it is expected that from the cells labeled by a pulse of Cre activity, only bona fide progenitors will carry the reporter to the tailbud. Indeed, labeled cells can be observed along the whole axis caudal to the region where the reporter was activated (see Fig. 1B or Aires et al., 2019). Both the enrichment of tailbud cells isolated with this strategy (Tail^Prog) in known markers for the NMC population, as well as the high concordance between our dataset and previously published data for that progenitor population obtained from embryos or in vitro-differentiated human or mouse ES cells (Dias et al., 2020; Gouti et al., 2017; Koch et al., 2017; Verrier et al., 2018; Wymeersch et al., 2019), support our conclusion that our lineage-tracing strategy is a non-biased reliable tailbud axial progenitor labeling method. However, it should be noted that the labeling conditions used to restrict the time frame of effective cell labeling might also have reduced the number of effectively labeled progenitors. Therefore, the tailbud cells isolated and analyzed using this approach might represent only a fraction of the actual progenitors and it is therefore possible that the tailbud NMC-population molecular fingerprint obtained in these experiments is not fully comprehensive.

Our initial goal was to identify cell-surface markers that could facilitate isolation of physiologically active cells from the NMC population. From the list of new genes that we found to be differentially expressed between NMC populations and their derivatives, only Epha1 showed expression restricted to the progenitor-containing areas. Whether any of the other candidate surface molecules, alone or in combination, could help in developing protocols for NMP isolation will require additional work.

Functional and expression criteria indicate that the NMC population is contained within the Epha1-positive cell compartment of the tailbud. On the one hand, the Epha1-positive cells were enriched in Sox2⁺/Tbxt⁺ cells, the standard signature for the NMC population. In addition, and most importantly, Epha1-expressing cells were able to generate both mesodermal and neural derivatives when cultured in vitro under the regular neuro-mesodermal differentiation conditions, whereas Epha1-negative cells obtained from the same tissue were mostly unable to grow and failed to enter any of those differentiation fates under the same conditions. From our data, it is not possible to determine whether all or only a subset of Epha1-expressing cells has NMC properties. This actually resembles the previous observation that only a subset of the Sox2⁺/Tbxt⁺ cells, probably those expressing moderate levels of both proteins, might be bona fide NMC cells (Wymeersch et al., 2016).

Our experiments also indicate that the levels of Epha1 are not uniform in the caudal epiblast or the tailbud cells expressing this molecule. This was clearer from the FACS profiles showing the existence of at least two populations of Epha1-positive cells in these embryonic regions. The observation that the Epha1^High population was particularly enriched in Sox2⁺/Tbxt⁺ cells suggested that it could represent the bona fide marker for the NMC population. However, the association between Epha1 expression levels and the NMC-population identity might be more complex, as a significant number of Sox2⁺/Tbxt⁺ cells were also found in cells with lower Epha1 expression levels, indicating that the NMC population could be distributed between the Epha1^High and Epha1^Low compartments. Consistent with this, the transcriptome profile of tailbud Epha1^Low cells shows considerable enrichment in genes known to be expressed in the NMC population. The issue, then, is to understand the identity of cells containing high and low Epha1 expression levels. The transcriptome profile of the Epha1^High compartment, in addition to the presence of markers for the NMC population, also contains a significant enrichment in mesoderm-associated genes, thus indicating that at least part of these cells could represent progenitors entering mesodermal fates. Indeed, scRNA-seq analysis of Theiler stage 12 (around E8.5) (Pijuan-Sala et al., 2019) indicates that Epha1 is strongly expressed in cells that were classified as belonging to NMC population and caudal mesodermal clusters, colocalizing with the Sox2⁺/Tbxt⁺/Nkx1-2⁺ domain. The distribution of higher Epha1 immunoreactivity in the caudal tissue of E8.5 embryos and in the tailbud is also consistent with this hypothesis, as it is associated with areas representing early steps of mesodermal formation (e.g. primitive streak) (McGrew et al., 2008; Wymeersch et al., 2016, 2019). Similarly, the general patterns observed within the Sox2⁺/Tbxt⁺ gate in the FACS plots from Epha1^High cells might indicate the presence of cells at early stages in the progression towards the mesodermal lineage, as the cell distribution in the Sox2 expression axis includes a significant proportion of cells with relatively low Sox2 content but comparatively high Tbxt expression. This pattern is not observed in the Sox2⁺/Tbxt⁺ cells from the Epha1^Low gate, further highlighting the differences between the Epha1^High and Epha1^Low compartments. Intriguingly, the transcriptomic profile of Epha1^Low cells indicates that they might also contain precursors for the notochord and the floor plate, which are clearly not present in the Epha1^high compartment. This finding is interesting in light of observations in chicken and zebrafish embryos that indicate the existence of a common progenitor for both structures (Row et al., 2016; Teillet et al., 1998). The presence of a similar progenitor in mice is not clear (Jeong and Epstein, 2003). Although our RNA-seq data cannot provide definitive proof for the existence of a bipotent progenitor contributing to both the notochord and floor plate, they might provide a new tool to investigate this interesting issue. Based on the above observations, we suggest that Epha1 is associated with axial progenitors, including the NMC population and those for the midline (notochord and floor plate), and that a transient increase in Epha1 expression is a hallmark of cells within the NMC population entering the mesodermal progenitor compartment that will generate paraxial mesoderm during axial extension. Interestingly, tailbud Epha1^High cells seem to contain a more-complete mesenchymal profile than the tailbud NMC population, which has been shown to contain an intermediate epithelial/mesenchymal signature derived from an incomplete EMT (Dias et al., 2020; Guibentif et al., 2021; Guillot et al., 2021). This finding might indicate that, in the tailbud, mesodermal differentiation from the NMC population includes completion of the partial EMT that is characteristic of the axial progenitors. Such an intermediate EMT state has also been recently described for axial progenitors in zebrafish embryos, acting as a checkpoint to guarantee proper control of their differentiation (Goto et al., 2017; Kinney et al., 2020), indicating that epithelial-to-mesenchymal plasticity might be a general feature of these cells in different vertebrate clades (Binagui-Casas et al., 2021).

In conclusion, Epha1 and its expression levels might be a useful marker with which to access cells within the NMC-region/population in a physiologically active state without resorting to previous incorporation of genetic reporters. However, whether this molecule plays a specific role in NMCs and/or NMPs behavior remains to be determined. Epha1 is a member of the ephrin receptor family (Lisabeth et al., 2013). These receptors interact with ephrins on the surface of adjacent cells to control a variety of differentiation and morphogenetic processes involving tissue compartmentalization. It is therefore possible that increased Epha1 levels could be part of the mechanism promoting progenitor exit from their niche and their segregation into the mesenchymal compartment. Whether or not this is the case will need further investigation. Epha1 inactivation in mice resulted in a rather mild phenotype in axial structures: kinked tails (Duffy et al., 2008). However, a full appreciation of the role of Epha1 in NMP biology might be obscured by possible functional redundancy with other ephrin receptors, most particularly Epha2, as it binds the same ligands as Epha1 and its inactivation also affects development of axial structures (Naruse-Nakajima et al., 2001).

Embryo staging was defined according to the standard timed mating approach, considering E0.5 the morning on which a mating plug was found. To isolate axial progenitors from developing embryos, matings were set up between Cdx2P-Cre^ERT transgenics (Jurberg et al., 2013) and ROSA26-YFP-Cre reporter (Srinivas et al., 2001) mice. Pregnant females were treated at E7.5 with a single intraperitoneal injection of 200 μl of 1 mg/ml tamoxifen (Sigma-Aldrich, T5648) in corn oil. Embryos were then collected at E10.5 by caesarean section, dissected in ice-cold phosphate-buffered saline (PBS) (Biowest, L0615) and processed for cell sorting (see below). Wild-type embryos for in situ hybridization or immunofluorescence were dissected in ice-cold PBS, fixed in 4% paraformaldehyde (PFA) (Sigma-Aldrich, P6148) in PBS and stored in methanol. Mice were genotyped as described by Aires et al. (2019), using primers listed in Table S4.

All animal procedures were performed in accordance with Portuguese (Portaria 1005/92) and European (directive 2010/63/EU) legislation and guidance on animal use in bioscience research. The project was reviewed and approved by the Ethics Committee of the Instituto Gulbenkian de Ciência and by the Portuguese National Entity Direcção Geral de Alimentação Veterinária (license 014308).

Whole-mount in situ hybridization was performed according to the protocol described by Aires et al. (2019). Probes for Efna1, Epha1, Ngfr, Cldn9, Nkd2 and Arl4d were prepared by amplifying cDNA fragments and cloning them into appropriate vectors for in vitro transcription. The sequences of all primers used to amplify these cDNAs are listed in Table S4. Whole-mount stained embryos were embedded in gelatin and sectioned using a vibratome following the protocol described by Dias et al. (2020). The β-galactosidase staining was performed as described by Aires et al. (2019).

Whole-mount immunostaining was carried out according to Dias et al. (2021). RapiClear 1.52 (SunJin Lab) was used for embryo clearing. Primary antibodies (used at 1:200 dilution) were as follows: rabbit anti-Sox2 (Abcam, AB92494) and goat anti-Epha1 antibody (R&D systems, AF3034). Secondary antibodies (all used at 1:1000 dilution) were as follows: donkey anti-goat 488 (Molecular Probes, A11055) and donkey anti-rabbit 568 (Molecular Probes, A10042). Images were acquired with a Prairie two-photon system, and the image dataset pre-processing was performed as described previously (Dias et al., 2021). No deconvolution or z-depth signal attenuation was performed in order to reduce possible pre-processing interferences that could misrepresent the real Epha1 expression levels in the tailbud.

Immunofluorescence staining of mouse sections and fixed sorted cells were performed as described previously (Aires et al., 2019). Primary antibodies (used at 1:200) were as follows: chicken anti-GFP (Abcam, AB13970), rabbit anti-Sox2 (Abcam, AB92494), mouse anti-Sox2 (Santa Cruz Biotechnology, sc-365823), rabbit anti-Tuj1 (Abcam, AB18207) and goat anti-Tbx6 (R&D, AF4744). Secondary antibodies (1:1000) were as follows: goat anti-chicken 488 (Thermo Fisher Scientific, A-11039), donkey anti-rabbit 568 (Molecular Probes, A10042) and donkey anti-goat 488 (Molecular Probes, A11055). Samples were analyzed using a Leica SP5 live or a Zeiss LSM 980 confocal microscope.

CJ7 mouse ES cells (Swiatek and Gridley, 1993) were maintained in ES cell medium [DMEM High Glucose (Biowest, S17532L0102), 15% defined fetal bovine serum (Hyclone, GE Healthcare, SH30070.03), 1% MEM non-essential amino acid solution (Sigma-Aldrich, M-7145), 2 mM L-glutamine, 1% EmbryoMax Nucleosides (Millipore, ES-008-D), 100 U/ml penicillin and 100 μg/ml streptomycin (Sigma-Aldrich, P7539), 0.1 mM β-mercaptoethanol and 1000 U/ml LIF (Millipore, ESG1107)] on mitomycin C-inactivated primary mouse embryo fibroblasts. To start differentiation, ES cells were divided into NMC-like populations using a protocol adapted from Gouti et al. (2014). Briefly, ES cells were removed from feeders by dissociation using 0.05% trypsin-EDTA solution (Sigma-Aldrich, 59417C) and seeded at a density of 5000 cells/cm² on CellBIND Surface dishes (Corning, 3294) in N2B27 medium [Dulbecco's Modified Eagle Medium/F12 (Gibco, 21331-020) and Neurobasal medium (Gibco, 21103-049) (1:1), 40 μg/ml BSA, 0.1 mM β-mercaptoethanol and supplemented with 1× N-2 (Gibco, LS17502048) and 1× B-27 minus vitamin A (Gibco, LS12587001)]. Cells were grown in N2B27 medium with 10 ng/ml bFgf (Peprotech, 100-18B) for 3 days (D1-D3). Neuromesodermal identity was induced by the addition of 5 μM CHIR99021 (Abcam, ab120890) from D2 to D3.

Wild-type mouse embryos around E10/10.5 were dissected in ice-cold M2 medium (Sigma-Aldrich, M7167) and their tailbuds, up until the hindlimbs, dissociated into single cells through mechanical pipetting in cold Accutase solution (Sigma-Aldrich, A6964). NMC populations from differentiated mouse ES cells were also dissociated to single cells following a similar protocol. After adding two volumes of PBS/10% donkey serum (DS) (Biowest, S2170), the single-cell suspension was centrifuged at 120 g and washed twice with PBS/10% DS. After that, the single cells were incubated in blocking solution (10% DS with 1:100 dilution of 2.4G2 anti-mouse Fc block in PBS) on ice for 30 min and then stained (also on ice) for 30 min with a 1:100 dilution of goat anti-Epha1 antibody. After two washes with PBS/10% DS, cells were incubated with a 1:1500 dilution of donkey anti-goat A647 antibody (Thermo Fisher Scientific, A-21447) for an additional 30 min on ice and subsequently washed twice with PBS/10% DS. Finally, the single cells were resuspended in PBS, filtered through a 100 μm cell strainer and sorted on a FACSAria IIu (BD Biosciences). Epha1-positive and -negative cells were identified using a 633 nm excitation laser with filter detection of 660/20 following the gating conditions illustrated in Fig. S5 and then further processed similarly. Cell sorting was performed using a purity mask of 16 to avoid sort non-target particles in our samples and ensure high purity (>80-90%) of the sorted samples. Sorting experiments were performed independently three times.

Epha1-positive and -negative FACS cells from embryonic tissue or from in vitro-generated NMC cells were collected in N2B27 medium supplemented with 10 ng/ml bFgf, 5 μM CHIR99021 and 10 μM Y-27632 inhibitor. To remove PBS added by the sorting machine, the single-cell suspension was centrifuged at 270 g and the culture medium replaced. Cells were plated in precoated 0.1% gelatin coverslips, at a density of around 3×10⁴/cm² and cultured at 37°C, with 5% CO₂, for 3 to 24 h (the best results for mesoderm differentiation were achieved with only 3 h). After this period, for differentiation in mesoderm, cells were cultured with N2B27+5-8 μM of CHIR99021. For neural differentiation, we used N2B27 medium supplemented with 100 nM of retinoic acid (Sigma, R2625) and 500 nM of Shh agonist SAG (Calbiochem, 364590-63-6). After 2 or 3.5 days of differentiation, cells were washed twice with PBS, fixed in 4% PFA for 20 min and stored in PBS at 4°C for immunohistochemistry.

Two regions were collected from E10.5 Cdx2P-Cre^ERT::ROSA26-YFP embryos: the tailbuds and a more-proximal region of the tail tip (up to the third somite). To obtain a single-cell suspension, tissue was incubated on ice for 15 min in Accutase. Digestion was terminated by adding two volumes of PBS/10% DS and washed twice with PBS/10% DS. Cells were then resuspended in PBS/2% DS and filtered through a 100 μm cell strainer. Cells were sorted according to their YFP-positive fluorescence in a MoFlo sorter (Beckman Coulter) using a 488 nm excitation laser with detector filter of 520/40. The YFP parameters were set using cells from Cdx2P-Cre^ERT tails dissected in parallel to serve as a YFP-negative control. YFP-positive cells were collected directly in TRI Reagent (Sigma-Aldrich, T9424) and stored at −80°C. The collected YFP-positive cells from the tail tip and the proximal tail were designated as Tail^Prog and Tail^Desc, respectively.

To isolate Epha1 cell populations for RNA-seq, single-cell suspensions were obtained from tailbuds and proximal tail regions of E10.5 wild-type embryos as mentioned above. Gating conditions of Epha1^High and Epha1^Low cell populations (illustrated in Fig. S3) were based on an apparent separation in total Epha1-positive cells histogram. For RNA-seq analyses, sorted cells were collected directly in TRI Reagent and kept at −80°C until further use. For all RNA-seq experiments, 15,000-20,000 purified-sorted cells per sample were collected to obtain a sufficient concentration of high-quality RNA.

Total RNA was isolated from the TRI Reagent suspension following the manufacturer's protocol, with the addition of 10 mg RNase-free glycogen (Roche, 10901393001) in the isopropanol step. RNA samples were then resuspended in RNase-free water. RNA concentration and purity were determined on an AATI Fragment Analyzer (Agilent). For the Tail^Desc and Tail^Prog samples, RNA-seq libraries were prepared from two biological replicates using TruSeq Stranded mRNA sample Prep Kit (Illumina, 20020594) and sequenced using Illumina HiSeq 2500 system at the CRG Genomics Unit (Barcelona, Spain). At least 25 million single end 50 bases reads were generated for each library. Read alignments were performed by TopHat2 v2.0.9 (Kim et al., 2013) with Bowtie2 v2.1.0.0. (Langmead and Salzberg, 2012). Differential expression analysis between Tail^Desc, Tail^Prog and Tail^Tot (Aires et al., 2019) was performed using CuffDiff v2.1.1 (Trapnell et al., 2013).

RNA-seq from Epha1^High and Epha1^Low cells was performed using two separate biological replicates. Libraries were prepared from total RNA using the SMART-Seq2 protocol (Picelli et al., 2014). Sequencing was performed on Illumina NextSeq500 at the IGC Genomics Facility, generating 20-25 million single-end 75 base reads per sample. Read alignments were performed as above. Read count normalization and differential expression between Epha1^High and Epha1^Low samples were analyzed using the DESeq2 R package (Love et al., 2014). To faithfully compare all samples (Tail^Contr, Tail^Prog, Epha1^High and Epha1^Low) from the above-mentioned RNA-seq independent experiments, normalization and differential expression were performed using the DESeq2 R package (Table S4). Representative heatmaps were created using BioVinci 2.0 data visualization software. The sequencing data of the RNA-seq experiments have been deposited in the NCBI trace and Short-read Archive (SRA) under accession numbers PRJNA527654 and PRJNA527619.

Quantitative RT-qPCR was carried out as described previously (Aires et al., 2019). The sequences of the primers used are given in Table S5.

Cells obtained from embryos and stained with the Epha1 antibody as described above were then washed twice with PBS/10% DS and processed for staining using True-Nuclear transcription factor Buffer Set (BioLegend, 424401) according to the manufacturer's instructions and following the protocol described previously (Aires et al., 2019). The following primary antibodies were used: rabbit anti-Sox2 (1:200; Abcam, AB92494) and mouse anti-Tbxt (1:100; Santa Cruz Biotechnology, sc-166962). Donkey anti-rabbit A568) (Molecular Probes, A10042; 1:1500) and donkey anti-mouse Alexa Fluor 488 (Abcam, ab150105; 1:1500) were used as secondary antibodies. For multicolor FACS analysis, single stain compensation controls were used to correct the spectral overlap between different fluorophores. Gating conditions of Epha1^High and Epha1^Low subsets were based on an apparent separation in a total Epha1-positive cells histogram (A647 channel). Quadrant gates were established according to fluorescence levels detected by the unstained controls and processed without primary antibodies; the negative control was determined by forelimbs cells exposed to the same conditions as the tail samples (Fig. S3). Flow cytometry data was analyzed using FlowJo 10 (BD Biosciences) software. Quadrant averages were calculated using at least three independent experiments and one-way analysis of variance ANOVA was used to determine statistical significance.

The public single-cell molecular map of mouse gastrulation and early organogenesis (Pijuan-Sala et al., 2019) was used to assess the gene expression of Epha1, Tbxt, Sox2 and Nkx1-2 in Theiler stage 12 wild-type mouse embryos. E8.5 CLE single cells were obtained from Gouti et al. (2017), following the procedure described previously (Dias et al., 2020) and using the single-cell consensus clustering (SC3) software (Kiselev et al., 2017). Gene expression visualization was obtained using SPRING (Weinreb et al., 2018), with the following parameters: minimum UMI total (for filtering cells)=1000; minimum number of cells with >=3 counts (for filtering genes) was set to 3; gene variability percentile (for filtering genes)=50; number of PCA dimensions (for building graph) was set to 20; number of nearest neighbors (for graph)=5.

We thank the IGC animal house, Flow Cytometry, Genomics and Imaging Facilities for their expert services, advice and assistance, and Daniel Sobral from the IGC Bioinformatics Unit for assistance with RNA-seq analysis; the SunJin laboratory for the RapiClear test sample; Val Wilson for sharing microarray processed data; and Anestis Tsakiridis for feedback on the culture and differentiation of FAC-sorted single cells. We also acknowledge all members of the Mallo laboratory for helpful discussions and comments throughout the course of the project, and Rita Aires and Ana Casaca for reading the manuscript.