Dynamic extrinsic pacing of the HOX clock in human axial progenitors controls motor neuron subtype specification

The patterning of bilaterian body is orchestrated by the differential expression of HOX transcription factors along their rostro-caudal axis. In vertebrates, HOX proteins are encoded by 39 genes, organized into four genomic clusters (HOXA, HOXB, HOXC and HOXD), that display anterior boundaries of expression co-linear to their chromosomal location. Genes located in the 3′ region of the clusters are expressed more anteriorly than their 5′ neighbors. The regionalized expression of HOX genes in post-occipital regions is initiated by their sequential, collinear, activation in axial progenitors, a dynamic cell population that fuels rostro-caudal extension of the embryo (Cambray and Wilson, 2002, 2007; Deschamps and Duboule, 2017; Forlani et al., 2003; Gouti et al., 2017; Henrique et al., 2015; Kondoh et al., 2016; Wymeersch et al., 2019). As axial progenitors contribute to progressively more caudal mesodermal and neuroectodermal structures, the 3′ to 5′ sequence of HOX gene activation is translated into a co-linear spatial pattern of expression in the progenies (Deschamps and Duboule, 2017; Henrique et al., 2015). Hence, the coupling between rostro-caudal extension of the body axis and the sequential activation of Hox genes over time, the HOX clock, is a central element of axial patterning. However, the mechanisms pacing the clock in axial progenitors are still unclear. Although extrinsic factors such as retinoic acid (RA), Wnt, fibroblast growth factors (FGFs) or growth differentiation factors (GDFs) control HOX patterns of expression, cell-intrinsic changes from a transcriptionally repressive to a permissive chromatin state occur within HOX complexes and correlate with the HOX temporal sequence of induction (Bel-Vialar et al., 2002; Deschamps and Duboule, 2017; Liu et al., 2001; Mazzoni et al., 2013; Narendra et al., 2015; Neijts et al., 2017; Noordermeer et al., 2014; Philippidou and Dasen, 2013; Soshnikova and Duboule, 2009). Whether the progressive opening of chromatin along the complexes serves as an internal timer actuated by extrinsic cues or whether sequences of secreted factors that activate progressively more caudal HOX genes largely defines the tempo of the induction remains unclear (Bel-Vialar et al., 2002; Del Corral and Storey, 2004; Deschamps and Duboule, 2017; Ebisuya and Briscoe, 2018; Lippmann et al., 2015; Mazzoni et al., 2013; Wymeersch et al., 2019). Inability to distinguish between these two hypotheses might underlie our limited ability to finely tune the rostro-caudal identity of caudal cell types derived from pluripotent stem cells [PSCs: embryonic stem cells (ESCs) and induced pluripotent stem cells (iPSCs); Diaz-Cuadros et al., 2020; Du et al., 2015; Duval et al., 2019; Faustino Martins et al., 2020; Frith et al., 2018; Gouti et al., 2014; Li et al., 2005; Lippmann et al., 2015; Matsuda et al., 2020; Maury et al., 2015; Ogura et al., 2018; Peljto et al., 2010; Sasai et al., 2014; Verrier et al., 2018]. The two models imply distinct strategies for in vitro differentiation: the ‘intrinsic model' predicts that the specification of posterior identities will require a precise synchronization between the timing of differentiation and the internal HOX timer. Alternatively, the ‘extrinsic model' implies that exposing axial progenitors to the relevant extrinsic cues should control the HOX clock to generate progenies of defined rostro-caudal identities. To approach this question and its consequences for cell engineering, we aimed to generate axial progenitors from human pluripotent stem cells (hPSCs) to study the mechanisms pacing the HOX clock during their differentiation into spinal motor neurons (MNs), a cell type that relies on a precise HOX code to acquire appropriate rostro-caudal subtype identities. Here, we report the generation of functional axial progenitors from hPSCs demonstrated by their ability to generate spinal MN subtypes located along the rostro-caudal axis of human embryos. As they become older, hPSC-derived axial progenitors generated MN subtypes born later in development and located more caudally in embryos, showing that they undergo a temporal shift in their rostro-caudal potential. Transcriptomic analysis indicated that this shift is linked to a temporal activation of HOX genes paralleled by a graded increase in FGF signaling. Blocking the FGF pathway stalled HOX temporal activation which resulted in the specification of early-born anterior MNs from late progenitors. In contrast, FGF precocious increase accelerated the temporal sequence of HOX activation and promoted the generation of more caudal, thoracic, MN subtypes. Combining FGFs with GDF11, another extrinsic factor, further accelerated the clock, favoring the specification of late-born lumbar MNs. Overall, our results indicate that extrinsic factors function as key pacers of the HOX clock tempo. Manipulations of the pacing pathways allow the efficient and synchronous engineering of human cell types of defined rostro-caudal identities for basic and translation applications.

