Single cell transcriptome dynamics from pluripotency to FLK1⁺ mesoderm

Previous studies have established that blood and endothelial cells differentiate from FLK1 (KDR)-expressing mesoderm (Ema et al., 2006; Kataoka et al., 1997; Kattman et al., 2006; Lugus et al., 2009). Flk1⁺mesoderm can also contribute to smooth muscle, cardiac lineages and paraxial mesoderm tissues (Ema et al., 2006; Kataoka et al., 1997; Kattman et al., 2006; Motoike et al., 2003; Yamashita et al., 2000). However, high plasticity, heterogeneity and rapid differentiation of early embryonic cells have hampered investigations on the lineage relationships among these cell types, the mechanisms underlying the lineage restriction steps regulating Flk1⁺mesoderm development, and the following commitment from Flk1⁺ mesoderm to its downstream lineages. Recent progress in single cell RNA profiling technology has advanced our understanding of the complex cell constitution at discrete early embryonic stages (Farrell et al., 2018; Gong et al., 2017; Ibarra-Soria et al., 2018; Lescroart et al., 2018; Mohammed et al., 2017; Moignard et al., 2015; Pijuan-Sala et al., 2019; Scialdone et al., 2016; Wagner et al., 2018) or in in vitro-differentiated embryonic stem (ES) cells (Chu et al., 2016; Han et al., 2018; Spangler et al., 2018). However, none of these studies has provided a complete process from pluripotent cells to Flk1⁺ mesoderm owing to missing intermediate states or batch effects. Therefore, there is still a gap in our current knowledge on the continuous sequence of Flk1⁺ mesoderm development and its segregation into the hemangiogenic lineage.

Mouse ES cells in culture differentiate in an asynchronous manner giving rise to embryoid bodies (EBs). Both undifferentiated pluripotent cells and differentiated progenies (Park et al., 2013; Spangler et al., 2018) co-exist in mouse EBs at early stages of differentiation. We reasoned that cells in the intermediate states also exist in the same culture, representing the full spectrum of differentiation, which would provide a unique opportunity to delineate a continuous differentiation process without batch effects. By taking advantage of the asynchronous nature of the differentiation of mouse ES cells and culture conditions that favor hemangiogenic lineage differentiation, we performed single cell RNA sequencing (scRNA-seq) of EB cells to elucidate a continuous developmental path from the pluripotent state to Flk1⁺ mesoderm, and to the hemangiogenic lineage. We found that the cardiac lineage branched out immediately after nascent mesoderm specification. Intermediate states positioned between ES and Flk1⁺ mesoderm were distinguished by a few temporally overlapping transcription factor modules, suggesting a plastic and combinational regulatory mechanism regulating Flk1⁺ mesoderm development. Notably, expression of most of the characteristic genes of Flk1⁺ mesoderm continuously increased in the smooth muscle branch, implying a ‘default’ route. Our data suggested that ETV2 turns on the hemangiogenic program from Flk1⁺ mesoderm by restraining the smooth muscle fate. We further identified that the SRC kinase potentially regulates ETV2-mediated hemangiogenic lineage specification, possibly through mediating cell adhesion signaling. These studies provide a comprehensive description of mesoderm and hematopoietic/endothelial cell lineage development at the single cell level.

In capturing the dynamics of transcriptome of the hemangiogenic differentiation, we performed scRNA-seq using day 4 EB cells in hemangiogenic differentiation. We reasoned that day 4 would allow access to a range of cells from pluripotent cells to those of the hemangiogenic cell lineage, as hemangiogenic cells extensively emerge and pluripotent stem cells are still present at this time. After filtering out low-quality cells, 1848 cells, which had 4373 genes/27,272 unique mRNAs detected per cell in average (Fig. 1A), were kept for further analysis. Day 4 EB cells were grouped into 18 populations (Fig. 1B, Fig. S1A) and 936 genes were differentially upregulated across the populations (Table S1). t-Stochastic neighbor embedding (t-SNE) analysis indicated that most clusters did not form a distinct, well-segregated population (Fig. 1B). Instead, distinct populations branched out gradually from the trunk, underpinning continuous developmental processes. We assigned different fates to the branches according to the expression of related marker genes (Table S1). SPRING (Weinreb et al., 2018), a different tool for visualization of high-dimensional scRNA-seq data, produced a similar clustering pattern of the cells (Fig. S1B). High expression of Sox2, Pou5f1 (also known as Oct4), Nanog and Zfp42 (also known as Rex1) highlights the maintenance of undifferentiated naïve pluripotent cells in the EBs (Fig. 1C,D). Flk1⁺ mesoderm and the hemangiogenic lineage, as marked by high expression of Flk1 and Etv2, and upregulation of the ETV2 target genes Tal1 and Gata2, clustered at the other end of the ES cell population (Fig. 1C,D). Expression of related genes showed a gradual shift in intermediate clusters from the pluripotent ES population to the Flk1⁺ mesoderm branch, revealing continuous intermediate states from ES to Flk1⁺ mesoderm that are accessible in EBs.

