Reconstruction of phrenic neuron identity in embryonic stem cell-derived motor neurons

Land vertebrates, including humans, use lungs to breathe air. The inspiratory and expiratory movements of the lungs are driven by a complex neural circuitry that consists of a central network in the brainstem that generates breathing rhythms and an output layer of motor neurons (MNs) that connect to respiratory muscles. These respiratory circuits develop prenatally and must become functional immediately after birth. Although significant progress has been made in understanding the central pattern generator itself (Champagnat et al., 2009; Fortin and Thoby-Brisson, 2009), relatively little is known about the formation of MNs that relay breathing rhythms from the CNS to the periphery. In mammals, respiration is driven by muscles that connect to the rib cage and thereby indirectly inflate and deflate the lungs. Arguably, the most important of these muscles is the diaphragm, which forms the boundary between the thoracic and abdominal cavities and contracts during inspiration. The diaphragm is innervated by the phrenic nucleus (PN), a population of MNs located in the mid-cervical spinal cord. During embryonic development, phrenic neurons emerge alongside other MNs from ventral progenitors (Arber et al., 1999; Briscoe et al., 2000; Thaler et al., 1999), send their axons through cervical ventral roots and then project caudally through the thoracic cavity to innervate the diaphragm muscle (Allan and Greer, 1997). We have only a partial understanding of the molecular cascade that establishes phrenic MN identity. It is, however, important that this pathway is defined as this would provide insights into how mammal-specific anatomical adaptations are patterned and allow us to model aspects of respiratory motor circuitry in neuronal cultures to study neuromuscular function and disease.

Spinal MNs segregate into distinct columns during embryogenesis. Each column connects to a specific set of muscles: the medial motor column (MMC) projects to epaxial muscles, the lateral motor column (LMC) innervates limb muscles (Lance-Jones and Landmesser, 1980; Tsuchida et al., 1994), and the hypaxial motor column (HMC) innervates body wall muscles (Dasen et al., 2003; Peljto and Wichterle, 2011). MNs acquire subtype transcriptional identities due to exposure to locally restricted morphogens (Marshall et al., 1992). For example, the expression of the MMC determinants Lhx3 and Lhx4 is sustained by floor plate-derived Wnt4 and Wnt5 (Agalliu et al., 2009). Likewise, brachial LMC fate depends on overlapping, segmentally restricted gradients of retinoic acid and Fgfs (Liu et al., 2001), which induce Hox6 paralogues and the accessory factor Foxp1 in presumptive LMC neurons in register with forelimbs (Dasen et al., 2008; Dasen et al., 2003). The HMC, by contrast, appears to lack specific determinants and might represent a ground state of MNs. Phrenic neurons are thought to be derivatives of the HMC (Rousso et al., 2008).

Some candidate determinants for early phrenic development have been identified, although their contribution to phrenic neuron specification is poorly understood: the transcription factor (TF) Pou3f1 is enriched in phrenic neurons, and Pou3f1 deficiency in mice leads to disorganisation of the PN (Bermingham et al., 1996). Absence of mid-cervical Hox5 paralogues affects the maintenance of phrenic neurons, but not their initial specification (Philippidou et al., 2012). Lastly, Foxp1 appears to negatively regulate the phrenic MN lineage, since Foxp1 mutants have increased numbers of phrenic MNs (Rousso et al., 2008). What is largely lacking at this point is an understanding of which potential effector genes are downstream of these factors, how these and other determinants interact and, in the case of Pou3f1, how the expression of the factor itself is initiated.

MN development can be recapitulated in vitro from mouse or human embryonic stem cells (ESCs), which will form functional spinal MNs under the appropriate culture conditions (Li et al., 2008; Miles et al., 2004; Wichterle et al., 2002). ESC-derivation of MNs depends on the same extrinsic and intrinsic cues that act during normal embryogenesis and has been repeatedly used to investigate subtype-specific developmental pathways in these cells (Jung et al., 2010; Peljto et al., 2010; Soundararajan et al., 2006). We set out to apply this approach to the acquisition of phrenic neuron identity.