The differential expression of Hox transcription factors along the rostro-caudal axis of the vertebrate spinal cord controls the formation of molecularly and functionally distinct neuronal populations, which therefore represent a functional product of Hox gene regulation. In particular, the Hox code orchestrates the acquisition of subtype-specific features in mouse and chick spinal motor neuron (MNs), controlling the formation of the distinct locomotor circuits (Dasen, 2017; Philippidou and Dasen, 2013). Whether this spinal HOX code and associated MN subtypes is conserved in human remained unknown, thus preventing faithful assessment of HOX regulation and its link with cell ‘rostro-caudal' identity during hPSC differentiation. We therefore mapped on human embryos, at 6.3 and 7.5 weeks of development, seven HOX transcription factors that display spinal collinear expression patterns and instruct MN subtype specification in mouse and chick (Philippidou and Dasen, 2013) (Fig. 1A, Fig. S1). As in animal models, human MNs expressed ISL1 or HB9 all along the spinal cord (Fig. 1A, circles in Fig. S1) (Amoroso et al., 2013). Within MNs, HOX displayed rostro-caudal patterns resembling those of mouse and chick (Dasen et al., 2003; Liu et al., 2001): cervical MNs expressed HOXA/C5, while brachial MNs expressed HOXC6, thoracic MNs expressed HOXC9 and lumbar MNs expressed HOXC10. Caudal brachial MNs co-expressed HOXC6 and HOXC8, and anterior thoracic MNs expressed HOXC8 and HOXC9. HOXD9 labeled caudal thoracic MNs as well as anterior lumbar MNs, which also expressed HOXC10 (Fig. 1A, Fig. S1). In amniotes, this Hox code instructs the formation of distinct motor columns that innervate common muscle groups, and motor pools that innervate a single muscle (Philippidou and Dasen, 2013). To be able to assess in vitro whether changes in HOX expression have functional consequences, we mapped columnar and pool markers within HOX expression domains (Fig. 1A, Figs S1 and S2). As in mouse and chick, MNs expressing high levels of FOXP1 (FOXP1^high) were observed at brachio-thoracic (HOXC6 and HOXC8/HOXC9) and lumbar (HOXC10) levels, where they formed a lateral motor column (LMC) matching the location of limb-innervating MNs (Amoroso et al., 2013; Routal and Pal, 1999) (Fig. 1A, Figs S1 and S2). Within the FOXP1^high LMC, MNs expressing the transcription factor SCIP (POU3F1) together with HOXC8 were observed in the caudal brachial spinal cord (Helmbacher et al., 2003). In addition, SCIP/HOXC8/HOXC9-positive MNs were detected in the anterior thoracic region, which might correspond to hand-innervating MNs (Fig. S2) (Bell et al., 2017; Mendelsohn et al., 2017).

Overall, HOX transcription factors are regionally expressed along the rostro-caudal axis of the human spinal cord. Within the HOX domains, distinct MN subtypes, identifiable by specific combination of transcription factors, are generated at stereotyped positions. These data provided readouts to assess the mechanisms regulating HOX expression in axial progenitors and their impact on cell type specification during hPSC differentiation.

We first sought to test whether axial progenitors (i.e. progenitors producing derivatives along the anteroposterior axis; Cambray and Wilson, 2002, 2007; Forlani et al., 2003; Kondoh and Takemoto, 2012) could be generated from hPSCs using MN subtype specification as a readout. We have previously reported sequences of extrinsic cues leading to the efficient generation of spinal MNs from hPSCs through a putative axial progenitor stage (Maury et al., 2015). Exposure of free-floating embryoid bodies (EBs) to a WNT pathway agonist (CHIR99021, CHIR) and inhibitors of BMP and TGFβ pathways (LDN193189 and SB431542) generated progenitors expressing CDX2, which encodes a transcription factor expressed in axial progenitors and a regulator of caudal HOX gene induction (Bel-Vialar et al., 2002; Bialecka et al., 2010; Gouti et al., 2014; Maury et al., 2015; Mazzoni et al., 2013; Neijts et al., 2017). Then, human spinal MN progenitors (83% of the cells on average) were generated upon exposure of these progenitors to: (1) retinoic acid (RA), which was previously shown to guide axial progenitors toward a neural fate; and (2) an agonist of the sonic hedgehog pathway (SAG), which induces a ventral MN progenitor fate (Briscoe and Novitch, 2008; Maury et al., 2015; Ribes et al., 2009). Finally, inhibition of Notch signaling on day 9 (D9) promoted the differentiation of MN progenitors in post-mitotic MNs (Maury et al., 2015). We thus used MN subtype markers, including HOX transcription factors, to assess the rostro-caudal identity of MNs produced under the previously reported conditions (Fig. 1B) (Maury et al., 2015). First, staining for HB9, ISL1 and the pan-neuronal marker neurofilament light chain (NEFL), together with quantification of ISL1⁺ cells, confirmed the efficient generation of spinal MNs as previously shown (Fig. 1C, Fig. S3A). Analysis of HOX expression showed that exposing cells to RA from D2 followed by SAG on D4 gave rise to some HOXC6 MNs corresponding to anterior brachial location in vivo. This identity was acquired by most MNs when RA was added 1 day later, on D3 (70.1% of HOXC6⁺ and 13.9% HOXC8⁺ MNs) (Fig. 1D,E). Addition of RA/SAG from D4 generated MNs with a caudal brachial identity (44.6% HOXC6/C8⁺) from which 37.1% expressed a high level of FOXP1, which identifies limb-innervating MNs in the spinal cord (Fig. 1D-F, Fig. S3B). Overall, these results suggested that the specification of more caudal MN subtypes was dependent on either: (1) the duration of exposure to RA, with a shorter RA exposure promoting caudalization; or (2) the time-point at which the progenitors received RA, as previously suggested (Bel-Vialar et al., 2002; Del Corral and Storey, 2004; Lippmann et al., 2015). To distinguish between these two possibilities, D3 progenitors that gave rise mostly to HOXC6 MNs received shorter RA treatments (D3-5 or D3-8 versus D3-9). These two conditions failed to promote more caudal identities (Fig. S3C-E). Instead, HOXC6 expression was lost in the shortest RA treatment (48 h, D3-5, Fig. S3D), indicating a potential role for RA in HOXC6 MNs specification as previously demonstrated in chick (Liu et al., 2001). These results indicate that the day at which progenitors are exposed to RA/SAG is the main trigger of caudalization (Fig. S3C-E). Further delaying RA/SAG then induced even more caudal MN subtypes. Treatment of cells with RA/SAG on D5 generated MNs with an anterior thoracic identity (59.8% HOXC8/C9⁺) (Fig. 1B-E) from which 27.6% MNs acquired a FOXP1⁺ limb-innervating identity (Fig. 1F). Accordingly, HOXC9/FOXP1/SCIP-positive MNs, which are located in the human anterior thoracic spinal cord and might correspond to hand-innervating MNs, were observed almost exclusively in this condition (Fig. 1F, Fig. S3G). Even later treatment, on D6 or D7, specified MNs acquiring a mid-thoracic identity, as demonstrated by the expression of HOXC9 and the loss of FOXP1^high MNs (Fig. 1B-F, Fig. S3F,G). The progressive caudalization of MN identity upon incremental delays of RA/SAG addition was confirmed with another hESC line and two iPSC lines (Fig. S4A,B). Interestingly, for one iPSC line, we had to increase CHIR concentration from 3 to 4 µM to observe this caudalization, suggesting line-to-line differences in WNT pathway integration (Fig. S4A-E). Indeed, and in agreement with the lack of caudalization, upon 3 µM of CHIR, this iPSC line exhibited a decreased induction of CDX2 at D3, a transcription factor induced in axial progenitors by WNT and FGF signaling pathways, and required for caudal HOX gene induction (Fig. S4F) (Amin et al., 2016; Bel-Vialar et al., 2002; Mazzoni et al., 2013; Nordström et al., 2006; Young et al., 2009).

