A Gene Regulatory Network Balances Neural and Mesoderm Specification during Vertebrate Trunk Development

Abstract

Transcriptional networks, regulated by extracellular signals, control cell fate decisions and determine the size and composition of developing tissues. One example is the network controlling bipotent neuromesodermal progenitors (NMPs) that fuel embryo elongation by generating spinal cord and trunk mesoderm tissue. Here, we use single-cell transcriptomics to identify the molecular signature of NMPs and reverse engineer the mechanism that regulates their differentiation. Together with genetic perturbations, this reveals a transcriptional network that integrates opposing retinoic acid (RA) and Wnt signals to determine the rate at which cells enter and exit the NMP state. RA, produced by newly generated mesodermal cells, provides feedback that initiates NMP generation and induces neural differentiation, thereby coordinating the production of neural and mesodermal tissue. Together, the data define a regulatory network architecture that balances the generation of different cell types from bipotential progenitors in order to facilitate orderly axis elongation.

Keywords: NMPs; dynamical systems modeling; gene regulatory networks; neuromesodermal progenitors; retinoic acid; single-cell transcriptome analysis; vertebrate development.

Graphical abstract

Figure 1

Single-Cell Transcriptome Analysis of In Vivo NMPs Defines the Molecular Signature of e8.5 and e9.5 NMPs (A) Strategy for single-cell transcriptional analysis of e-NMPs dissected from the CLE region of e8.5 and e9.5 mouse embryos. (B) Hierarchical clustering of embryo-derived single cells, using the genes in the first three modules (Table S1), columns represent cells and rows correspond to genes. This separates cells into two large groups correlating with developmental age (e8.5 and e9.5). Within the e8.5 group, three smaller clusters could be distinguished, an e8.5 NMP identity (expressing genes in module 1), MPCs (expressing genes in modules 1 and 3) and mesodermal cells (Meso) (expressing genes in module 3). The e9.5 group was subdivided into two groups, one associated with e9.5 NMP identity (module 2) the other with MPC fate (modules 2 and 3). The cluster identity of each cell from the e8.5 and e9.5 embryos is indicated in orange (e8.5 NMPs), purple (e9.5 NMPs), brown and pink (MPCs), and red (Meso). (C) Pseudotemporal ordering of cells (right) obtained via the associated cell state graph (left) (expression levels are indicated as normalized counts per million reads). The white-to-red colors indicate NMP-to-mesodermal fractional identity of each cells defined by Sox2, Nkx1.2, Msgn1, Tbx6, and Meox1 levels. (D) Developmental trajectories of single cells from e8.5 and e9.5 mouse embryos reveals three distinct populations, NMP (T/Bra⁺/Sox2⁺), MPC (T/Bra⁺/Msgn1⁺/Tbx6⁺), and PSM (T/Bra⁻/Msgn1⁺/Tbx6⁺). (E) Molecular signature of e-NMP cells identified by differential expression analysis of the e8.5 (top) and e9.5 (bottom) NMPs with PSM cells. (F) The analysis identified 68 genes that were associated with both e8.5 and e9.5 NMP identity and 31 genes associated with mesodermal differentiation. For clarity we have shown only the 31 genes most enriched in e-NMPs. The complete gene lists are given in Table S2. (G) Differential expression analysis of e8.5 NMPs and e9.5 NMPs defines the molecular signature of e-NMPs at different developmental stages. Log CPM, logarithmic counts per million; Log f.c., logarithmic fold change; NMPs, neuromesodermal progenitors; MPC, mesodermal progenitors; PSM, presomitic mesoderm; SB, somite border.

Figure 2

Single-Cell Analysis of In Vitro Derived NMPs Indicates that They Resemble Their In Vivo Counterparts (A) Schematic of the differentiation conditions used for the generation of in vitro NMPs from mouse pluripotent stem cells. (B) Hierarchical clustering of all in vitro derived D3 cells using the genes contained in the three modules (identified from the in vivo population) reveals three distinct clusters. An e8.5 NMP identity characterized by genes in module 1, an e9.5 NMP identity characterized by genes in module 2, and a mesodermal characterized by genes in module 3. (C) More than 60% of D3 NMP cells had a profile similar to e8.5 NMP cells, and 10% had a profile similar to e9.5 NMP cells. (D) tSNE projection of e8.5, e9.5 and D3 NMPs using the genes comprising the three modules reveals that D3 cells appear more similar to e8.5 NMPs than e9.5 NMPs. (E) Pseudotemporal ordering of in vitro generated NMP cells at D3 identifies a similar developmental trajectory to in vivo derived cells. The cell state graphs were obtained using the 99-gene signature of e-NMP cells partially shown in Figure 1F. The white-to-red colors indicate NMP-to-mesodermal fractional identity of each cell.