To address how phrenic neuron fate is established in the developing spinal cord, we first identified candidate determinants in primary MNs sorted from mouse embryos, and then used a systematic in vitro gain-of-function (GOF) screening approach to test whether any given candidate approximates phrenic neuron transcriptional patterns when ectopically expressed in ESC-derived MNs (ESC-MNs). The aim was to define modules of effector genes downstream of the key determinants, as well as to understand how the determinants interact with each other. We found that the TFs Pou3f1, Hoxa5 and Notch intracellular domain (NICD) combine to regulate distinct sets of effector genes, which together comprise a large fraction of all phrenic neuron-specific genes. Moreover, expression of the receptors Cdh10, which is downstream of Pou3f1, or Pcdh10, a gene coordinately regulated by Hoxa5 and Pou3f1, is sufficient to mediate like-like clustering of MNs into aggregates in vitro, mimicking nucleus formation in vivo. Our findings suggest that local Notch-Delta interaction in the ventral spinal cord might be an early event in phrenic neuron specification and that a defined combination of intrinsic and/or extrinsic factors may be used to emulate phrenic neuron transcriptional identity and morphological features in MNs derived from pluripotent stem cells in vitro.

In the ventral embryonic spinal cord, phrenic neurons are generated alongside other MNs from ventral progenitors. To determine the transcriptional profile of embryonic phrenic neurons, and how they differ from other MNs, we isolated phrenic neurons from E11.5 Hb9::GFP mouse embryos. The genetic reporter labels axons and cell bodies of all spinal MNs (Wichterle et al., 2002). In cervical and thoracic trunk explants, the cell bodies of phrenic neurons were retrogradely labelled with the tracer TMR-dextran (Fig. 1A,B). Then, spinal cord segments C3-C5 were excised from the trunks, dissociated, and phrenic neurons were sorted as GFP⁺ TMR⁺ cells by flow cytometry (Fig. 1C). Immunohistochemistry confirmed that only Pou3f1⁺ phrenic MNs in the medial-dorsal part of the ventral horn were TMR positive (Fig. 1B). As a control population, non-phrenic MNs were isolated from the same cell suspension. Given their segmental origin, these MNs are likely to represent a mix of MMC, HMC and LMC neurons (Peljto et al., 2010). A second control population of pure LMC neurons was isolated by retrograde labelling through the radial nerve (Fig. 1D,E). Following isolation by flow cytometry, genome-wide transcriptional patterns of phrenic neurons and the two control populations of MNs were determined using Affymetrix arrays.

We next devised a scheme to systematically catalogue genes into groups based on their expression pattern in radial LMC MNs, phrenic MNs and non-phrenic control MNs derived from segments C3-C5, respectively. A logical value was attributed to each gene in each cell population: genes expressed significantly above the mean across the three MN populations received the value +1, genes expressed significantly below the mean -1, and genes with expression levels above the 5th and below the 95th percentile around the mean received the value 0. Hoxa1, Cdh10 and Klf5 are examples of phrenic neuron-enriched genes with an expression value of -1 in radial MNs, +1 in phrenic MNs, and 0 in control MNs (Fig. 2A). This pattern can be described with the row vector [-1 +1 0].

A number of known MN subtype markers found in different groups of neurons suggest that our approach yields genuine expression patterns (Table 1). The MMC determinants Lhx3 and Lhx4 follow the pattern [-1 -1 +1] and are enriched in the C3-C5 control (Table 1), the only cell population with an MMC component. Pou3f1, one of the few established phrenic neuron markers (Bermingham et al., 1996; Rousso et al., 2008), is found both in phrenic neurons and C3-C5 control MNs. Subsequent expression analysis revealed that, at E11.5, Pou3f1 is also expressed in a subset of MMC neurons (Fig. 2B; supplementary material Figs S3, S4), which explains the pattern that we observed. A second factor required for PN development, Hoxa5 (Philippidou et al., 2012), is enriched in phrenic neurons as well. We were able to confirm phrenic neuron-specific expression for several candidate genes by in situ hybridisation (Fig. 2C-U). These genes fall into different functional categories and include cadherins (Fig. 2C-F), Ig superfamily receptors (Fig. 2G,H), BMP inhibitors (Fig. 2K,L), synaptic proteins (Fig. 2P) and TFs (Fig. 2R-U). The identification of phrenic MN-specific cadherins is significant, as adhesion molecules of this protein family are expressed in specific combinations on many, if not all, MN subsets. Cadherins mediate like-like clustering of MN cell bodies and drive the formation of topographic maps in the ventral spinal cord (Bello et al., 2012; Demireva et al., 2011; Price et al., 2002). In summary, we were able to identify a large number of genes potentially enriched in primary phrenic neurons and confirmed most of those analysed by histology.