Overall, WNT activation combined with TGF-β/BMP pathway inhibition converts hPSC into cells that can generate progenies of distinct rostro-caudal identities, a hallmark of axial progenitors. The time window between the initial exposure to WNT agonist and later exposure to RA defines HOX expression patterns in the generated progenies.

We then aimed at defining more precisely the identity of these hPSC derived-axial progenitors and the molecular changes paralleling their temporal switch in rostro-caudal potential, including HOX gene induction patterns. We first wondered whether the progenitors generated between D2 and D4 share a molecular identity with other vertebrate axial progenitors. In mouse, axial progenitor activity is carried in part by cells located in the caudal lateral epiblast (CLE) and the node-primitive streak border, which express Sox2, Nkx1.2, Cdx2 and TbxT (Brachyury) (Albors et al., 2018; Amin et al., 2016; Cambray and Wilson, 2002, 2007; Edri et al., 2019; Forlani et al., 2003; Gouti et al., 2014, 2017; Henrique et al., 2015; Koch et al., 2017; Kondoh et al., 2016; Mathis and Nicolas, 2003; Tzouanacou et al., 2009; Wymeersch et al., 2016, 2019). Single-cell lineage-tracing experiments showed that some of these axial progenitors are dual-fated neuromesodermal progenitors (NMPs) that generate progenies both in the spinal cord and paraxial mesoderm (Tzouanacou et al., 2009). The relative levels of Sox2 and TbxT transcription factors correlate with a preponderant differentiation of the cells toward the mesodermal (higher TbxT) or neural lineages (higher Sox2). We thus monitored the relative level of expression of SOX2, TBXT and NKX1.2 transcripts in D2, D3 and D4 progenitors in comparison with both hESCs and hiPSCs (Fig. S5A,B). The three markers were present from D2 to D4. SOX2 displayed levels similar to those in hPSCs, while NKX1.2 and TBXT were highly induced (Fig. S5A,B). We then stained for SOX2, TBXT and CDX2 (Fig. 2A, Fig. S5C). They were all co-expressed during the timeframe, yet TBXT progressively decreased over time to be weakly expressed by D4. This decrease in TBXT expression likely indicates a progressive transition towards preneural plate progenitors, which lose TbxT expression and are committed to the neural lineage (Del Corral and Storey, 2004; Delfino-Machín et al., 2005; Olivera-Martinez et al., 2014; Wymeersch et al., 2016).

Whole-transcriptome analysis confirmed that, beside these markers, D2, D3 and D4 hPSC-derived progenitors shared many similarities with mouse axial progenitors found in the caudal epiblast (Fig. 2, Fig. S5 and Tables S1, S2). Functional enrichment analysis of genes increased in progenitors in comparison with hESCs indicated an activation of HOX genes and the WNT pathway in D3 and D4 progenitors, accompanied by a decrease of key pluripotency markers, such as NANOG or KLF2 and KLF4, as observed in CLE cells (Fig. 2C, Fig. S5D-F) (Del Corral and Storey, 2004; Deschamps and Duboule, 2017; Olivera-Martinez and Storey, 2007; Olivera-Martinez et al., 2014; Wymeersch et al., 2019). More strikingly, the most common genes defining mouse CLE cells (Edri et al., 2019; Gouti et al., 2017; Olivera-Martinez et al., 2014; Wymeersch et al., 2019) were present from D2 to D4 (Fig. 2C, Table S1), and D2 and D3 progenitors expressed 122 out of 142 genes defining mouse E8.5 and E9.5 NMPs in comparison with mesodermal progenitors (RPM≥1.0, Fig. S5G, Table S4). Out of these 122 genes, the CLE genes CDX1, FGF17, RXRG, SP8, WNT5A, WNT8A and FST all belonged to the 50 most enriched genes in D2 and D3 progenitors compared with hESCs (Fig. S5H, Tables S1 and S4). There was also a strong overlap with genes found in recently reported hPSC-derived axial progenitors (Fig. S5J, Table S3) (Frith et al., 2018). Overall, the transcriptomic signature of the progenitors, together with the co-expression of SOX2, TBXT and CDX2, demonstrate the WNT-induced specification of human axial progenitors resembling mouse CLE progenitors (in which NMPs reside) (Edri et al., 2019; Gouti et al., 2014, 2017; Henrique et al., 2015; Koch et al., 2017; Liu et al., 2001; Wymeersch et al., 2016). However, the degree of similarity with the different states of mouse axial progenitors, including early and late NMPs and whether individual cells are bipotent, would require additional analysis and experiments. Furthermore, the progressive decrease of TBXT indicates a transition toward more neural potent progenitors, as observed in the mouse anterior CLE and in regions anterior to the node where pre-neural axial progenitors are located (Del Corral et al., 2003; Diez del Corral et al., 2002; Forlani et al., 2003; Wymeersch et al., 2016).