To determine in further detail the developmental relationships among the clusters between the ES and Flk1⁺ populations, we first examined the dynamics of the loss of pluripotency and the bifurcation of mesoderm and endoderm fates. Pluripotency can be maintained in two states: ‘naïve’ and ‘primed’ (Kalkan and Smith, 2014). The former represents the primitive phase and is characterized by high expression of Zfp42, whereas the latter loses Zfp42 expression and harbors an unstable pluripotency gene network, which marks a nascent stage of gastrulation/mesendoderm (Kalkan and Smith, 2014). We picked out clusters that were located in the endoderm path (4, 5, 8, 9, 11 and 17 in Fig. 1B), and re-clustered the cells based on the 163 most differentially expressed genes among these cells (Fig. 2A). A narrow ‘path’ from the Utf1-expressing pluripotent population leads to the Sox17-expressing endoderm (the region surrounded by a red dashed line), next to the sprouting Mesp1-expressing nascent mesoderm. We next ordered cells in the endoderm path into a pseudo-time developmental line and examined the kinetics of the 163 most differentially expressed genes (Fig. 2B-D, Fig. S2A). At the beginning of the path, we observed an exit from the naïve pluripotency (Fig. 2C, flanking the blue dashed line). Unexpectedly, expression of a small group of genes (represented by Dppa3; Fig. 2C and red box in Fig. 2B) already displayed strong heterogeneity in the Zfp42-expressing population. Gradual loss of the expression of these genes preceded the exit from naïve pluripotency. This may suggest an extra layer of guarding mechanisms of pluripotency beyond the naïve circuit (as represented by the genes in the blue box in Fig. 2B, which have relatively constant expression in the naïve pluripotency stage). Primitive streak markers such as T and Fgf5 were upregulated in cells exiting from the naïve pluripotent state and entering the mesendoderm stage (Fig. 2A). Shortly after, the endoderm and mesoderm lineages bifurcated. Fgf5 expression became more enriched in endoderm, whereas T expression was exclusively enriched in mesoderm (Fig. 2A). Gata6 and Gsc were expressed in both mesoderm and endoderm, although they were expressed more extensively in the latter (Fig. 2A). The majority of cells in the endoderm population were separated from the trunk, with very few cells expressing Foxa2 at a high level in the intermediate state. When we lowered the stringency and included more variable genes to cluster the cells, the connection between endoderm and the mesendoderm cluster became more obvious (Fig. S2B,C). These results suggest that the core endoderm genes form a strong self-enforcing circuit, rendering the intermediate specification state relatively unstable.

The immediate upstream population of Flk1⁺ mesoderm expressed Pdgfra and Mesp1 (Fig. 2E, region marked by a blue dashed line). These cells also expressed low levels of Tbx6, implying that, to some extent, they are ‘primed’ for paraxial mesoderm, which generates tissues including skeletal muscle and cartilage. As such, they might represent an early ancestral state for both paraxial mesoderm and Flk1⁺ lateral plate mesoderm. Consistently, cells in this region had an expression pattern similar to embryonic day (E) 7.0 posterior primitive streak, where nascent mesoderm emerges (Fig. 2F). Furthermore, these cells were transcriptionally similar to Mesp1-labeled early mesoderm cell population in mouse E6.75 embryos in recent scRNA-seq data from Lescroart et al., (2018). Similar to EB-derived progenitors, this embryonic population of cells still expressed high levels of the pluripotency gene Pou5f1, and already expressed high levels of Mesp1 and Pdgfra. Moreover, although they expressed low levels of Tbx6, they have not yet upregulated Flk1 (cluster 1 in Fig. S2D). The cardiac lineage branched out mainly from the Pdgfra⁺Flk1⁻ early mesoderm population, and a small group of Flk1⁺ mesoderm also expressed the cardiac transcription factor Isl1 (Fig. 2G). Notably, smooth muscle and hemangiogenic lineages directly bifurcated from an Flk1-expressing population, suggesting a common origination (Fig. 2H). RNA velocity estimation of the variable genes, based on the abundance of unspliced and spliced mRNAs (La Manno et al., 2018), revealed persistent overall cell state dynamics from naïve pluripotency to differentiation branches (Fig. S2E). These results suggested a continuous process from pluripotency to Flk1⁺ mesoderm and its derivate tissues (Fig. 2I).