To assess the efficacy of these genes to assign PN identity, we adapted the ESC-MN culture system such that we could perform GOF screens for multiple factors in purified MNs. We chose a mouse ESC clone that allows the efficient insertion of DOX-inducible transgenes by Cre/loxP-mediated cassette exchange (Iacovino et al., 2011) and equipped it with an additional MN-specific reporter transgene that allows magnetic enrichment of this cell type (Fig. 3A). During the initial phase of the project, we used ESCs carrying an Hb9::CD2-GFP MN-specific reporter (supplementary material Fig. S5, Table S2). To optimise the sorting procedure, we generated an ESC subclone that carries the Hb9::CD14-IRES-GFP transgene (supplementary material Fig. S1), which resulted in an improved signal-to-background ratio (Fig. 3B) and sorting efficiency. The GFP MN reporter tightly correlates with the endogenous embryonic MN marker Isl1/2 (Fig. 3C,D), and, following MACS, yielded large numbers (10⁵-10⁶) of viable, enriched ESC-MNs (Fig. 3B,E-G). Furthermore, we were able to induce the expression of candidate transgenes with DOX in culture, as shown here for the candidate factor Pou3f1 (Fig. 3H,I).

We next examined the subtype composition of parental ESC-MNs to determine the baseline identity in our assay. Primary mid-cervical MNs belong to three distinct subsets: HMC, MMC and PN (supplementary material Fig. S3A-G). When we compared the phenotypes of ESC-MNs with in vivo MNs, we found that they mostly belonged to the Pou3f1^- Lhx3^- HMC and the Pou3f1^- Lhx3⁺ or Pou3f1⁺ Lhx3⁺ MMC subsets, with only a small percentage of Pou3f1⁺ Lhx3^- phrenic-like neurons (supplementary material Fig. S3H-L), consistent with an earlier study (Peljto et al., 2010). The percentage of ESC-MNs with an MMC phenotype declines between day 5 and day 6 (supplementary material Fig. S3L). This might be explained by the fact that all MN progenitors express Lhx3 initially, but only MNs committed to the MMC fate sustain its expression (Sharma et al., 1998). The MN subtype composition does not depend on the ESC clone used in this study, as we have detected similar ratios of MN subtypes in cultures differentiated from a second, independently derived ESC clone (supplementary material Fig. S3M).

Our analysis of genes enriched in primary phrenic neurons pinpointed a cluster of known and putative TFs that would be good candidates as determinants of PN identity (Table 1; supplementary material Table S1). In cases in which two closely related genes were isolated, only one gene was investigated further. Our transcriptional profile also identified Notch signalling as a potential key player in this process, as we found that Notch targets, such as Hey1 and Hey2 (Table 1), are enriched in phrenic neurons. Furthermore, we included a dominant-negative mutant form of Lhx3 (DNLhx3) (West et al., 2004) because DNLhx3 suppresses Lhx3 target genes and might inhibit this key MMC determinant (Agalliu et al., 2009). Vwc2 was the only non-TF chosen for analysis, because it is a chordin-like BMP inhibitor and may modulate neural patterning. Finally, we selected Hoxc6 and wild-type Lhx3 as control TFs, as they represent known LMC and MMC determinants, respectively. We established 19 sets of ESC subclones carrying single candidate transgenes (supplementary material Table S2) and confirmed inducible candidate gene expression in embryoid bodies (EBs) (supplementary material Fig. S6).

To address whether candidate determinants promote transcriptional profiles associated with phrenic neurons when ectopically expressed in ESC-MNs, we differentiated transgenic ESCs into MNs in EB cultures. MNs were magnetically sorted from dissociated EBs on day 5 and induced to express the transgene with DOX during the subsequent 30 hour culture period. In some experiments with Notch constructs, the ESC-MNs were only induced for the last 10 hours to limit the effect of the transgenes. RNA of ESC-MNs derived from two independent subclones for each candidate gene was isolated, and gene expression profiles determined using Affymetrix arrays. The expression data are represented as a scatter plot in Fig. 4, in which the number of PN-enriched genes (supplementary material Table S4) repressed/induced (y-axis) is plotted against all genes repressed/induced by a given candidate factor (x-axis). The mean value for each transgene is represented by two data points, one for repression (lower left quadrant) and one for induction (upper right quadrant). We would like to point out that the terms repression and induction do not necessarily imply a direct genetic interaction, as we measure mRNA levels and not promoter binding.