Next, we aimed to define the molecular changes paralleling the temporal change in rostro-caudal potential of the progenitors. A comparative analysis of axial progenitor transcriptomes indicated that D2 progenitors already activated most HOXB genes (except HOXB13), which are also activated early in development in mouse and chick and are not differentially expressed along the spinal cord (Dasen et al., 2003; Denans et al., 2015; Deschamps and van Nes, 2005) (Fig. 2D). Similarly, members of the HOXA, HOXC and HOXD complexes expressed early in development with anterior borders of expression in the hindbrain or anterior spinal cord were already activated at D2: HOXA4 and HOXA5 were expressed at low levels (Fig. 2D, Table S1) whereas HOXC4 and HOXC5 were below the detection level but detectable using real time PCR (see Fig. 3B,C). HOX genes regionally expressed along the brachial and thoracic spinal cord (Fig. 1A, Fig. S1) then displayed a temporal co-linear activation (Fig. 2D). At D2, the progenitors expressed low levels of HOXC6 that increased by D3; HOXC8 was expressed at higher levels than HOXC9 on D3 (Fig. 2D). The expression of HOXC8 and HOXC9 further increased by D4 (by 3.8- and 3.4-fold, respectively). In agreement, this progenitor state will give rise to HOXC8 and some HOXC9 MNs when exposed to RA/SAG (Fig. 1C,D). Non-expressed HOX genes corresponded to the 10 to 13 paralog groups – the latest and most caudally expressed genes (Gaunt, 1991; Izpisúa-Belmonte et al., 1991; Philippidou and Dasen, 2013; Fig. 1A, Fig. S1). Hence, WNT activation induced a temporal co-linear activation of HOX genes, which prefigured the generation of MNs of progressively more caudal identities upon incremental delays of RA/SAG treatment (Fig. 1).

To define the pathways activated in parallel to the sequential induction of HOX genes, we performed functional enrichment analysis on the upregulated genes between D2 and D3 when brachial and thoracic spinal HOX genes start to be induced [log2(FC)≥1.0, P-value≤0.05, 232 genes]. We detected a significant enrichment for target genes associated with an activation of the mitogen-activated protein kinase (MAPK) pathways (Fig. 2E). Indeed, typical target genes of the MAPK ERK1/2 pathway (such as DUSP4, DUSP6, SPRED2 and SPRY2) gradually increased from D2 to D4 (Fig. 2E,F), which was confirmed by real-time PCR in hiPSC-derived progenitors (Fig. S5K). MAPKs are classical mediators of FGF signaling (Lunn et al., 2007). In agreement, we observed a temporal increase of FGF8 and FGF17 expression, two secreted FGF ligands, previously described to increase over time in the caudal epiblast of chick and mouse embryos (Fig. 2F,G) (Liu et al., 2001; Wymeersch et al., 2019).

Hence, the sequential co-linear activation of HOX genes within axial progenitors is paralleled by an increase in FGF ligands and MAPK target genes, suggesting a temporal increase in FGF pathway activity. As FGFs modulate the expression of caudal Hox genes in other systems (Bel-Vialar et al., 2002; Dasen et al., 2003; Liu et al., 2001), we hypothesized that graded paracrine or autocrine FGF signaling might be triggering the sequential induction of HOX genes.

To test whether endogenous FGF signaling was necessary for the temporal sequence of HOX gene induction and caudal MN subtype specification, we first exposed progenitors to PD173074, a selective inhibitor of FGFR1/3, and monitored the impact on HOXA and HOXC genes, which undergo temporal activation in axial progenitors (Fig. 3A,B, Fig. S6A). Inhibition of FGFR1/3 for 24 h blocked the temporal increase of HOXC6, HOXA7, HOXC8 and HOXC9 but did not impact the expression of anteriorly expressed HOXA4, HOXA5 and HOXC5, which are not temporally activated in this time-frame (Fig. 3A,B). We then tested whether this effect was mediated by MEK1/2, a downstream effector of the FGFR pathway. Exposing D3 progenitors to PD0325901, a selective MEK1/2 inhibitor, for 24 h or 48 h also blocked the temporal increase of HOXC6, HOXA7, HOXC8 and HOXC9 (Fig. 3C, Fig. S6A). MEK inhibition even induced a drop in expression of all HOX genes, including HOXA4, HOXA5 and HOXC5, that was less pronounced at 48 h than at 24 h (Fig. 3C, Fig. S6B). This result suggests a stronger impact of MEK inhibition on HOX genes, potentially reflecting a more pronounced and rapid inhibition of pathway activity when targeting downstream effectors.

To determine the functional consequences of the stalled HOX temporal induction, we monitored MN subtype specification after FGFR1/3 and MEK1/2 inhibition (Fig. 3D-I, Fig. S6). First, the efficiency of MN generation was not impacted by the two inhibitors in the different tested conditions, even though we collected a reduced number of EBs in the D3-D7 PD173074 condition (Fig. S6D,F). Exposure to RA/SAG on D7 normally leads to the specification of HOXC8 and C9 MNs (Figs 1B-D and 3D,F,H). Addition of PD173074 or PD0325901 from D3 to D7 followed by RA/SAG at D7 fully prevented the generation of HOXC8/9 MNs (Fig. 3E,F). Instead, HOXC6 MNs were specified, an identity normally obtained when D3 progenitors are exposed to RA (Figs 1B,C and 3D-F). A similar result was obtained with only a 48 h inhibition from D3 to D5 followed by RA on D5, which also gave rise to HOXC8/9 MNs in absence of inhibitor (Fig. 3E,G-I, Fig. S6B,C). These results agreed with the view that FGF pathway inhibition stalled the progression of the HOX clock, resulting in the specification of MNs normally obtained from younger progenitors. To further verify this hypothesis, we inhibited the pathway 1 day later, at D4, followed by RA on D5. In contrast to D3 inhibition, this later inhibition did not prevent HOXC8 MN specification but blocked selectively HOXC9 MNs specification (Figs 1B and 3G-I). Altogether, these results demonstrate that autocrine and/or paracrine FGF signaling is necessary for HOX temporal induction within axial progenitors and the subsequent specification of caudal brachial and anterior thoracic MNs.