Figure 3

Single-Cell Analysis of the Differentiation Route of NMPs toward the Neural Lineage (A) Schematic of the differentiation conditions used for the generation of in vitro NMPs and neural progenitor cells (NPCs) from mESCs. (B) Immunohistochemistry of cultures at D3 indicates that most cells co-express Brachyury (T/Bra) and Sox2, characteristic of NMPs, and a small percentage of cells expressed Tbx6. Removal of CHIR after D3 and culture until D4 in NB (neurobasal) conditions induces the generation of NPCs that express Sox2 in the absence of T/Bra, and a few mesodermal cells expressing Tbx6. (C) Quantitation of T/Bra⁺, Sox2⁺, Tbx6⁺, or T/Bra⁺/Sox2⁺ signal⁺ area normalized to DAPI area at D3 (bFGF/CHIR) and D4 (NB) of differentiation. Error bars indicate SD of four randomly selected independent fields. (D) Pseudotemporal ordering of D3 and D4 cells identifies four different populations and two distinct differentiation trajectories that lead to either neural or mesodermal identity. Expression of Cdx1, Cdx2, Cdx4, T/Bra, and Nkx1.2 is high in NMP cells, the expression of mesoderm specific genes Msgn1 and Tbx6 is high in Mesoderm cells. By contrast, there is induction of Sox1, Irx3, and Zic2 along the neural trajectory. The expression of RA signaling pathway components is differentially regulated in each developmental trajectory. RXRγ and RARγ are co-expressed in D3 NMP cells (similar to e8.5 NMPs). Expression of Aldh1a2 correlated strongly with Msgn1 and Tbx6 expression. See also Figure S1. (E) Hierarchical clustering of D3 and D4 single cells partitions the cells into four major groups. An NMP cluster (NMP, cerulean), a transitioning NMP (t-NMP, cyan) that express markers of developmentally older NMPs, a mesodermal (Meso, red), and a neural progenitor cell cluster (NPC, green). D3 cells are indicated in gray and D4 cells in black. See Tables S3 and S4. (F) Diagram illustrating the mesoderm or neural progenitor fate choice made by NMP cells. (G) Schematic of the posterior part of an e8.5 mouse embryo. NMP cells (cerulean) expressing Bra/Sox2/Cdx genes are located in the CLE region, close to the NSB in the anterior part of the primitive streak. As cells leave the NMP zone they differentiate to MPC progenitors (cerulean/red) expressing Bra/Msgn1/Tbx6, which results in upregulation of Aldh1a2. Thus, increased levels of RA produced in close proximity to the niche promote Sox2 expression and the differentiation of NMPs to PNT cells (green/yellow) then NPCs (green). The transcriptional network that controls the cell fate decision of NMP cells toward neural or mesodermal identities is summarized adjacent to the embryo model. NMP, neuromesodermal progenitor; t-NMP, transitioning neuromesodermal progenitor; NPCs, neural progenitor cells; PSM, presomitic mesoderm; MPC, mesodermal progenitor cells; PNT, pre-neural tube cells; NSB, node-streak border; PS, primitive streak; NB, neurobasal conditions.

Figure 4