Thus, the lower the y/x ratio for repression (R-ratio) and the higher the y/x ratio for induction (I-ratio) the more the transcriptional pattern evoked by a transgene approximates that of primary phrenic neurons. The P-values obtained using a hypergeometric distribution for over-representation of induced or repressed genes (supplementary material Table S5) provide a very similar picture to that offered by the R- and I-ratios, except that these results are more difficult to visualise.

Of all candidate determinants tested, ESC-MNs overexpressing Pou3f1 showed the most specific induction of PN-enriched genes (I-ratio=0.3; Fig. 4, upper right quadrant), whereas not a single specific gene is repressed (R-ratio=0). The genes upregulated by Pou3f1 include the confirmed PN-specific genes Cdh9, Cdh10 and Pcdh11x (Fig. 2C,D,F). NICD, when DOX-induced for the last 10 hours of ESC-MN culture, was identified as a second PN-candidate factor with a high I-ratio (0.102) and low R-ratio (0.05) (Fig. 4; supplementary material Fig. S7A). The set of PN-specific genes upregulated by iNICD includes Vwc2l (Fig. 2K) and is largely non-overlapping with that observed in inducible (i) Pou3f1 ESC-MNs (supplementary material Fig. S8A, Table S6), although a few target genes, such as Hs6st2, are shared. The timing and signal intensity of Notch activity appear to be critical, as iNICD expressed for 30 hours, or the attenuated mutant isoform iNERT (Schroeder et al., 2003) induced for 10 hours or 30 hours, do not evoke a PN-like pattern. iHoxc6, despite the absence of Hoxc6 from primary phrenic neurons (Table 1), has a high I-ratio (0.109), yet also a high R-ratio (0.163) (supplementary material Fig. S7B). However, although iHoxc6 does induce 20 PN-enriched genes, including Hoxa5 (supplementary material Table S6), it also upregulates two key LMC determinants, Foxp1 (Dasen et al., 2008) and Aldh1a2 (Sockanathan and Jessell, 1998) (supplementary material Table S7). Ectopic iHoxa5 expression in ESC-MNs, by contrast, does not induce these LMC markers, and positively regulates only three phrenic neuron-enriched genes, including Ptn (Fig. 2J; supplementary material Table S6). Thus, when activated in isolation, iHoxa5 does not emulate a PN transcriptional programme. Finally, iDNLhx3 does not evoke PN-like patterns (Fig. 4), but surprisingly shows the best I-ratio for LMC-like patterns (supplementary material Fig. S9, Tables S5, S7). This suggests that repression of Lhx3 targets unlocks parts of the LMC but not of the PN transcriptional programme in ESC-MNs.

To explore whether any of the candidate determinants interact to elicit PN-like transcriptional patterns, we combined the most promising candidate, Pou3f1, with Hoxa5, Hoxc6, Hoxa1 and DNLhx3 in an expression system that allows us to simultaneously induce two candidate genes in the same ESC-MNs (Bondue et al., 2011). In addition, we combined Pou3f1 with Lhx3 to investigate negative interactions, as we observed that Pou3f1 target genes, such as Cdh9 and Cdh10 (Fig. 2C,D), are restricted to the PN, whereas Pou3f1 itself is also expressed by a subpopulation of the MMC (supplementary material Figs S3, S4) (Rousso et al., 2008). Array analysis of double-transgenic ESC-MNs revealed that iPou3f1/iHoxa5 and iPou3f1/iHoxc6 DOX induction resulted in nearly identical I-ratios (0.158 versus 0.163; Fig. 4; supplementary material Fig. S7C). By contrast, the R-ratio for iPou3f1/iHoxa5 (0.032) indicates higher specificity than that for iPou3f1/iHoxc6 (0.069). As in the iHoxc6 single-transgenic ESC-MNs, iPou3f1/iHoxc6 expression upregulates the LMC determinants Foxp1 and Aldh1a2, whereas iPou3f1/iHoxa5 expression does not (supplementary material Table S7). Interestingly, the number of PN-specific genes activated above the twofold threshold in iPou3f1/iHoxa5 ESC-MNs (22 genes) is higher than the sum of the two target gene sets in single transgenic MNs (12+3 genes), which suggests that several PN target genes, such as Pcdh10, require the combined activity of both determinants (supplementary material Table S6). Comparison of transcriptional patterns of iPou3f1 and iPou3f1/iLhx3 ESC-MNs shows that four of 12 iPou3f1-induced PN-specific targets, including Cdh9, are downregulated more than twofold by co-expression of Lhx3 (supplementary material Fig. S7D, Table S6). These findings provide evidence that Pou3f1 can interact with other determinants to regulate the activation of PN-specific target genes, and that both positive (Hoxa5) and negative (Lhx3) combinatorial interactions can be observed in ESC-MNs. Furthermore, Hoxa5 and Hoxc6, despite considerable overlap in their target genes (supplementary material Fig. S8B), differ in that Hoxc6 initiates LMC-like transcriptional patterns in ESC-MNs, whereas Hoxa5 does not.