Our transcriptomic analysis indicated a temporal increase in FGF ligand gene expression, suggesting that an increase in FGF concentration and/or duration of exposure might be pacing the HOX clock. To test whether modulating the level of environmental FGF was sufficient to modify HOX induction, we exposed early D3 progenitors to RA/SAG, together with FGF2 or FGF8 at different concentrations and for different durations. In all conditions that received FGF2, more caudal MN identities were induced without the need for delaying RA addition (Fig. 4A-D, Fig. S7). In addition, the extent of caudalization varied with FGF2 concentration, with increasing concentrations promoting more caudal identities: caudal brachial HOXC8⁺ MNs were induced at 15 ng/ml (68.2%, 8.6-fold increase) and HOXC9⁺/HOXC6⁻ thoracic MNs at 60 ng/ml (58.9%, 825-fold increase), while 120 ng/ml further reduced HOXC6⁺ MNs (Fig. 4A-D, Fig. S7A-G). Similar results were obtained with FGF8 (Fig. S7B). FGFs acted directly and rapidly on axial progenitors, as addition of 120 ng/ml FGF2 for 24 h or 48 h induced MNs of a caudal brachial or mid-thoracic identity (after 24h, 49.7% of the MNs were HOXC6⁺, 79.4% HOXC8⁺, 14.7% HOXC9⁺ and 47.9% FOXP1^high/SCIP⁺. After 48 h, 7.5% were HOXC6⁺, 76.4% HOXC9⁺ and 13.2% FOXP1^high/SCIP⁺; Fig. 4A-C, S7C-I).

To determine whether this caudalizing effect was linked to an accelerated induction of brachial and thoracic HOX genes, we performed real-time PCR analysis 24 h and 48 h post-FGF2 treatment. We indeed observed an earlier increase of HOXC8, HOXC9 and HOXD9 mRNA expression compared with RA/SAG controls (Fig. 4F).

Hence, a precocious rise in FGF signaling accelerates the induction of caudal HOX genes in early axial progenitors, resulting in the specification of more caudal cell types within the same differentiation timeline (14 days). These results demonstrate that the levels or duration of FGF signaling dynamically regulate the pace of HOX co-linear activation in axial progenitors.

FGF2 accelerated HOX induction to generate MNs up to the mid-thoracic level, suggesting that early axial progenitors might lack competency to generate the most caudal segments. Alternatively, other extrinsic factors might be required to further accelerate the induction of HOX genes and promote more caudal identities. GDF11 is a member of the TGFβ family implicated in the control of axial elongation and MN subtype specification, and is required for the expression of the most 5′ HOX gene starting from the 10 paralogs in vivo and in vitro in late NMP-like cells (Aires et al., 2019; Gaunt et al., 2013; Lippmann et al., 2015; Liu, 2006; Liu et al., 2001; McPherron et al., 1999; Peljto et al., 2010). We exposed D3 progenitors to a combination of RA and GDF11 (25 ng/ml; Liu et al., 2001) for 24 h or 48 h. A 24 h GDF11 treatment induced caudal brachial and anterior thoracic MNs (64.7% HOXC8⁺ and 46.4% HOXC9⁺), the latter increasing after a 48 h treatment (Fig. 4A,B,E, Fig. S8A,B). However, as with FGF2, none of the more caudal identities was observed. In chick, exposure of spinal cord explants to combination of FGF2 and GDF11 promoted more caudal MNs than the two factors separately (Liu et al., 2001). We thus tested whether this combination might accelerate the HOX clock and promote caudal thoracic or lumbar identities. 24 h after exposure of D3 axial progenitors to FGF2 and GDF11, HOXC9, HOXD9 and HOXC10 mRNAs were strongly induced (respectively 65.4-, 2774.77- and 3329.44-fold compared with controls) and further increased after 48 h (Fig. 4F). In agreement, a 24 h (D3-D4) treatment combined with RA led to the specification of caudal thoracic MNs (73.9% HOXC9⁺ and 45% HOXD9⁺) whereas lumbar HOXC10⁺ MNs were observed after a 48 h FGF/GDF11 treatment (10.1% HOXC10⁺, 58% HOXD9⁺) (Fig. 4B,E, Fig. S8B-F). Importantly, adding FGF2/GDF11 prior RA (FGF/GDF11 at D3 followed by RA D5) largely increased the generation of lumbar MNs (81% HOXC10⁺, 89% HOXD9⁺) (Fig. S8G-I). This result suggests that RA antagonized the caudalizing activity of FGF2/GDF11, potentially by repressing FGF8 gene expression as previously shown (Del Corral and Storey, 2004; Del Corral et al., 2003) and/or by promoting a differentiation stage at which HOX expression is stalled. Hence, a combination of FGF2 and GDF11 further accelerates the HOX clock in young D3 progenitors to induce late and caudally expressed HOX10 genes, resulting in the generation of lumbar MNs.

In this study, using axial progenitors derived from hPSCs, we demonstrate that the parameters of exposure to two extrinsic factors, FGFs and GDF11, dynamically regulate the speed at which the HOX clock proceeds. Extrinsic control of the clock results in the synchronous specification of progenies of distinct rostro-caudal identities, which were born at different stages of embryogenesis. These results support the view that the pace of the HOX clock is largely regulated by extrinsic factors rather than by an intrinsic timer.