To investigate whether the finding that pluripotent and early gastrulation populations co-exist in day 4 EBs is due to differential cell proliferation, we examined expression of cell cycle-related genes. The expression pattern of typical S-phase marker genes showed no obvious bias in different cell populations (Fig. 3A, Fig. S3A). Using the cell cycle analysis tool in the Seurat package (Butler et al., 2018), we assigned each cell to a specific phase in the cell cycle and found that the majority of cells in the mesoderm developmental route proliferated robustly (Fig. 3B, Fig. S3B). Analysis of cell proliferation using the violet proliferation dye 450 (VPD450) also confirmed that EB cells with different levels of differentiation (indicated by the reporter of the hemangiogenic marker gene Etv2) proliferated similarly in our system (Fig. 3C, top). As a control, differentiation block and cell cycle delay by the ERK inhibitor PD0325901 was observed in the VPD450 dilution assay (Fig. 3C, bottom).

Next, we retrieved a continuous process of the hemangiogenic lineage development by manually choosing the populations corresponding to the shortest path from the naïve pluripotent state to the hemangiogenic branch and reconstructed a fine developmental route (Fig. 4A). A group of Flk1-expressing cells achieved the highest level of Etv2 expression, followed by dramatic upregulation of the ETV2 target genes Tal1 (also known as Scl) and Gata2, confirming that a threshold level of Etv2 expression is necessary for the hemangiogenic lineage specification (Zhao and Choi, 2017) (Fig. 4A). We ordered the cells into a one-dimensional pseudo-time line and assigned the populations to specific stages according to the dynamics of marker gene expression (Fig. 4B). The 68 most variable transcription regulation-related genes were clustered based on their expression along the pseudo-time line (Fig. 4C). They aggregated into a few obvious modules covering different developmental stages (blue dashed line boxed regions in the heatmap). Notably, the naïve pluripotency module (Fig. 4C, box II) and the hemangiogenic module (box III) were relatively independent, whereas each of the intermediate gastrulation modules were shared by adjacent stages (Fig. 4C boxes I, IV and V). The gene network constructed based on the co-expressed genes also highlighted the collaboration among transcription factors of individual modules (Fig. 4D). This stepwise combinatorial usage of different regulatory circuits fits the diverse lineage specifications in gastrulation. Consistent with this, the intermediate populations were less committed, largely interchangeable (Sakurai et al., 2006), and more similar to each other in transcriptome (Fig. 4E, Fig. S4A,B). These results suggested that the intermediate states constitute a ‘relay’-like multistep specification process (Fig. 4F). Exit from the naïve pluripotency and activation of the hemangiogenic program might be the two most crucial transitions in hemangiogenic lineage development.