In order to validate the transcriptional patterns observed by array analysis, we tested the expression of selected candidate genes in iYFP, iHoxa5, iNICD (10 hour pulse), iPou3f1 and iPou3f1/iHoxa5 ESC-MNs by qRT-PCR (supplementary material Table S8). The qRT-PCR expression data are largely consistent with the array data, although some additional target genes were detected. Crucially, the qRT-PCR analysis showed that NICD upregulates the PN determinant Pou3f1.

Do the three positive phrenic determinants we identified simply act in parallel or is there a transcriptional hierarchy, i.e. does one of the candidate factors initiate phrenic neuron specification? As the qRT-PCR data suggested that NICD activity upregulates Pou3f1, we further examined a possible hierarchical relationship between the two factors. We compared EBs derived from the parental ESC line with those derived from iNERT ESCs, which express an attenuated, tamoxifen (4HT)-inducible version of NICD (Fig. 5A-F). Even in the absence of induction with DOX and 4HT, iNERT EBs contained a significantly higher proportion of both phrenic-like Pou3f1⁺ Lhx3^- MNs and MMC-like Pou3f1⁺ Lhx3⁺ MNs (Fig. 5G), without a change in total MN numbers. This increase is likely to be caused by low-level transgene expression/activity in the absence of added inducers, which is a common phenomenon seen in inducible genetic systems (Howe et al., 1995; Schroeder et al., 2003). Full induction of the iNERT transgene led to a decrease in the number of MNs (Fig. 5E,F), consistent with previous findings when NICD was overexpressed in spinal progenitors in vivo (Dias et al., 2012). We next examined whether Notch signalling is required for Pou3f1 expression in a subset of normal ESC-MNs. We treated developing EBs with γ-secretase inhibitors, which block ligand-dependent Notch activation (De Strooper et al., 1999). Both of the inhibitors that we tested decreased the percentage of Pou3f1⁺ ESC-MNs by more than 90% (Fig. 5H-N), whereas they did not affect MN development in general. To test if Notch GOF upregulates additional PN markers, we examined the effect of a 10 hour pulse of DOX on the transcriptional patterns of iNICD EBs and cultured ESC-MNs, compared with normal controls. We found that NICD induces the PN-enriched genes Synpr (supplementary material Fig. S10), Cdh9, Cdh10, Pcdh11x and Ptn (supplementary material Table S8).

Given that Notch activity controls Pou3f1 in MNs in vitro, does the presence of a Notch ligand correlate with the appearance of Pou3f1⁺ MNs in vivo? Analysis of mid-cervical ventral spinal cord revealed that at E10.5, which is when MNs emerge from the pMN domain, a cluster of Dll4-positive cells is located just dorsal of nascent Pou3f1⁺ Isl1/2⁺ MNs (Fig. 6A,B). Previous studies have mapped Dll4 expression to the p2 progenitor domain (Del Barrio et al., 2007; Peng et al., 2007). Furthermore, Notch1 protein is detectable on the most immature, medial MNs, and Pou3f1 expression in MNs overlaps with that of the Notch target gene Hey1 (supplementary material Fig. S11). At E11.5, Pou3f1⁺ MNs have migrated laterally to settle in the PN and the MMC (Fig. 6D), a pattern that remains unchanged at E12.5 (Fig. 6G). Dll4 expression, by contrast, is transient and not seen at later time points (Fig. 6E,H), except in blood vessels. The expression of the Pou3f1 target gene Cdh10 is delayed compared with Pou3f1 itself and is first detected in the PN at E11.5 (Fig. 6C,F,I). The juxtaposition of Pou3f1 and Dll4 in vivo is mirrored by the pattern that we see in EBs, where patches of Pou3f1⁺ Isl1/2⁺ MNs are localised close to Dll4⁺ non-MNs (Fig. 6J,K).