To reach these conclusions, we used 3D, embryoid body-based differentiation, of hPSC to selectively generate axial progenitors, an identity demonstrated by the co-expression of SOX2, TBXT and CDX2, as observed in mouse, chick and human embryos, complemented by detailed transcriptomic analysis and, more importantly, by functional analyses showing their ability to generate cell types found along the rostro-caudal axis of the human embryonic spinal cord (Figs 1, 2 and Figs S1-S5). Axial progenitors feed axial elongation in the embryo by generating diverse neural and mesodermal derivatives. They appear to form a heterogeneous dynamic population of progenitors comprising progenitors differentiating into diverse neural and mesodermal progenies (lateral, axial and paraxial), bi-fated NMPs that generate both spinal and paraxial mesoderm, as well as mono-fated progenitors that give rise either to mesodermal or spinal progenies (Albors et al., 2018; Forlani et al., 2003; Henrique et al., 2015; Tzouanacou et al., 2009; Wymeersch et al., 2016, 2019). Defining the precise degree of similarity between the reported hPSC-derived axial progenitors with these distinct states or with previously reported in vitro-generated axial progenitors and, in particular, whether individual progenitors are mono-fated or bi-fated NMPs will require further investigations (Wymeersch et al., 2021).

With efficient access to axial progenitors, we demonstrated that modulating the concentration, duration and combination of the caudalizing factors FGFs and GDF11 controls the speed at which the temporal activation of HOX genes occurs. The pace of the HOX clock is thus dynamically controlled by exposure parameters to extrinsic cues. Hence, the sequential changes in chromatin structure occurring along HOX complexes during their activation might be a consequence of changes in extrinsic signals rather than in the main mechanism intrinsically controlling the timing of HOX induction. (Bel-Vialar et al., 2002; Del Corral and Storey, 2004; Kimelman and Martin, 2012; Kmita and Duboule, 2003; Lippmann et al., 2015; Mazzoni et al., 2013; Narendra et al., 2015; Neijts and Deschamps, 2017; Noordermeer et al., 2011, 2014; Soshnikova and Duboule, 2009; Tschopp et al., 2009). As the sequence of activation from 3′ to 5′ genes was maintained under the different experimental conditions, the extrinsic cues seem to modulate the speed but not the directionality at which a repressive chromatin state is cleared from HOX clusters. This clearance could be temporally progressive or, alternatively, individual extrinsic factors might induce domain-wide, saltatory remodeling of repressive chromatin marks (up to HOXC9 in presence of FGFs, up to HOXD9/HOXC10 or more caudally in presence of FGF and GDF11) (Kmita and Duboule, 2003; Mazzoni et al., 2013; Narendra et al., 2016; Soshnikova and Duboule, 2009). Individual HOX genes within these transcriptionally permissive domains might be activated in a 3′ to 5′ gene direction due to their reliance on different combinations of transcription factors or differential strength of enhancer-promoter interactions (Kmita and Duboule, 2003; Mazzoni et al., 2013; Neijts and Deschamps, 2017). Determining the changes in histones marks, chromatin structure and the binding of signaling effectors within HOX complexes should help determine how extrinsic cues signal to the genome to control rostro-caudal patterning.

Dynamic environmental control of Hox expression was already suggested by grafting experiments in chick (Ensini et al., 1998; McGrew et al., 2008). In particular, it was demonstrated that heterochronic grafting of ‘old' axial progenitors to a ‘younger' caudal stem zone reverted their HOX profile to the ‘young' one (McGrew et al., 2008). Although our results might provide a potential molecular basis for environmental changes in HOX expression in caudal progenitors, we did not observe a reversion of the HOX profile upon FGF pathway inhibition but a stalled activation (Fig. 3E-I). This observation suggests that the FGF pathway regulates the tempo at which HOX activation proceeds rather than the maintenance of the HOX profile.

Importantly, the pacing of the HOX clock by secreted factors might ensure a community effect to synchronize HOX expression between neighboring progenitors, allowing the emergence of expression domains at tissue level (Durston, 2019). Furthermore, FGFs and GDFs are also implicated in embryo axial elongation (Aires et al., 2019; Boulet and Capecchi, 2012; Jurberg et al., 2013; Mallo et al., 2009; McPherron et al., 1999). A common mechanism to control rostro-caudal extension of the body axis together with Hox gene induction would be a parsimonious way to couple morphogenesis and patterning (Denans et al., 2015; Young et al., 2009). Hence, determining the mechanisms controlling the temporality of FGF and GDF expression onset, levels and duration should provide a better understanding of post-occipital tissues development and evolution.

In a bioengineering perspective, the extrinsic control of the HOX clock provides a means with which to manipulate HOX expression for the synchronous and efficient engineering of neuronal subtypes found at distinct rostro-caudal identities. In particular, we provide evidence for the controlled generation of putative hand- and leg-controlling MNs (Mendelsohn et al., 2017). Importantly, MN subtypes display differential vulnerabilities in MN diseases. The targeted differentiations in MN subtypes reported here should thus help in the more accurate study of these incurable diseases or allow the access of more controlled sources of cells for putative cell therapy approaches (Abati et al., 2019; An et al., 2019; Baloh et al., 2018; Nijssen et al., 2017; Ragagnin et al., 2019; Sances et al., 2016; Steinbeck and Studer, 2015; Tung et al., 2019). Considering the role of HOX transcription factors in instructing cell diversity in other derivatives of axial progenitors, such as control over their expression, might have broader applications for cell engineering beside spinal cell types (Deschamps and Duboule, 2017; Frith et al., 2018; Helmbacher et al., 2003; Iimura et al., 2009). More broadly, the temporal generation of distinct types of neurons or glia from the same progenitor domain is a widely used strategy to increase cell diversity in the nervous system (Dias et al., 2014; Kohwi and Doe, 2013; Oberst et al., 2019a; Rossi et al., 2017). Extrinsic cues play important roles in the unfolding of these temporal sequences (Kawaguchi, 2019; Oberst et al., 2019a,b; Syed et al., 2017; Tiberi et al., 2012). Extrinsic manipulation of the temporality of these lineages might help improve the generation of early and late-born cells for both basic research, disease modeling and cell therapy.