Owing to intracellular noises and microenvironmental influences, gene expression within individual cells can fluctuate. scRNA-seq analysis offers an opportunity to utilize this luxuriant gene expression heterogeneity to assess a regulatory network. We previously reported that Etv2 expression above a threshold level is essential to activate hemangiogenic genes (Zhao and Choi, 2017). We confirmed this at a single cell endogenous mRNA level (Fig. 5A). We noticed that the kinetics of upregulation of ETV2 target genes soon became extremely steep after Etv2 expression reached a specific level, suggesting an Etv2-mediated sensitive ‘switch’ mechanism in activating its downstream target genes. We explored potential ETV2 cooperating genes for triggering this sensitive switch by examining potential upstream regulators of Tal1 and Lmo2, two ETV2 direct target genes that are crucial for establishing the hemangiogenic lineage (Liu et al., 2015; Porcher et al., 1996). We assumed that if more gene A-expressing cells simultaneously expressed gene B, it is likely that gene A expression depends on gene B (Fig. 5B; see Materials and Methods for more details). From the population in which ETV2 starts to activate hemangiogenic genes (population 12 in Fig. 1B), we assessed the possibility of a given gene being required for expression of Tal1 or Lmo2 (Fig. 5C,D). Etv2 and Flk1 were at the top of the list of genes required for Tal1 or Lmo2 expression (Fig. 5C,D, Table S2), consistent with the finding that VEGF-FLK1 signaling ensures Etv2 threshold expression to activate hemangiogenic genes (Zhao and Choi, 2017). The comparison of the top 100 genes required for Tal1 expression and the top 100 required for Lmo2 expression generated 66 in common, of which 22 were annotated to cell adhesion-related signaling (Fig. S5A, highlighted in Table S2). The majority of these 22 genes were already highly expressed in populations preceding the hemangiogenic stage (Fig. S5B). To test if the cell adhesion signaling plays a role in ETV2-mediated hemangiogenic gene activation, we utilized a doxycycline (DOX)-inducible Etv2-expressing ES cell line (Liu et al., 2012). Because extracellular matrix/cell adhesion signaling converges on the SRC kinase (Cabodi et al., 2010), we added the SRC inhibitor PP2 to day 3 EBs, while simultaneously inducing exogenous Etv2 expression. After 24 h, as expected, DOX-induced Etv2 expression greatly enhanced the FLK1⁺PDGFRα⁻ population that is enriched in the hemangiogenic lineage cell population (Liu et al., 2012), whereas inhibition of SRC using PP2 was sufficient to reduce Etv2 overexpression-induced hemangiogenic lineage skewing (Fig. S5C). Given that cell adhesion is necessary for gastrulation as well (Adams and Watt, 1993), it was possible that PP2 affected earlier stages rather than directly affecting the function of ETV2. To minimize this possibility, we induced ETV2-mCherry fusion protein expression ectopically in undifferentiated ES cells, in which the endogenous Etv2 level is negligible (Zhao and Choi, 2017) and exogenous Etv2 can induce Flk1 expression. Etv2 is also downstream of FLK1 signaling and the reciprocal activation between ETV2 and FLK1 signaling is an important mechanism for hemangiogenic fate determination (Zhao and Choi, 2017). Exogenous ETV2 induced more than 8% of FLK1⁺ cells in undifferentiated ES cells, and this ratio was significantly reduced by PP2 treatment (Fig. 5E). Unexpectedly, PP2 treatment also dramatically reduced ETV2 protein level, with only a limited number of cells achieving threshold expression of Etv2 (Fig. 5E). To exclude the possibility that PP2-mediated SRC inhibition affected DOX-induction system or general gene expression processes, we tested DOX-induced expression of another transcription factor, EOMES, and the constitutively active promoter EF1α-driven expression of the red fluorescent protein mCherry. We found that DOX-induced EOMES protein level only mildly decreased in the presence of PP2, while EF1α-driven mCherry expression was insensitive to this treatment (Fig. S5D). These results suggest that SRC signaling contributes to ETV2-mediated hemangiogenic fate specification possibly through upstream cell adhesion signals.

Based on the observation of the direct bifurcation of the smooth muscle and hemangiogenic lineages from Flk1⁺ mesoderm (Fig. 2I), we analyzed the cell populations corresponding to the developmental route of the smooth muscle lineage (Fig. 6A, Fig. S6A), and examined expression of the 427 most variable genes along the pseudo-time of differentiation (Fig. 6B). Many of the variable genes aggregated as modules. The modules corresponding to mesendoderm, nascent mesoderm, and Flk1-expressing mesoderm stages displayed largely overlapping ‘relay’-like patterns. Three groups of genes were most extensively expressed in the smooth muscle branch (Fig. 6B, gene clusters a, b and c in the black dashed line boxed region). Of the three groups of genes, groups b and c were already upregulated prominently in Flk1⁺ mesoderm, with cluster c being upregulated even earlier from the mesendoderm stage. The transcription factors Hand1, Tbx20 and Foxf1 were found in gene cluster b (Fig. 6B,C). Hand1 and Foxf1 have been reported to be important for smooth muscle development (Mahlapuu et al., 2001; Morikawa and Cserjesi, 2004). Expression of the genes in clusters b and c remained similar in the hemangiogenic branch. Only a small cluster of genes was exclusively upregulated after cells entered the smooth muscle branch (gene cluster a). Importantly, all the smooth muscle transcription factor genes, except Tbx2, were already extensively expressed in the Flk1⁺ mesoderm and constituted a major characteristic of the Flk1⁺ mesoderm (Fig. S6B, in the black dashed box). These results suggest that smooth muscle might be the default fate of the Flk1⁺ mesoderm.