To examine whether Dll4 protein is capable of inducing Pou3f1 expression in primary MNs, we used whole embryo culture (WEC) (Osumi and Inoue, 2001) to ectopically express this Notch ligand in cervical spinal cords of intact mouse embryos. Owing to the timing of the WEC cultures, many non-MNs in close proximity to Isl1/2⁺ MNs and ventral progenitors express the co-transfected reporter nGFP (Fig. 6L,N), whereas MNs themselves are mostly negative. Hence, any effect seen in MNs is likely to be non-cell-autonomous. Embryos developed normally in culture and formed a Pou3f1⁺ PN (Fig. 6L-O). In DLL4-transfected embryos (Fig. 6L,M), but not in controls (Fig. 6N,O), the transfected side of the spinal cord contains significantly more MNs of phrenic phenotype (Fig. 6P). At the same time, the total number of MNs is unchanged (Fig. 6P), suggesting that ectopic expression of DLL4 does not expand MN progenitors. Taken together, these findings show that Pou3f1 induction in embryonic MNs is controlled by Notch signalling both in ESC-MNs and in vivo, and suggest that Delta-Notch interaction is involved in the acquisition of phrenic neuron identity.

A key event in early MN patterning is the aggregation of MNs of the same subtype into clusters/nuclei within the ventral horn. This like-like recognition of cell bodies is mediated by receptors of the cadherin family. Expression profiling of primary MNs suggested that the cadherins Cdh10 and Pcdh10 are highly enriched in the PN (Fig. 2D,E). To examine whether these potential effector genes are sufficient to drive specific MN clustering, we ectopically expressed Cdh10 or Pcdh10 in ESC-MNs and mixed them with RFP-labelled control ESC-MNs (Fig. 7A-C). Near-neighbour analysis revealed that iCdh10 MNs and iPcdh10 MNs, but not RFP^- control MNs, segregated from RFP⁺ control MNs in culture and formed nucleus-like aggregates (Fig. 7D-J).

The development of phrenic neurons, which drive respiration in mammals, requires the integration of positional signals to restrict their specification to mid-cervical segments of the spinal cord, where they emerge in register with myoblasts destined to form their target muscle, the diaphragm (Babiuk et al., 2003). We performed expression profiling in primary embryonic neurons followed by the reconstruction of phrenic-like transcriptional patterns in transgenic ESC-MNs to identify factors that combine to establish PN identity along the dorsoventral and anteroposterior axes of the spinal cord. We found that three factors, Pou3f1, Hoxa5 and Notch, interact to control sets of putative PN-specific effector genes: Hoxa5 only evokes PN-like transcriptional patterns when combined with Pou3f1, and Notch signalling is upstream of Pou3f1. The latter interaction was observed both in ESC-MNs and in spinal cords of intact mouse embryos. Thus, our study suggests that Delta-Notch signalling at the dorsal margin of the nascent motor columns has to intersect with Hox5-dependent mid-cervical segmental identity for phrenic neurons to develop (Fig. 8). Given that NICD controls Pou3f1, why is there only a limited overlap between their target genes in ESC-MNs (supplementary material Fig. S8A)? The most likely explanation for this is the short time frame of iNICD induction (10 hours), which is sufficient to upregulate Pou3f1 itself (supplementary material Table S8) but not second-order Pou3f1 targets. Based on co-expression assays in ESC-MNs, we also predict that Lhx3 is a negative regulator of phrenic identity in Pou3f1⁺ MMC neurons. In addition to transcriptional determinants, we identified and validated several sets of putative PN effector genes downstream of these regulators of PN identity. Many of these genes are regulated in a modular fashion by single factors, whereas others require the combined activity of two determinants.