Human fetal embryos at 6.3 weeks of gestation were obtained from pregnant women referred to the Department of Gynecology and Obstetrics at the Antoine Béclère hospital (Clamart, France) for legally induced abortions in the first trimester of pregnancy, as previously described (Lambrot et al., 2006). All women provided written informed consent for scientific use of the fetal tissues. None of the abortions was due to fetal abnormality. The fetal age was determined by measuring the length of limbs and feet (Evtouchenko et al., 1996). The project was approved by the local Medical Ethics Committee and by the French Biomedicine Agency (reference number PFS 12-002). Alternatively, human embryonic spinal cords (n=2) at stage 7.5 weeks of gestation were collected in accordance with the national guidelines of the USA (National Institutes of Health, US Food and Drug Administration) and the State of New York and under Columbia University institutionally approved ethical guidelines relating to anonymous tissue. The material was obtained after elective abortions and was classified on the basis of external morphology according to the Carnegie stages. Gestational age was determined by the date of the last menstrual period of the patient or by ultrasound, if the ultrasound estimated date differed by 1 week (as indicated by the obstetrician). In all cases, the spinal cord was removed as intactly as possible before fixation with fresh, ice-cold 4% PFA for 1.5 h, washed abundantly with PBS and then cryoprotected overnight in 30% sucrose. Post-fixation, the cord was measured and cut into anatomical sections to accommodate embedding in OCT compound (Leica) and stored at −80°C before cutting on a cryostat. Sections (16 µm) were cut along the full length of the cord.

Human SA001 embryonic stem cell (ESC) line (male, RRID: CVCL_B347) was obtained from Cellectis and H9 ESC line (female, RRID: CVCL_9773) was obtained from Wicell. Both cell lines were used according to the French current legislation for hESC (Agency of Biomedicine, authorization number AFSB1530532S). The induced pluripotent stem cell (iPSC) line WTSIi002 (male, RRID: CVCL_AH30, alternative name HPSI0913i-eika_2) was obtained from the European Bank for Pluripotent Stem Cells (EBISC). WTC-mEGFP-Safe harbor locus (AAVS1)-cl6 produced by the Allen Cell Institute was obtained from Coriell (AICS-0036-006, male, RRID: CVCL_JM19). Experiments with iPSCs were approved by relevant ethic committees (declaration DC-2015-2559). All PSC lines were cultured at 37°C on Matrigel (Corning) in mTSER1 medium (Stem Cell Technologies) and amplified using EDTA (Life Technologies) clump-based passaging. They were tested for potential mycoplasma contamination every other week (MycoAlertTM Mycoplasma Detection Kit, Lonza, LT07-118). No contamination was detected during the study. PSCs were thawed in presence of Y-27632 (10 μM, Stemgent or Stem Cell Technologies) and the culture medium was changed every day.

Human PSC embryoid body-based differentiation was performed as previously described (Maury et al., 2015). hPSC were dissociated with cold Accutase (Life Technologies) for 3-5 min at 37°C and resuspended in differentiation medium N2B27 [Advanced DMEM F12, Neurobasal vol:vol (Life Technologies)], supplemented with N2 (Life Technologies), B27 without Vitamin A (Life Technologies), penicillin/streptomycin 1%, β-mercaptoethanol 0.1% (Life Technologies), Y-27632 (10 μM, Stemgent or Stem Cell Technologies), CHIR-99021 (3 µM or 4 µM Selleckchem), SB431542 (20 μM, Selleckchem) and LDN 193189 (0.1 μM, Selleckchem). Cells were seeded in ultra-low attachment 6 well plates (Corning) (2×10⁵ cells ml⁻¹) to form embryoid bodies (EBs). All conditions of differentiation received the same medium, without Y-27632, at day 2 and at day 3, but SB431542 was removed at day 3. The differentiation then proceeded according to the schematics presented in Fig. 1B, Fig. 3A,D,G, Fig. 4A. SAG (smoothened agonist, Merck millipore), FGF2 (recombinant human FGF basic, Peprotech), RA (Sigma-Aldrich), GDF11 (recombinant human/murine/rat GDF11, Peprotech), PD0325901 (Selleckchem), PD173074 (Selleckchem) and DAPT (Stemgent or Stem Cell Technologies) were added at indicated time points. For concentrations used, see Table S5. Media were changed every other day unless specified.