We investigated further populations 0, 3, 10, 12 and 13 (as labeled in Fig. 1B), corresponding to the bifurcation of hemangiogenic and smooth muscle fates (Fig. 6D). Expression of Bmp4, which can enhance vascular smooth muscle (VSM) development (Cheung et al., 2012), was dramatically upregulated in the smooth muscle branch [from population 0 (black) to population 13 (yellow)], but not in population 12 (green), where Etv2 initiates specification of the hemangiogenic fate (Fig. 6D). Hand1 and Foxf1 were dramatically downregulated in population 12 compared with the upstream population 3, implying an active repression of them by Etv2 (Fig. 6E). Consistent with this, by analyzing the work of Wareing and colleagues (Wareing et al., 2012) we found that re-expression of Etv2 in Etv2-knockout EBs inhibited expression of Hand1 and Foxf1 (Fig. 6F). These results suggest that Etv2 represses the default smooth muscle program when specifying the hemangiogenic lineage from Flk1⁺ mesoderm.

To understand better the relationship between Flk1⁺ mesoderm and the smooth muscle cell program, we further analyzed the overall expression score of the group b genes in Fig. 6B (listed in Table S3). The group b gene score reflected the dynamics of Flk1⁺ mesoderm identity in day 4 EB data (Fig. S7A). The cardiac lineage progenitors also had a high group b gene score; as such, we termed it the lateral plate mesoderm score, ‘LPM-score’ (both Flk1⁺ mesoderm and Flk1⁻ cardiac tissues arise from lateral plate mesoderm). If the Flk1⁺ mesoderm's default fate is indeed smooth muscle lineage in vivo, it was expected that Flk1⁺ mesoderm identity would persist in the smooth muscle lineage. To determine whether this is the case, we performed scRNA-seq of a mixed cell pool of E9.5 and E10.5 mouse yolk sac, an extra-embryonic tissue that generates the first-arising hemangiogenic lineages and vasculature structure (Fig. 7A). We chose E9.5 and E10.5 yolk sac to ensure sufficient representation of fully committed smooth muscle lineage, because published single cell analysis of early mouse embryos from implantation to E8.5 (Ibarra-Soria et al., 2018; Pijuan-Sala et al., 2019) did not have good coverage of extra-embryonic smooth muscle tissues. We found that the LPM-score was particularly high in smooth muscle cells of the mouse yolk sac [VSM; Fig. 7B], exemplified by the co-expression of the smooth muscle actin gene Acta2 and the transcription factors Foxf1 and Hand1, and the absence of the expression of the hematopoietic and endothelial transcription factor Tal1 (Fig. 7C, Fig. S7B). It has been reported that the extra-embryonic amnion expresses smooth muscle genes (Fleury et al., 2015). In recently reported scRNA-seq of E8.25 mouse embryos (Ibarra-Soria et al., 2018), the amnion population had a high LPM-score and expressed smooth muscle markers (Fig. 7D-F, Fig. S7C,D). Consistent with this, the 72 marker genes shared by day 4 EB smooth muscle, yolk sac VSM, and amnion were enriched with smooth muscle function annotation (Fig. 7G,H). Collectively, smooth muscle lineage might be the default fate of Flk1⁺ mesoderm.

It has been reported that single Flk1⁺ mesoderm cells can form both VSM and endothelial cells in culture and in vivo (Yamashita et al., 2000). Moreover, sorted FLK1⁺ cells are able to form cardiomyocytes and even skeletal muscle (Liu et al., 2012; Sakurai et al., 2006). However, the lineage relationships among cardiac, VSM, skeletal muscle and hemangiogenic tissues concerning Flk1⁺ mesoderm allocation remain elusive. Our single cell transcriptome profiling of EB cells revealed a continuous differentiation trajectory. We showed that the cardiac branch mainly segregates out from Pdgfra⁺Flk1⁻ nascent mesoderm. Unexpectedly, smooth muscle and hemangiogenic lineages have the closest lineage relationship and directly bifurcate from a common group of Flk1-expressing cells. Notably, most Flk1⁺ mesoderm stage markers were further upregulated in the smooth muscle branch, suggesting that smooth muscle might be the default fate of the Flk1-expressing mesoderm. We propose that ETV2 drives the branching of the hemangiogenic lineage from the smooth muscle fate. Expression of cell adhesion signaling coincides with ETV2 target gene activation and hemangiogenic lineage specification. Consistent with this, inhibition of SRC, a kinase important for cell adhesion signaling, was sufficient to repress hemangiogenic lineage skewing induced by Etv2 overexpression. Meanwhile, crucial smooth muscle transcription factors, Hand1 and Foxf1, become actively downregulated. Notably, neither Hand1 nor Foxf1 are direct ETV2 target genes (Liu et al., 2015). It will be important in the future to explore how this repression is achieved.