The aggregation of MNs that connect to the same target muscle into nuclei and motor pools is a key event in early MN development and depends on MN subtype-specific expression of cadherins (Bello et al., 2012; Demireva et al., 2011; Price et al., 2002). This establishes a topographic map of MN cell bodies within the ventral horn and ensures that functionally related MNs receive the same sensory synaptic input (Sürmeli et al., 2011), as well as facilitating the formation of electrical synapses between MNs to coordinate activity (Personius et al., 2007). We found that the phrenic neuron-enriched receptors Cdh10, which is regulated by Pou3f1, and Pcdh10 (Hirano et al., 1999), which is regulated by Hoxa5 and Pou3f1, mediate like-like clustering of ESC-MN cell bodies in vitro (Fig. 7, Fig. 8C) and may contribute to the aggregation of phrenic neurons into the PN, probably in concert with other adhesion molecules.

Our findings on PN specification in ESC-MNs suggest that the coordinated expression of cadherins and other effector genes in phrenic neurons depends on the combinatorial activity of at least three positive and one negative determinant. Nevertheless, our understanding of PN gene regulation remains incomplete. For example, we do not know how the PN-enriched receptor genes Alcam and Plxnc1 are regulated. There are a number of possible explanations for the limitations to our approach: (1) additional, unknown transcriptional determinants might positively regulate additional sets of PN-selective genes; (2) ESC-MNs might contain more negative regulators, which mask aspects of induced PN identity; (3) peripheral cues encountered by outgrowing phrenic axons might control some PN target genes, and these cues would not be present in vitro; (4) some of the PN-enriched genes that are not yet confirmed by histology could be false positives, in particular those in the least stringent [0 +1 0] group (supplementary material Table S4).

The two cervical Hox genes Hoxa5 and Hoxc6 evoked strikingly similar transcriptional patterns when overexpressed in ESC-MNs as single transgenes (supplementary material Fig. S8B) or combined with Pou3f1 (supplementary material Table S7), consistent with the fact that different cervical Hox proteins can bind to the same consensus DNA site (Pellerin et al., 1994). Co-expression of iPou3f1 plus iHoxa5 revealed combinatorial effects, and the overall PN specificity of the pattern is higher than that of iPou3f1 plus iHoxc6. Crucially, iHoxc6 and iPou3f1 plus iHoxc6, but not iHoxa5 and iPou3f1 plus iHoxa5, induce the TF Foxp1. This key LMC determinant, which largely depends on Hox6 expression in vivo (Dasen et al., 2008; Lacombe et al., 2013), suppresses PN identity when ectopically expressed in cervical MNs (Rousso et al., 2008). Thus, the combination of iPou3f1 plus iHoxc6 induces an LMC-type programme that is similar to that of the PN, yet at the same time is incompatible with PN identity. An in vivo correlate to the Foxp1-positive iPou3f1/iHoxc6 ESC-MNs indeed exists: ulnar motor pools, which are caudal LMC subpopulations that innervate distal limb muscles in mice (Lacombe et al., 2013), and anatomically related flexor carpi ulnaris MNs in chick (Dasen et al., 2005) co-express Pou3f1 and Hoxc6 during its initial specification. Despite their different identities, aspects of the subtype-specific genetic programmes are shared between the PN and ulnar MNs, as both populations selectively express the Pou3f1 target genes Cdh9 and Cdh10 (Fig. 2B-D; supplementary material Fig. S12).

MN subtype identity has been linked to Notch signalling: in Gde2 (Gdpd5) null mutant mice, which have increased levels of activated Notch, MNs show a marked shift in columnar subtype composition (Sabharwal et al., 2011). Intriguingly, in their study, Notch affects MN fate at the level of progenitor cells. We first observe Pou3f1 upregulation in Isl1/2⁺ postmitotic MNs (Fig. 6A). Nevertheless, we do not know whether Notch induction of this key PN determinant occurs before or after cell cycle exit. Delta-Notch interaction as a mechanism to diversify V2a/V2b spinal interneurons appears to be conserved between fish, birds and mammals (Batista et al., 2008; Del Barrio et al., 2007; Peng et al., 2007). There are two scenarios of how the ancestral Dll4⁺ territory in the p2 progenitor domain, which drives this process, might relate to the emergence of the PN during mammalian evolution. In the first scenario, Dll4 expression might have acquired a novel role in establishing the mammal-specific PN. This would imply that phrenic neurons arose from cervical HMC neurons, which were respecified by Delta-Notch signalling, as the diaphragm evolved from cervical hypaxial muscle (Perry et al., 2010). Alternatively, the PN might be an evolutionary modification of an ancestral LMC subset that predates mammals. Hirasawa and Kuratani recently proposed that the diaphragm originated from limb muscle associated with the pectoral girdle (Hirasawa and Kuratani, 2013). If this scenario is correct, then phrenic neurons might have evolved from a Pou3f1⁺ LMC subset that lost its brachial Hox code due to a duplication of two cervical segments coupled with a caudal shift in brachial identity, an event that might have occurred during early mammalian phylogeny (Hirasawa and Kuratani, 2013). This view is supported by the fact that, in mammals, there are two Pou3f1⁺ MN subsets spaced by about two segments at cervical levels, one rostral of the LMC (the PN) and one caudal within the LMC (ulnar MNs).