EBs were collected, rinsed with PBS then fixed with 4% PFA for 5 min at 4°C and rinsed with PBS three times for 5 min. EBs were cryoprotected with 30% sucrose and embedded in OCT (Leica) prior to sectioning with a cryostat. Alternatively, day 2, 3 and 4 progenitors were plated on Matrigel (Corning, diluted according to the manufacturer's recommendation) -coated coverslips and allowed to adhere between 30 and 60 min prior to fixation with 4% PFA for 5 min at 4°C.

All immunostainings were performed as follows: cells or sections were incubated with a saturation solution (PBS, 10% FBS, 0.2% Triton) for 10 min. Primary antibodies (Table S5) were diluted in a staining solution (PBS, 2% FBS, 0.2% Triton) and incubated overnight at 4°C in a humidified chamber. After four PBS washes (10 min each), secondary antibodies (Alexa488, Alexa555 and Alexa647, Life Technologies, 1:1000) were added for 1 h at room temperature. After three PBS washes, DAPI was added to cells (Invitrogen, 1:3000) for 5 min. Cells or slices were then mounted in Fluoromount (Sigma-Aldrich or Cliniscience).

Samples were visualized and imaged using either a Zeiss LSM 880 Confocal Laser Scanning Microscope controlled by Zen black software (Zeiss), a confocal microscope TCS SP5 II (Leica) or a DM6000 microscope (Leica) equipped with CoolSNAP EZ CDD camera, controlled by MetaMorph software. Alternatively, images were acquired using the automated microscope Cell Discoverer 7 (Zeiss), equipped with an Axiocam 506m camera, with Zen black software (Zeiss).

Total RNA was extracted (RNeasy Plus Mini Kit, Qiagen) and cDNA synthesized using SuperScript III (Invitrogen). Quantitative real-time PCR was performed using a 7900HT fast real-time PCR system (Applied Biosystems) with Sybr Green PCR Master Mix (Applied Biosystems) or performed using QuantStudio 5 Real-Time PCR System (Thermofisher Scientific) and a mix with qPCR Brilliant II SYBR MM with low ROX (Agilent). Primers are listed in Table S5. All expression data were normalized to cyclophilin A mRNA (ΔCt). All analyses were performed with three technical replicates per plate. In Fig. 3 and Fig. S5, mRNA expression levels are represented as a percentage of the highest expressed gene (ΔCt_max) among all the conditions [% relative to max gene expression=(ΔCt_max/ΔCt)×100]. In Fig. 4F, expression levels are represented as a percentage relative to the maximal expression level of each gene (ΔCt_geneMax) among all conditions [% relative to max gene expression=(ΔCt_geneMax/ΔCt)×100]. In Fig. S4, data are presented as relative expression levels determined by calculating 2^−ΔΔCt.

hESCs were exposed to CHIR-99021 together with LDN193189 and SB431542 while forming EBs. Progenitors were collected on days 2, 3 and 4 of differentiation, and processed for total RNA preparation. For each of the eight samples, total RNA was reverse transcribed using the Ion AmpliSeq Transcriptome Human Gene Expression kit (Thermofisher Scientific). The cDNA libraries were amplified and barcoded using the Ion AmpliSeq Transcriptome Human Gene Expression core panel and Ion Xpress Barcode Adapter (Thermofisher Scientific). The amplicons were quantified using Agilent High Sensitivity DNA kit before the samples were pooled in sets of eight. Emulsion PCR and enrichment was performed on the Ion OT2 system Instrument using the Ion PI Hi-Q OT2 200 kit (Thermofisher Scientific). Samples were loaded on an Ion PI v3 Chip and sequenced on the Ion Proton System using Ion PI Hi-Q sequencing 200 kit chemistry (200 bp read length; Thermofisher Scientific). The Ion Proton reads (FASTQ files) were imported into the RNA-seq pipeline of Partek Flow software (v6 Partek) using hg19 as a reference genome. The number of reads per sample ranged from 7.5 million to 12 million reads. To determine genes that are differentially expressed between groups, mapped reads were quantified using Partek E/M algorithm normalized by the Total count/sample (the resulting counts represent the gene expression levels on reads/millions for over 20,800 different genes present in the AmpliSeq Human Gene Expression panel). Genes with an average of reads ≤1.0 in all the time points were excluded from the analysis. The evaluation of the differential expression between two conditions was performed using the EdgeR package under R. Pathway enrichment analyses were performed on upregulated genes (FC⩾2.0, P-value<0.05) between two time points by interrogating Reactome database. Hierarchical clustering was performed using Ward's method with the ward.d2 algorithm. The Euclidean distance is represented. Significant enrichments were calculated using a hyper-geometrical test and Benjamini-Hochberg correction for multiple comparisons. The enrichment score was calculated as described by Wang et al. (2017).

The normalized transcriptomic data are provided in Table S1. The GEO accession number for the raw data is GSE153519. Lists of genes differentially expressed between time points are available in Table S2. Lists of common genes between D3 enriched genes and NMPs-like genes have been identified by Frith et al. (2018) and are available in Table S3. Expression of the genes in our model that have been identified in NMP cells at E8.5 and E9.5 by Gouti et al. (2017) are listed in Table S4.

All statistics were computed using Graphpad Prism software. One-way analysis of variance (ANOVA) with a Kruskal-Wallis post-hoc analysis was performed following normality tests provided by Prism. Number of samples (n), dispersion measures and P-values are indicated in figure legends. In all figures, n are independent differentiations started from independent newly thawed hPSCs vials. For each condition at least four independent EB sections were imaged in which all the cells were quantified by automated image analysis (four EBs were quantified in one replicate of the FGFRi 3-7 conditions; in all the other experiments 6 to 9 EBs were quantified). The images were exported and saved as TIFF with Fiji if needed (Schindelin et al., 2012). Quantitative analyses on images were performed using the CellProfiler software (Carpenter et al., 2006) (Broad Institute open source at www.cellprofiler.org). DAPI-stained nuclei were segmented into primary objects using the CellProfiler segmentation pipeline and the nuclear mask was used to define objects on the target channels. The threshold for defining positive nuclei for a given target was obtained using EB sections negative for the target of interest. All images across conditions were then automatically analyzed in batch to ensure unbiased analysis. The percentage of colocalization between markers (ISL1/HOXC6/C8, ISL1/C8/C9, ISL1/C9/C10) was obtained after co-immunostaining for the markers of interest followed by quantification in Cell profiler. The analyses of the FOXP1 and SCIP immunostaining intensities were performed by combining CellProfiler with the software FCS express 7 (DeNovo Software). Nuclei were segmented into primary objects as described above, and FOXP1 and SCIP fluorescence intensities were calculated for each primary object. Fluorescence intensity plots for FOXP1 and SCIP were then generated by FCS express 7 software to visualize the intensity levels of the different markers for each individual cell and determine the percentage of cells above a given threshold. Cell profiler pipelines for quantification are available upon request. For real-time PCR analysis, one-way analysis of variance (ANOVA) with a Kruskal-Wallis post-hoc analysis and correction of multiple comparisons by controlling the FDR (Benjamini and Hochberg) was performed following normality tests provided by Prism. Number of samples (n), dispersion measures and P-values are indicated in figure legends.

We thank Frédéric Causeret for providing sections of mouse embryos; Susan Morton, Thomas Jessell and Jeremy Dasen for the generous gift of antibodies; and F. Causeret, P. Gilardi, A. Rebsam, F. Giudicelli and E. Mazzoni for critical reading of the manuscript. Differentiation were performed at the ‘Cell engineering facility' of the Institut du Fer à Moulin (IFM) and at I-STEM. Imaging experiments were carried out at the Imaging facility of the Institut du Fer à Moulin (IFM) and at I-STEM.

The peer review history is available online at https://dev.biologists.org/lookup/doi/10.1242/dev.194514.reviewer-comments.pdf