Hemangiogenic lineage development can be clearly annotated by a series of transcription factor modules, therefore providing strong molecular support for the current definition of early embryo developmental stages (Fig. 4C). Cells in the gastrulation stages, after exiting pluripotency and before entering the hemangiogenic state, are overall similar to each other in transcriptome, and the underlying transcription factor modules show a ‘relay’-like highly overlapping pattern. Consistent with this, FLK1⁻PDGFRα⁻, FLK1⁻PDGFRα⁺ and FLK1⁺PDGFRα⁺ populations in early differentiation of ES cells are reversible/interchangeable (Sakurai et al., 2006). These results suggest that exit from the naïve pluripotent state and activation of the hemangiogenic program might be the two rate-limiting steps in hemangiogenic lineage development, and that the intermediate gastrulation stages are plastic and adaptable to multiple lineage fates through combinational usage of limited regulatory circuits (Fig. 7I).

In summary, our work describes a continuous process of a mesoderm lineage tree leading to hemangiogenic and smooth muscle fates and the underlying molecular networks, thereby providing comprehensive information on early embryo and blood/blood vessel development. These findings are fundamental for the study of basic cell fate determination and gene network structures, and for designing more effective strategies to generate hematopoietic and endothelial cells for regenerative medicine.

Mouse ES cells were maintained and differentiated as previously reported (Liu et al., 2012; Zhao and Choi, 2017). Briefly, ES cells were maintained on mouse embryo fibroblast (MEF) feeder cells in Dulbecco-modified Eagle medium containing 15% fetal bovine serum, 100 units/ml LIF, 1× MEM Non-Essential Amino Acids Solution (Gibco), 1× Glutamax Supplement (Gibco) and 4.5×10⁻⁴ M 1-thioglycerol (MTG, Sigma-Aldrich). For feeder-free culture, ES cells were maintained in a gelatin-coated dish in Iscove's modified Dulbecco medium (IMDM) with the same supplements used for maintenance. For EB differentiation, ES cells were first cultured in a feeder-free condition for 2 days. After culture, single cell suspensions were prepared and 8000 cells were added per ml of differentiation medium containing IMDM, 15% pre-screened fetal calf serum, 1× Glutamax, 50 µg/ml ascorbic acid and 4.5×10^–4 M MTG. Differentiation was carried out on a bacterial Petri dish. DOX-inducible ETV2 (Liu et al., 2012) or ETV2-mCherry (Zhao and Choi, 2017) ES cells have been previously described in detail. Briefly, coding sequences of ETV2 or ETV2-mCherry were targeted into the Rosa26 locus by homologous recombination using A2Lox ES cells carrying the TRE-loxP-Δneo inducible target locus (Iacovino et al., 2009). The SRC inhibitor PP2 (Sigma-Aldrich) and DOX (Sigma-Aldrich) treatment were performed as indicated in the text.

Yolk sacs were dissected from equal numbers of mouse E9.5 and E10.5 embryos, mixed together and dissociated with 0.25% collagenase at 37°C for 30 min. The cells were stored at −80°C in fetal bovine serum with 10% DMSO. When the 10x Genomics instrument was ready, the cells were thawed and washed, and stained using antibody against Ter119 (Ly76) (1:200, 12-5921-82, eBioscience). Dead cells and Ter119⁺ cells were excluded by sorting to enrich live non-erythroid cells.

Day 4 EBs were dissociated with Accutase solution (Sigma-Aldrich) for 5 min into single cells. The single cell suspension was washed and re-suspended in PBS.

Single cell suspension at 300 cells/ µl in PBS were subjected to Chromium 10x Genomics library construction and HiSeq2500 sequencing (The Genome Technology Access Center, Washington University in St. Louis). The sequenced reads were mapped to the GRCm38 assembly using Cell Range 2.0.1 (10x Genomics).