In this study, we have demonstrated that MN derivation from ESCs is a viable screening tool to systematically dissect developmental pathways. Our experimental strategy was based on initially identifying candidate determinants in sorted primary embryonic neurons by transcriptional profiling and then validating the factors in the ESC-derived equivalent of these ex vivo cells. Both a highly reductionist approach (genome-wide expression analysis in ESC-MNs) and the more complex, heterogeneous EB culture system provided insights into how a network of TFs interact to establish PN identity. Similar ESC-based screens could be undertaken for other cell types that are difficult to access in large numbers. We think that the derivation of MNs from pluripotent stem cells will mature into a powerful tool to investigate not only prenatal developmental programmes, but also subtype-dependent features of adult MN function and connectivity, complementing existing animal models. Our long-term aim is to assemble adult-like phrenic ESC-MNs, as well as other defined cell types, into artificial neuromuscular circuits to simulate aspects of respiratory motor function in health and disease in vitro. The progress towards understanding embryonic phrenic neuron identity reported here will allow us to extend the investigation to later steps of neuronal maturation and circuit formation.

Embryonic phrenic neurons were labelled by injecting a retrograde tracer (TMR-dextran) into the phrenic nerve in E11.5 trunk explants derived from mouse embryos carrying the Hb9::GFP transgene (Wichterle et al., 2002). Explants were cultured in an oxygenated bath for 2 hours and then dissociated with papain, and TMR⁺ GFP⁺ phrenic neurons, as well as GFP⁺ control MNs, isolated by flow cytometry (Fig. 1A-C). In separate experiments, radial LMC neurons were retrogradely labelled and isolated using the same methodology. Gene expression profiles of sorted primary MNs were determined by Affymetrix array analysis.

MNs were directly differentiated in vitro from mouse ESCs as described (Peljto et al., 2010). Briefly, A2.lox ESCs (Iacovino et al., 2011) carrying the MN-specific, magnetically sortable reporter genes Hb9::CD2GFP or Hb9::CD14-IRES-GFP (supplementary material Fig. S1) were grown for 2 days as EBs, induced with 1 μM retinoic acid and 0.5 μM smoothened agonist (SAG), and dissociated on day 5. MNs were then isolated by magnetically activated cell sorting (MACS). MACS-purified ESC-MNs were cultured for 30 hours on Matrigel, and the expression of PN-specific candidate genes (supplementary material Tables S1, S2) was induced with doxycycline (DOX). Then, transcriptional profiles of ESC-MNs were determined by Affymetrix array analysis and qRT-PCR (Spandidos et al., 2010) (supplementary material Table S3). ESC-MNs expressing the neutral transgene YFP were used as the baseline, as the DOX-inducible Cre recombinase present in the parental ESC line triggers a DNA damage response in MNs (supplementary material Fig. S2).

Detailed protocols are available in supplementary Materials and Methods.

We thank Tom Jessell and Andrew Lumsden for support and encouragement; Noriko Osumi and Yoko Arai for advice on WECs; Ira Schieren, Kaity Miao and Fazal Oozeer for expert technical assistance; Jeremy Dasen, Fiona Watt, Robert Knight, Linda Greensmith, Christopher E. Henderson, QueeLim Ch’ng, Stephen Price and Anthony Graham for comments on the manuscript; Chris G. Tan for sharing unpublished observations; and Austin Smith, Esteban O. Mazzoni, Hynek Wichterle, Koichi Kawakami, Cedric Blainpain, Dafe A. Uwanogho, Shinji Hirano, Ursula Just, Malcolm Logan, Linda Greensmith, Matthew Grist, Uwe Drescher and Simon J. Rhodes for providing reagents.