Induction of NMPs Requires Low Levels of RA Signaling (A) Schematic of in vitro differentiation conditions and diagram illustrating the different cell fate choices of epiblast cells in the presence of bFGF/CHIR signaling and variable levels of RA. (B) Immunohistochemistry at D3 of differentiation indicates that Aldh1a2^−/− ESCs exposed to bFGF/CHIR downregulate Sox2 and express the mesodermal markers T/Bra and Tbx6. WT ESCs differentiated under the same conditions co-express T/Bra and Sox2. (C) Exposure of cells to increased levels of RA from D2 to D3 eliminates NMP induction and instead induces an NPC identity, evident by the expression of Sox2 in the absence of T/Bra and Tbx6. (D) qRT-PCR analysis of the expression of Sox2, T/Bra, Tbx6, Msgn1, Sox1, and Wnt3a at D3 in Aldh1a2^−/−cells and WT ESCs treated with bFGF/CHIR or bFGF/CHIR/RA (RA 10 or 100 nM). Mesodermal markers are induced in Aldh1a2^−/− cells, whereas the expression of Sox2 and Sox1 is abolished. By contrast, increasing RA concentrations induce neural fate identity, characterized by the upregulation of Sox2 and Sox1, whereas the expression of the mesodermal genes, T/Bra, Msgn1, and Tbx6 is absent. Expression of Wnt3a is significantly dowregulated under RA conditions. See also Figure S2. (E) Schematic of differentiation conditions used for Cdx^1,2,4−/−cells. (F) Immunohistochemistry at D3 indicates the induction of cells that co-express T/Bra and Sox2. Also, Tbx6 expression is initially induced in the Cdx^1,2,4−/−cells. (G) Although continuing exposure to CHIR results in WT cells predominantly adopting Tbx6-expressing mesodermal identity, Cdx^1,2,4−/−cells acquire a Sox2-expressing NPC identity. (H) Inhibition of RA signaling (with 1 μM BMS) partially restores mesodermal differentiation, revealed by the upregulation of T/Bra in the Cdx^1,2,4−/− cells. (I) qRT-PCR analysis of mesodermal genes T/Bra and Tbx6, RA signaling pathway components Cyp26a1, Aldh1a2, and Wnt and Fgf signaling ligands Wnt3a and Fgf8 in Cdx^1,2,4−/− and WT cells at D3 (NMP conditions), D5 (CHIR conditions), or D5 CHIR conditions with RA inhibition (1 μM BMS) from D3 to D5. In the Cdx^1,2,4−/−cells the expression of T/Bra is induced at D3 but at lower levels, and the expression of Aldh1a2 is substantially increased compared with WT ESCs. At D5, expression of mesodermal markers T/Bra and Tbx6, as well as Wnt3a and Fgf8, is downregulated in Cdx^1,2,4−/−cells. RA inhibition in the presence of Fgf signaling results in partial restoration of T/Bra and Fgf8 in the Cdx^1,2,4−/−cells. qRT-PCR data were normalized against β-actin. Error bars indicate SD of three biological replicates. See also Figure S3. (J) Quantitation of Bra⁺, Sox2⁺, Tbx6⁺, or Bra⁺/Sox2⁺ signal⁺ area normalized to DAPI area at D3 (bFGF/CHIR) and D5 (CHIR or CHIR/BMS) of differentiation. Error bars indicate SD of four randomly selected independent fields. PNT, pre-neural tube; PS, primitive streak.

Figure 5

NMP Cells Are Generated in the Absence of Msgn1 but Cannot Efficiently Differentiate to PSM Identity (A) Schematic of in vitro culture conditions used for assaying Msgn1^−/− cells. (B) Immunohistochemistry of Msgn1^−/− cells at D3 of differentiation reveals most cells co-express T/Bra and Sox2. The expression of Tbx6 is also evident in some cells at D3. (C) At D5 under mesodermal conditions WT ESCs predominantly differentiate to PSM, whereas Msgn1^−/− cells are maintained in an NMP state co-expressing T/Bra with Sox2. Tbx6 is also expressed in some cells, but these have not downregulated T/Bra. (D) Downregulation of T/Bra is delayed in the Msgn1^−/− cells under NB conditions, as is evident by the presence of T/Bra-expressing cells at D5. See also Figure S4. (E) Quantitation of T/Bra⁺, Sox2⁺, Tbx6⁺, or T/Bra⁺/Sox2⁺ signal⁺ area normalized to DAPI area at D3 (bFGF/CHIR) and D5 (CHIR) of Msgn1^−/− ESC differentiation. Error bars indicate SD of four randomly selected independent fields. (F) qRT-PCR analysis of T/Bra, Tbx6, Cdx2, Sox2 and Nkx1.2 at D3 (NMP conditions), D5-CHIR (mesodermal conditions) and D5 NB (neural conditions) in Msgn1^−/− and WT cells. At D3 Msgn1^−/− cells express high levels of T/Bra and lower levels of Tbx6 compared with controls. Cdx2 and Nkx1.2 expression is not affected, whereas Sox2 is downregulated. At D5 CHIR conditions Msgn1^−/− cells express high levels of NMP markers, T/Bra, Sox2, Cdx2, Nkx1.2 and lower levels of Tbx6, compared with WT cells that have acquired a PSM identity. (G) The expression of Wnt3a and Fgf8 is significantly higher in the Msgn1^−/− cells at D5 CHIR, whereas the expression of Aldh1a2 is downregulated. qRT-PCR data were normalized relative to β-actin. Error bars indicate SD of three biological replicates.