The scRNA-seq output was imported into Seurat 1.4 (Butler et al., 2018), and genes expressed in at least three cells were kept for analysis. Cells with more than 5% mitochondria reads or fewer than 2000 unique genes detected were filtered out; 1848 cells, which had 4373 genes/27,272 unique mRNAs detected per cell in average, were kept for further analysis. Gene expression counts were normalized and log transformed.

According to the average expression and dispersion gene expression of each gene, we selected the 427 most variable genes for further analysis of the remaining cells. nUMIs (number of detected UMIs in each cell) were re-scaled to regress the effects of detection depth of each cell and mitochondria reads. The re-scaled data were run for principal components analysis (PCA) to reduce dimensions. According to the PCA results, the cells were clustered into 18 populations, and also used for t-SNE presentation. Cells in selected populations were imported into Monocle 2 (Trapnell et al., 2014) for re-clustering, PCA plotting and pseudo-time ordering.

For hierarchical clustering and heatmap plotting of selected genes or cells, the R package ‘pheatmap’ was used. All scRNA-seq data analyses were finished in R environment.

For cell cycle score and phase assignment with Seurat, the ‘CellCycleScoring’ function was invoked.

For LPM-score gene analysis, group b genes in Fig. 6B were processed by invoking the ‘AddModuleScore’ function in Seurat.

For SPRING visualization, a matrix of expression counts of the 427 most variable genes versus the filtered cells was uploaded (https://kleintools.hms.harvard.edu/tools/spring.html) per the guidance (Weinreb et al., 2018). Cell cluster information from Seurat analysis was also loaded for viewing.

For RNA velocity estimation, the bam file in the 10x Genomics output data was first re-counted with the velocyto counting pipeline in Python according to the manual (La Manno et al., 2018). The generated loom file was loaded to velocyto.R. Only the counts of the 427 most variable genes were used for RNA velocity estimation. Finally, coordinates from Seurat's t-SNE analysis were used to embed the velocity results.

For gene regulatory network analysis, the R package GENIE3 (Huynh-Thu et al., 2010) was used. To simplify the results, only transcription-related genes in the 427 most variable genes were used to calculate the gene-gene link matrix, and only the top 250 weighted link were adopted for visualization. To view the network, the R package igraph was used.

For Zipcode mapping, the mean expression of all the genes in chosen cells was uploaded to www.itranscriptome.org/ according to the developer's guidance (Peng et al., 2016).

Cells within cluster 12 were divided into two groups, one with Tal1 expression above the mean Tal1 level across the whole population, and the other with Tal1 expression below the mean level. Cells in the first group (with gene x expression above its mean level) were counted separately from cells in the second group (with gene x expression below its mean level), and the sum of the two counts was divided by the total cell number of cluster 12. This ratio was considered as the possibility of gene x following Tal1. In fact, close gene-gene correlation can be attributed to three possibilities (with Tal1 used as an example): (1) gene x lies upstream of Tal1 and is required for its expression, (2) gene x is TAL1's target, (3) both gene x and Tal1 are ETV2's target and thus respond to ETV2 in a similar way. However, in the population in which expression of Tal1 has just been initiated, it is not likely that TAL1 activates its own target genes yet, reducing the likelihood of possibility 2. We can refer to a gene's expression pattern in the differentiation trajectory to test possibility 3. In fact, most of the genes tightly following Tal1 and Lmo2 were already extensively activated before ETV2 started to function.

Etv2-tdTomato ES cells (Zhao and Choi, 2017) were dissociated into single cell suspension and incubated with 1 µM of VPD450 (BD Horizon) in PBS for 10 min at room temperature. Cells were then aggregated for normal differentiation.

Wareing and colleagues' microarray data (Wareing et al., 2012) were re-analyzed using the Expression Console software (Affymetrix).

Gene functional annotation was performed using Metascape (http://metascape.org/gp/index.html) (Tripathi et al., 2015).

Single cells were incubated with anti-mouse FLK1 (1:400, 136406, BioLegend) and anti-mouse PDGFRα (1:500, 135908, BioLegend). Data were acquired on LSR-Fortessa flow cytometer (BD Biosciences) and analyzed using the FlowJo (Treestar) software.

The results of flow cytometry were analyzed using Student's t-test. P<0.05 was considered significant.

We thank members of the Choi lab for helpful discussion and support.