Figure 6

Maintenance of NMP Cell Identity in the Absence of Tbx6 (A) Schematic of conditions to assay Tbx6^−/− ESCs. (B) Tbx6^−/− cells co-expressing T/Bra and Sox2 at D3 under NMP conditions (bFGF/CHIR). (C) Tbx6^−/− cells are maintained as NMPs characterized by the co-expression of T/Bra and Sox2 at D5 (CHIR conditions), whereas WT cells mostly downregulate T/Bra and Sox2 and instead express Tbx6. (D) At D5 in CHIR conditions, Tbx6^−/− cells acquire an identity more similar to e9.5 NMPs characterized by the co-expression of T/Bra⁺/Sox2⁺ with Hoxc10⁺. In WT cells, few late NMP cells could be detected, as most Hoxc10-expressing cells were T/Bra negative. (E) qRT-PCR analysis of NMP, mesodermal, and neural markers at D3 (NMP conditions), D5 CHIR (mesodermal conditions), and D5 NB (neural conditions) shows that NMPs are induced at D3 in the Tbx6^−/−cells and maintained during exposure to CHIR. Expression of Aldh1a2 is upregulated at D3 and D5 (CHIR) in WT ESCs, whereas Tbx6^−/− cells express low levels of Aldh1a2 and higher levels of Fgf8 and Wnt3a at D5 (CHIR). (F and G) The Tbx6^−/− cells maintained NMP identity for 9 days (two passages) under bFGF/CHIR conditions (F) and progressively express more posterior Hox genes (G). (H) tSNE projection of the in vitro WT and Tbx6^−/− passaged NMPs with e8.5 and e9.5 NMPs revealed that D3 WT and Tbx6^−/− in vitro NMPs are similar to e8.5 NMPs, whereas the passaged D6 Tbx6^−/− in vitro NMPs closely resemble e9.5 NMPs. (I and J) Taking the intersection of differentially expressed genes identified between in vitro D3 Tbx6^−/− and D6 Tbx6^−/− cells, and in vivo e8.5 and e9.5 NMPs showed similar changes in gene expression of early (I) and late NMP signature genes (J).

Figure 7

Reverse Engineering the NMP Gene Regulatory Network (A) Six dynamical patterns summarize the experimental observations used as objectives to identify the best-fit network topology. Targeted time points are squares positioned at observed gene levels for Sox2 (orange), T/Bra (blue), and Msgn1/Tbx6 (red). The x axis represents the simulated time between D2 and D5. The y axis represents the protein concentration in a.u. All simulations are performed with a deterministic ordinary differential equation model from the initial conditions: Sox2^HIGH, T/Bra^LOW, and Msgn1/Tbx6^LOW. In addition to WT simulations, T/Bra and Msgn1/Tbx6 mutant simulations were performed by setting the relevant protein production rate to zero. Dashed lines represents the four RA/Wnt signaling conditions and solid lines a typical best-fit solution. (B) Overview of the 64 parallel parameter explorations evaluated. Each bar plot represents the average score of the best-fit solution cluster for one topology (a perfect score is 0). Bar colors identify the six objectives shown in (A). Underneath, circles describe the associated topology (blue circles for repression, red circles for activation) and the radius of each circle is proportional to the binding affinities in the best-fit solutions. (C) The best-fit topology identified by the parameter exploration. (D) Results of stochastic simulations. Colored pie charts show the ratios of cell states observed at D5 of stochastic simulations using the optimal topology under various RA/Wnt conditions between D3 and D5. Prior to D3, all trajectories are simulated with the same RA^LOW/Wnt^HIGH signaling condition. (E) Stochastic trajectories obtained with four RA/Wnt signaling conditions that are representative of the four conditions in (D). The thick solid lines indicate the average trajectories for neural (green) and mesodermal (red) fates. Circles represent the attractors of the dynamical system. (F) Experimental results. The ratios of different cell types obtained by co-measuring Sox2, T/Bra, and Tbx6 levels by flow cytometry. Assays were conducted using six RA/Wnt signaling conditions between D3 and D4. Grayscale panels show 2D kernel density estimations for each cell obtained from the protein fluorophore intensities. Contour plots document the highest density regions of the four cell populations identified by k-means clustering. Axes follow the “logicle” scale (Parks et al., 2006). Associated bar plots indicate the ratios of cell types in each condition.