Transcriptome Dynamics of Hematopoietic Stem Cell Formation Revealed Using a Combinatorial Runx1 and Ly6a Reporter System.

Studies of hematopoietic stem cell (HSC) development from pre-HSC-producing hemogenic endothelial cells (HECs) are hampered by the rarity of these cells and the presence of other cell types with overlapping marker expression profiles. We generated a Tg(Runx1-mKO2; Ly6a-GFP) dual reporter mouse to visualize hematopoietic commitment and study pre-HSC emergence and maturation. Runx1-mKO2 marked all intra-arterial HECs and hematopoietic cluster cells (HCCs), including pre-HSCs, myeloid- and lymphoid progenitors, and HSCs themselves. However, HSC and lymphoid potential were almost exclusively found in reporter double-positive (DP) cells. Robust HSC activity was first detected in DP cells of the placenta, reflecting the importance of this niche for (pre-)HSC maturation and expansion before the fetal liver stage. A time course analysis by single-cell RNA sequencing revealed that as pre-HSCs mature into fetal liver stage HSCs, they show signs of interferon exposure, exhibit signatures of multi-lineage differentiation gene expression, and develop a prolonged cell cycle reminiscent of quiescent adult HSCs.

Introduction

Developmental hematopoiesis occurs in three phases: first, primitive erythroid cells and macrophages develop in the yolk sac (YS) starting around embryonic day 7.5 (E7.5) (Ferkowicz et al., 2003, Palis et al., 1999). By E8.5, the YS has begun to produce definitive erythroid and myeloid progenitors (McGrath et al., 2015, Palis et al., 1999) as well as the first lymphoid progenitors (LPs) (Godin et al., 1995, Nishikawa et al., 1998). The third wave produces definitive hematopoietic stem cells (HSCs) that are able to engraft, self-renew, and produce all hematopoietic lineages in myeloablated adult recipients (Orkin and Zon, 2008). HSCs appear around E10–E11 in the aorta-gonad-mesonephros (AGM) region (Medvinsky and Dzierzak, 1996, Muller et al., 1994), the umbilical and vitelline arteries (AUV) (de Bruijn et al., 2000), YS (Huang and Auerbach, 1993, Yoder et al., 1997), and the placenta (PL) (Gekas et al., 2005, Ottersbach and Dzierzak, 2005) (reviewed in Dzierzak and Bigas, 2018).

During definitive hematopoiesis, hemogenic endothelial cells (HECs) turn into hematopoietic cluster cells (HCCs) in a process called endothelial-to-hematopoietic transition (EHT) (Bertrand et al., 2010, Boisset et al., 2010, Jaffredo et al., 1998, Kissa and Herbomel, 2010). By E10.5, the hematopoietic clusters comprise hundreds of Kit+ cells in the dorsal aorta (DA) and the AUV (Yokomizo and Dzierzak, 2010); however, the vast majority of these cells still lack HSC activity (Muller et al., 1994). HSC-fated cells emerge from HECs as pre-HSCs that subsequently mature into HSCs (Ganuza et al., 2017, Henninger et al., 2017, Rybtsov et al., 2016, Taoudi et al., 2005).

Pro-HSCs (VE-cadherin [VEC]+CD41loCD43−CD45−), detectable as early as E9.5, progress through pre-HSC I (VEC+CD41loCD43+CD45−) and pre-HSC II (VEC+CD41loCD43+CD45+) stages before they acquire HSC potential (Rybtsov et al., 2011, Rybtsov et al., 2014, Taoudi et al., 2008). However, CD41, CD43, CD45, and VEC are broadly expressed by all cells that undergo EHT, irrespective of their potential. Alternatively, definitive HECs and hematopoietic stem/progenitor cells (HSPCs) can be marked using reporter genes that utilize the +23 enhancer of Runx1 (Bee et al., 2010, Ng et al., 2010, Nottingham et al., 2007), a transcription factor gene critical for definitive hematopoietic development (Chen et al., 2009). While the Runx1+23 enhancer is active in all emerging definitive hemogenic and hematopoietic cells, the Ly6a (Sca-1) gene promoter specifically marks pre-HSC-producing HECs (Chen et al., 2011, de Bruijn et al., 2002). Ly6a-GFP is not expressed in YS blood islands during the first wave of hematopoiesis (Chen et al., 2011). However, Ly6a-GFP+ cells are also present in non-hemogenic tissues (de Bruijn et al., 2002), and the reporter continues to be expressed in many lineage-committed blood cells (Ma et al., 2002). Together, these findings show that no single marker suffices to track HSC specification, highlighting the importance to improve currently available tools.

We now report that by combining a Runx1+23 enhancer regulated mKO2 reporter (Runx1-mKO2) with a Ly6a-GFP reporter we were able to accurately mark HECs and HCCs and follow their maturation into (pre-)HSCs and hematopoietic progenitor cells (HPCs). Erythro-myeloid HPCs are found in the Runx1-mKO2+ compartment (irrespective of Ly6a-GFP activity), whereas LPs and functional HSCs are restricted to the reporter double-positive (DP) compartment. We found HECs capable of producing DP pre-HSC-like cells in both the YS and E9.5 para-aortic splanchnopleura (PSp)/E10.5 AUV. However, robust HSC activity emerged later (E11.5), and most prominently in the PL. Using single-cell analyses of pre-HSC I, pre-HSC II/HSC, and fetal liver (FL) HSC transcriptomes we identified transcription factors, receptors, and processes whose expression correlates with this HSC development, including downregulation of cell-cycle genes, upregulation of interferon-induced genes, and regulators of multi-lineage differentiation. Thus, our data suggest that interferon exposure plays a critical role in pre-HSC maturation and that cycling FL HSCs are already primed to enter the quiescent state typical of adult long-term (LT) HSCs.

Results

The Runx1-mKO2 and Ly6a-GFP Dual Reporter System Specifically Marks HECs, HCCs, and HSPC during Definitive Hematopoiesis

We developed a new reporter construct (Figure S1A) in which the Runx1+23 enhancer drives expression of an mKO2 reporter fused to H2B to stabilize and enrich the signal in the nucleus. Two independent transgenic mouse lines with similar expression patterns and intensity were identified for further analysis (Figures S1B–S1D). Since HSCs emerge from a subpopulation of endothelial cells in which the Ly6a promoter is active (Chen et al., 2011, de Bruijn et al., 2002), Runx1-mKO2 mice were bred with Ly6a-GFP mice (Ma et al., 2002) to create Runx1-mKO2 and Ly6a-GFP dual reporter mice.

Consistent with endogenous Runx1 expression (Tober et al., 2013, Yzaguirre and Speck, 2016), mKO2 fluorescence was observed at E8.5 in Kdr-GF+ YS blood islands (Figure S1E); at this time, Ly6a-GFP expression is still absent (data not shown). By E9.5, mKO2+ HCCs have developed in the umbilical artery (UA) and vitelline artery (VA) (see Figures S1F and S1G). The HCCs of the E9.5/E10.5 VA were often large enough to allow detection by stereo fluorescence microscopy (Figures S1F and S1H). GFP was expressed by a fraction of the E9.5/E10.5 endothelial cells of the UA, VA, and YS (Figures S1F–S1H) and by some non-hemato-endothelial cells in the tail region (Figures S1F and S1H).

To study the activity of the reporter genes in the embryo proper and in more detail, we performed multi-color high-resolution 3D confocal microscopy of E9.5/E10.5 embryos stained for CD31 (marking endothelial cells and HCCs), or Kit (marking HCCs and HSPCs) (Figures 1A–1E, S2A, and S2B). We found that mKO2+ cells in the endothelium first form HCCs in the VA in “early” E9.5 (22 somite pairs [sp]) embryos (Figure 1A). Rare mKO2/GFP DP “bulging” endothelial cells (Yzaguirre and Speck, 2016) were present (Figure 1Aʹ); however, most of the HCCs remained Runx1-mKO2 single positive (R-SP) (Figure 1A″). By “mid” E9.5 (25 sp), the number and cellularity of these mKO2+Kit+ HCCs has increased rapidly, particularly near the base of the VA where mesenteric blood islands emerged (Figure 1B) (Garcia-Porrero et al., 1995). Some of the cells at the base of these clusters were DP for mKO2 and GFP (Figure 1Bʹ). HCC formation activity was also apparent in the UA. At this stage, the FL only contained a few R-SP, and virtually no DP cells (Figure 1B). By “late” E9.5 (28 sp), the number of mesenteric blood islands has decreased, small HCCs have begun to appear in the DA, and the FL has become colonized by mostly R-SP-positive cells (Figure 1C). Around E10.5, when HCC activity is at its peak (Yokomizo and Dzierzak, 2010), we observed sizable intra-aortic clusters in the AGM. Virtually all DA, UV, and VA HCCs showed strong Runx1-mKO2, CD31, and Kit staining (Figures 1D, 1E, S2A, and S2B). However, not all Runx1-mKO2+ cells co-expressed GFP, and the GFP expression levels as well as the frequency of GFP-expressing cells varied within and between individual clusters.

Figure 1

Analysis of Runx1-mKO2 and Ly6a-GFP Reporter Expression during the Emergence of Hematopoietic Cluster Cells

Whole-mount confocal Z-projection of E9.5–E10.5 embryos showing the merged image of expression of Runx1-mKO2 (red), Ly6a-GFP (green), and CD31 (A, C, and E) or Kit (B and D) (blue). Arrows point to HCCs. Arrowheads point to mKO2+ cells. Green arrowheads point to GFP+ cells. (A) Low magnification Z-projection image of E9.5 (22 somite pairs [sp]) showing dorsal aorta (da), umbilical artery (ua), and vitelline artery (va) region (40 sections/z = 2.78 μm). Boxed areas are magnified in A′ and A″. (A′) Bulging mKO2/GFP DP endothelial cells (ECs). (A″) Runx1-mKO2 single positive (R-SP) ECs and HCCs, and a flat DP EC (20 sections/z = 0.54 μm) (B) Low magnification Z-projection image of E9.5 (25 sp) embryo (63 sections/z = 2.53 μm). Arrows point to the prominent mKO2+ and Kit+ HCCs form along the vitelline artery delta fusing with dorsal aorta. Runx1-mKO2+ cells (asterisk) start to colonize in fetal liver (fl). Boxed area is magnified in Bʹ. (B′) Large Kit+ HCCs emerging into the lumen from the underlying endothelium (20 sections/z = 0.54 μm). (C) Low magnification Z-projection image of E9.5 (28sp) embryo (80 sections/z = 2.53 μm). Significant increase in Runx1-mKO2+ cells (asterisk) in fl. Boxed area is magnified in C’ (20 sections/z = 0.54 μm). (D and E) High magnification Z-projection image of HCCs in the da of E10.5 (36 sp) embryo (20 sections/z = 0.54 μm).

At no time did we detect any HCCs in the FL where HSPCs expand but do not emerge (Houssaint, 1981, Moore and Metcalf, 1970). We observed mKO2 expression in the para-aortic mesenchyme (Figures S2A and S2B), consistent with the expression of endogenous Runx1 in sub-endothelial aortic mesenchyme (North et al., 1999, Tober et al., 2013) and stromal cells of AGM region (Nef et al., 2005, North et al., 1999). Compared to previous studies of the Runx1+23 enhancer (Table S1) (Bee et al., 2010, Ng et al., 2010, Nottingham et al., 2007), our reporter construct appears to better avoid non-specific expression in the somites, possibly due to the inclusion of insulators.

The Reporter DP LK Population Marks Immunophenotypic HECs, Pro-HSCs, Pre-HSCs I/II, and FL HSPCs

To quantitatively assess reporter activity during definitive HSPC development, we performed flow cytometric analyses of the major hematopoietic tissues. In the E9.5 PSp, non-HECs that express VEC but not Kit, CD41 or CD45 lacked mKO2 expression altogether. Runx1-mKO2 reporter emerged in VEC+Kit+CD45− cells, with CD41 levels correlating well with mKO2 expression levels, thereby corroborating the notion that the Runx1-mKO2 reporter faithfully marks HECs and EHT activity. mKO2 reporter expression was maintained in immunophenotypic E9.5 pro-HSCs (VEC+Kit+CD41loCD45−) (Figure 2A) and in E10.5–E11.5 type-I (Lin−VEC+Kit+CD41lo(CD43+)CD45−) and type-II (Lin-VEC+Kit+CD41lo(CD43+)CD45+) pre-HSCs (Figure 2B) (Rybtsov et al., 2014). Of note, the pro- and pre-HSC compartments were markedly heterogeneous for Ly6a-GFP expression, with the ratio of DP to R-SP cells increasing as embryos mature (Figures S3A–S3C).

Figure 2

Runx1-mKO2 and Ly6a-GFP Expression in Pro-/Pre-HSCs and FL SLAM HSCs

(A) Flow cytometric analysis of Runx1-mKO2 expression in VEC+Kit−CD41−CD45− pre/non-hemogenic endothelial cells, co-induction of mKO2 and CD41 in the VEC+Kit+CD45− cell compartment (middle), and Runx1-mKO2/Ly6a-GFP expression in pro-HSCs (right) of the E9.5 para-aortic splanchnopleura (two independent experiments, n = 3).

(B) Flow cytometric analysis of Runx1-mKO2/Ly6a-GFP expression in the E10.5 (left) and E11.5 (right) immunophenotypic pre-HSC I/II compartments of AUV, yolk sac, placenta, and fetal liver (more than two independent experiments, n = 5 or more).

(C) Flow cytometric analysis of Runx1-mKO2/Ly6a-GFP expression in E13.5/E14.5 fetal liver LSK-SLAM HSCs and progenitors. Live (DAPI−) Lin−Sca-1+Kit+-gated cells are subdivided into four populations based on their expression of CD48 and CD150: HSCs (CD150+CD48−), MPPs (CD150−CD48−), HPC-1 (CD150−CD48+), and HPC-2 (CD150+CD48+). Expression of mKO2 in each population are analyzed on the right. (C′) Relative mean fluorescence intensity of mKO2 expression in each HSPC populations (five independent experiments, n = 16). Asterisks indicate significant differences (∗∗∗p < 0.001) compared with HSC according to Dunnett's multiple comparison test. Bars and values represent mean, ± σ (FACS plots) or SEM (bar graphs).

Next, we analyzed reporter expression in highly enriched E13.5 FL stem- and progenitor populations using the Lin−Sca-1+Kit+-signaling lymphocytic activation molecule (LSK-SLAM) purification scheme (Kent et al., 2009, Kiel et al., 2005, Oguro et al., 2013). LSK cells were gated on HSCs (CD150+CD48−), multi-potent progenitors (MPPs) (CD150−CD48−), and heterogeneous restricted progenitors (CD150−CD48+ HPC-1 and CD150+CD48+ HPC-2 cells). Consistent with the Sca-1+ gating, all four populations were Ly6a-GFP+ (Figure 2C). mKO2 expression levels were noticeably higher in the immunophenotypic SLAM HSC compartment compared with SLAM MPPs and SLAM HPCs (Figure 2Cʹ). Under conventional culture conditions, HSCs quickly differentiate into progenitor cells (Dykstra et al., 2006, Kent et al., 2008), resulting in a rapid decrease in their HSC frequency. Consistent with the expected quick differentiation of HSCs into progenitor cells under conventional in vitro culture conditions, Runx1-mKO2 levels in cultured FL SLAM HSCs rapidly dropped down to those seen in SLAM HPC-I progenitors (Figure S3D).

Erythro-Myeloid and Lymphoid Potential Progenitor Cells Are Enriched in Reporter DP HPCs

We used in vitro erythro-myeloid colony-forming unit cell (CFU-C) assays to investigate how hematopoietic potential correlates with the activity of the reporter genes. E10.5–E11.5 AUV, YS, PL, and FL Lin−Kit+ (LK) cells were separated into four subpopulations: DP, R-SP, L-SP, and double-negative (DN) by fluorescence activated cell sorting (FACS). (Figures 3A and S3E; Table S2). The immunophenotype of both R-SP and DP LK cells resembles that of pre-HSC I/II and HSCs (Figures S3C and S3E). Independent of the tissue source, both DP and R-SP populations were highly enriched for cells with erythroid and myeloid potential, with 10%–20% of the cells forming hematopoietic colonies in vitro (Figures 3B and 3C).

Figure 3

Runx1-mKO2 and Ly6a-GFP DP cells purified from different major hematopoietic tissues show similar erythroid-myeloid and lymphoid potential.

(A) Scheme for isolating Lin−Kit+ mKO2+GFP+ DP, mKO2+GFP−, (R-SP), GFP+mKO2− (L-SP), and mKO2−GFP− (DN) populations from E10.5 and E11.5 embryos for hematopoietic potential assays, the frequency of subpopulations are indicated in the flow plots (mean ± σ) and cell number collected per embryo (mean, SEM) are shown in the histograms and Table S2. Number of erythroid-myeloid CFUs (CFU-Cs) per 1,000 sorted cells of the indicated tissue from E10.5 (B) and E11.5 (C). E, burst-forming units– erythroid; GM, granulocyte/monocyte colony-forming units; Mix, granulocyte/monocyte/erythroid colony-forming units. Numbers are compiled from five to seven E10.5 or E11.5 conceptuses per population, five to six independent experiments. Error bars represent SEM. Asterisks indicate significant differences (∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001) compared with DN according to Dunnett's multiple comparison test. Differences between DP and R-SP were determined by unpaired two-tailed Student's t test. n.s., not significant. (D) Frequency of LPs in the DP and R-SP fractions of cells sorted from E11.5 embryos (mean ±95% confidence interval). B cells were identified as CD19+B220+ (left panel) and T cells as CD25+Thy1.2+ (right panel). Progenitor frequencies are shown above each column. Asterisks indicate significant differences (∗∗∗p < 0.001) between pairwise tests. Data are from three independent experiments using pooled cells.

To determine the frequencies of LPs with B or T cell potential we performed in vitro limiting dilution assays on OP9 or OP9-DL1 stromal cells (Schmitt et al., 2004). Low total live cell numbers after sorting precluded definitive analyses of the AUV and YS (Table S2). Nevertheless, with E11.5 PL and FL cells we observed LPs at very high frequencies in the DP population, with approximately one in three DP cells having B and/or T cell potential (Figure 3D). Of note, R-SP-sorted populations did occasionally become DP in vitro and were also able to give rise to lymphoid cells, albeit at much reduced frequencies. In contrast, L-SP and DN populations were completely devoid of lymphoid activities (data not shown). The strong enrichment of functional progenitors in the mKO2 fraction, combined with their near-absence in mKO2− fractions, validates the specificity of our reporter system; moreover, the correlation of lymphoid (but not myeloid) potential with GFP+ expression supports previous reports showing that LPs are confined to the Ly6a-GFP+ compartment and depend on Runx1/CBFβ function (Chen et al., 2011, Li et al., 2014).

HSC Activity Is Restricted to DP LK Cells, with the Most Robust HSC Activity Present in Placenta by E11.5

Next, we assessed the HSC activity of the four FACS-segregated populations (as described in Figures 3A and S3E) by transplanting them into irradiated adult recipient mice. Recipients with more than 5% reconstitution and contribution to both lymphoid and myeloid lineages were scored as successfully engrafted (Figure 4A). Using this stringent threshold, we did not detect engraftment with cells isolated from E10.5 embryos, presumably due to their extreme rarity, consistent with prior reports (Arora et al., 2014, Medvinsky et al., 2011). However, by E11.5, functional HSCs capable of generating robust (>10%) multi-lineage donor contribution were present in the DP compartments of YS, PL, and FL (Table S3), while R-SP, L-SP, and DN cells contained little if any HSC activity. We failed to detect donor cell contribution above 5% after direct transplantation of cells from any of the E11.5 AUV compartments; however, this was not entirely unexpected as HSCs are still rare in the E11.5 AGM (Taylor et al., 2010), with those that exist often generating contribution levels in adult recipients that fall below our 5% threshold (Arora et al., 2014). Nevertheless, we were able to detect E11.5 AUV-derived HSC activity after allowing the cells to maturate/expand ex vivo before transplantation (Table S3).

Figure 4

Placenta Runx1-mKO2 and Ly6a-GFP DP Cells Contain the Most Robust HSC Activity in E11.5 Embryos

(A) A representative peripheral blood assay of a recipient with robust multi-lineage engraftment by donor-derived (CD45.2+CD45.1−) HSCs.

(B) Lin−Kit+ (LK) DP, R-SP, L-SP, DN cell populations were sorted from PL and FL of E12.5 embryos and (1 embryo equivalent) transplanted to irradiated adults. Engraftment was evaluated by flow cytometric detection of multi-lineage donor-derived cells (>5%, dotted red line) in recipient peripheral blood at least 16 weeks after transplantation. Percent donor-derived blood is plotted; each dot is an individual transplant recipient. Error bars represent SEM (five independent experiments, n = 6–7).

(C) Limiting dilution analysis (LDA) for E11.5 placenta and fetal liver Lin−Kit+ DP.

(D) LDA of HSC frequency in E13.5 fetal liver LT HSC (LSK-SLAM) (two independent experiments, n = 3–11) (blue), E11.5 placenta (red), and fetal liver (black) (four independent experiments, n = 8–10). The slope of the line represents the log-active cell fraction. The dotted lines indicate 95% confidence interval. Engraftment failures are represented by down-pointing triangles. Calculated HSC frequencies and confidence intervals are shown above the plot.

Functional HSCs of E12.5 FL and PL remained almost uniformly positive for both reporters, with HSCs from PL engrafting much more robustly (higher average donor cell contribution per embryo equivalent) compared with FL (Figure 4B). Furthermore, limiting dilution analysis showed that HSC frequencies are much higher in the DP LK compartment of PL compared with FL at E11.5 (Figures 4C and 4D), contrasting with the much more similar frequencies of erythro-myeloid and lymphoid potential cells (Figures 3C and 3D). In summary, shortly after they emerge, functional HSCs were found primarily in the E11.5–E12.5 DP LK fraction of the PL, both on a per embryo and a cell frequency basis.

AUV/YS Pre-HSCs, PL HSCs, and FL LT HSCs Transcriptomes Are Highly Similar in Gene Regulatory Networks-Based Transcriptome Analysis

To gain insights into the transcriptional dynamics during hematopoietic commitment and maturation, we performed microfluidic-based single-cell RNA sequencing (RNA-seq). We analyzed DP LK cells isolated from E10.5 AUV and YS (emerging HCCs containing pre-HSC I), E11.5 PL (enriched for pre-HSC II/HSCs), as well as E13.5 FL HSCs (enriched using the LSK-SLAM HSC purification scheme). After applying stringent quality controls (see Experimental Procedures), 32 AUV, 28 YS, 36 PL, and 31 FL cells were kept, with read counts ranging from 1 to 3 × 106 per cell. Each population was profiled using CellNet, a computational platform that combines known gene expression data and gene regulatory networks to quantitatively deduce the cell or tissue type of a query sample from gene expression data (Cahan et al., 2014, Radley et al., 2017). CellNet classified virtually all DP LK cells and HSCs as “HSPCs” (Figure 5A). Unsupervised clustering of the cells using DBSCAN (Ester et al., 1996) with t-distributed stochastic neighbor embedding analysis (van der Maaten and Hinton, 2008) identified three subgroups of cells (Figures 5B and 5C): group 1 (primarily E10.5 AUV and YS DP LK cells), group 2 (mostly E11.5 PL DP LK cells), and group 3 (E13.5 FL LT HSCs). To determine whether there were differences in hematopoietic lineage identity between the clusters, we used SingleCellNet (SCN), which was trained on a compendium of mouse cell types (Tan and Cahan, 2019). We found that most cells in all three groups had high SCN scores for common lymphoid progenitor (CLP), granulocyte monocyte progenitor (GMP), and/or hematopoietic precursor cell (Figures 5D and 5E). However, cluster three had the highest proportion of cells classified as hematopoietic precursor cells, whereas clusters 1 and 2 had higher numbers of cells classified as CLPs and GMPs (Figure 5F). A small number of cells were classified as endothelial cells in all three clusters.

Figure 5

CellNet Analysis of Single-Cell Transcriptomes

(A) The classification heatmap of E10.5 AUV, YS, and E11.5 PL DP LK cells and E13.5 FL LT HSCs using CellNet.

(B) t-distributed stochastic neighbor embedding (t-SNE) analysis of single-cell transcriptomes.

(C) 127 single-cell transcriptomes were segregated into three cell groups by DBSCAN clustering.

(D) The classification heatmap of the three cell groups using SingleCellNet (SCN).

(E) The four major cell types SCN scores of the three cell groups. (F) The proportions of the major classified cell types of the three cell groups.

Waterfall Pipeline Analysis Reveals the Interferon Signaling and Prolonged Cell-Cycle Signatures during Fetal (Pre-)HSC Maturation

To gain further insights into the transcriptional dynamics during HSC emergence and maturation, we next applied Waterfall pipeline analysis (Shin et al., 2015), which identifies a developmental pseudotime trajectory and aligns the individual cells along this axis. Five main clusters emerged (Figures 6A and 6B), with the main pseudotime trajectory moving from cluster 5 through 4 to cluster 1 (Figure 6C). The asterisks in Figure 6C mark two E11.5 PL “outlier” cells that both have a megakaryocyte-like gene expression profile (Figure S4) and that may thus be related to the previously described megakaryocyte lineage biased HSCs (Sanjuan-Pla et al., 2013). The main pseudotime trajectory (constructed in an unsupervised manner), closely mirrors a developmental progression from E10.5 AUV and YS through E11.5 PL to E13.5 FL LT HSC (Figure 6D). Spearman correlation analysis identified genes and gene sets whose expression changes most significantly correlate with this pseudotime trajectory (Table S4). Among 8,672 genes detected, only 190 (2%) genes showed statistically significant positive correlation with the pseudotime trajectory, while the expression of 5,143 (59%) genes correlated negatively (Figure 6E upper). A heatmap of the expression of the 10 most positively and the 10 most negatively correlated genes is shown in Figure 6D.

Figure 6

Waterfall Analysis of Single-Cell Transcriptomes

(A) Visualization of cell-cell relationships using principal component analysis (PCA).

(B) The first three principal components were modeled as multivariate normal mixtures and automatically determined the number of clusters by using the model parameters that maximize the Bayesian Information Criteria, as implemented in the mClust R package.

(C) Waterfall builds minimum-spanning tree connecting the cluster centers in the PCA space and predicts the existence of one main trajectory (red line) that encompasses clusters 5, 4, and 1 and projects cells along the pseudotime coordinate. Asterisks indicate two outlier cells showing a megakaryocytic expression signature (see Figure S4).

(D) Genes directly involved with the trajectory progression were ranked by computing a Spearman rank correlation between the expression of each genes and the trajectory progression. The heatmaps show the top 20 genes most highly correlated or anti-correlated with the trajectory progression.

(E) Proportions of genes whose expression positively or negatively correlates with the pseudotime trajectory (upper) and Mouse Genome Informatic classification of the positive correlated genes (lower).

(F) Top 10 gene ontology biological process terms of correlated and anti-correlated genes along the pseudotime.

(G) Top 10 Reactome pathways analysis of the anti-correlated genes along the pseudotime.

(H) Cell-cycle phase assigned using CellRouter calculated cell-cycle phase scores.

(I) Positively correlated cell surface molecules and expression profiles of two representative examples, Robo4 and Esam. Each data point represents the gene expression level of a single cell with a color scheme shown in legends. Data points are fitted with local polynomial regression fitting (blue lines) with 95% confidence interval (gray area).

(J) Positively correlated transcription factors and representative expression profiles of two representative examples, Pou3f2 and Mecom.

The most significantly upregulated gene sets identified by gene ontology (GO) and Reactome pathway analysis (Figures 6F and 6G; Tables S5 and S6) were ribosome biogenesis (a potential Runx1 target [Cai et al., 2015]) and nonsense-mediated decay (essential for HSPC maintenance [Weischenfeldt et al., 2008]). Other signatures detected include negative correlations of pre-HSC maturation with cell cycle, Rho-GTPase signaling, organelle biosynthesis, RNA expression, and metabolism-related gene sets. We used CellRouter to calculate cell-cycle phase scores and assign a cell-cycle phase to each cell (Lummertz da Rocha et al., 2018). Projecting these cell-cycle phase assignments onto the Waterfall trajectory showed that about half of the pre-HSCs/HSCs in cluster 1 are in G1 phase, while the remainder is in a more proliferative state. This is in contrast with cells at the beginning of this trajectory (clusters 5 and 4), which appear to be much more proliferative (Figure 6H). These overall trends are potentially reflective of the downregulation of these pathways that is thought to occur as HSCs acquire a quiescent phenotype, typically associated with postnatal HSCs (Manesia et al., 2015). Indeed, we found >95% of E11.5 FL pre-HSC-I and pre-HSC-II cells to be in G1 or S/G2/M phase while E13.5 FL and adult BM SLAM HSCs include sizable G0 populations of 23% and 89%, respectively (Figures S5A and S5B). Furthermore, we saw significantly higher ex vivo division rates among mKO2high cells from E11.5 AUV (21.3%) and E11.5 PL (20.6%) compared with E14.5 FL SLAM HSCs (9.94%) (Figure S5C).

Upregulation of Surface Marker and Transcription Factor Genes during HSC Emergence and Maturation

More than one-third (72/190) of the genes whose expression we found to positively correlate with (pre-)HSC maturation have known hematopoietic or immune functions, including regulation of multi-lineage potential (Figure 6E, lower). Cell surface molecules that become upregulated (Figure 6I; Table S7) include Ly6a, Robo4 (involved in HSC homing and engraftment [Smith-Berdan et al., 2011]), and the HSC markers Cd93, Notch2, and ESAM (Petrenko et al., 1999, Varnum-Finney et al., 2011, Yokota et al., 2009). Furthermore, Abcg1, Ptprj, and Ptprs are known regulators of HSC proliferation and mobilization, Flk2/Flt3 signaling, and repopulating capacity, respectively (Arora et al., 2011, Quarmyne et al., 2015, Westerterp et al., 2012, Yvan-Charvet et al., 2010). The list of upregulated genes also includes several interferon-induced genes (Ifitm1, Ifitm2, Ifitm3, and Ly6a). Interferons have the ability to regulate HSPC proliferation (Essers et al., 2009, Raefsky et al., 1985) in the adult bone marrow and have more recently been implicated in promoting the developmental emergence and maturation of HSCs (Kim et al., 2016, Li et al., 2014, Samokhvalov et al., 2007).

Sixteen of the 190 upregulated genes are transcription factors (Figure 6J; Table S7), several of which have known HSC function in the adult and/or embryo, including Gfi1 (Hock et al., 2004), Smarca2 (Holmfeldt et al., 2016), Foxo3 (Miyamoto et al., 2007, Tothova et al., 2007), Pbx1 (Ficara et al., 2008), and Mecom (Goyama et al., 2008, Kataoka et al., 2011, Zhang et al., 2011). In summary, these findings validate the specificity and confirm the utility of the Runx1 and Ly6a dual reporter approach in the tracking, isolation, and analysis of the transcriptional dynamics of pre-HSC to HSC maturation.

Dual Reporter Mouse-Derived Embryonic Stem Cells Show Hemato-Endothelial Potential but Fail to Generate DP HSC-like Cells

Next, we generated dual reporter embryonic stem cells (ESCs) to enable time course imaging studies. Embryoid body-based hematopoietic in vitro differentiation readily generated VEC+CD31+ endothelial cells, including R-SP and L-SP cells (Figure S6A). Time course imaging confirmed that R-SP putative HECs rapidly undergo EHT to form mKO2bright HCCs that bud off from the endothelial cell layer (Figures S6B and S6C). However, none of the Runx1-mKO2+ HECs or HCCs co-expressed Ly6a-GFP (Figures S6A–S6D), consistent with the lack of production of cells with HSC potential under these in vitro differentiation conditions. As is the case in vivo (Figure 2A), Runx1-mKO2 expression correlated well with the induction of CD41 expression (Figure S6E), whereas Runx1-mKO2− cells did not express any CD41 regardless of their Ly6a-GFP expression status (Figures S6E and S6F). Taken together, the dual reporter ESC model confirms that the Runx1-mKO2 reporter marks HECs, and the paucity of DP HEC/HCC formation illustrates the inability of conventional protocols to generate ESC-HSCs, making it a potentially useful system for screening and identification of missing factors or pathways.

Discussion

We used two of the best-studied HEC/HSPC regulatory sequences, the Runx1+23 HSPC enhancer (Bee et al., 2010) and the Ly6a promoter (de Bruijn et al., 2002), to generate a novel dual reporter mouse line and carefully validated the specificity and utility of this system using a battery of tests, including 3D confocal immunofluorescence, flow cytometry, progenitor cell assays, transplantation studies, and single-cell transcriptome analysis.

The observed activity of the reporter genes during EHT, HCC formation, and HSC maturation is consistent with the model shown in Figure 7. The first sign of a forthcoming emergence of definitive-type HSPCs is the onset of Runx1+23 enhancer activity in HECs of the hemogenic vessels around E8.5 (Swiers et al., 2013). Shortly thereafter, EHT commences and the first Kit+ HCCs begin to bud off (Yokomizo and Dzierzak, 2010). HCCs only emerge from Runx1+23 enhancer positive HECs (Ng et al., 2010). However, while myeloid progenitor cell potential can be found in all HCCs populations, lymphoid and HSC fates are largely restricted to the Runx1+23/Ly6a DP fractions. Our system has several advantages over other existing approaches for visualizing and isolating specific cell populations (1) emerging fetal HSCs can be distinguished from contaminating maternal HSCs; (2) it enables antibody-free live imaging studies; (3) the strong nuclear mKO2 reporter provides a brighter and more focused signal that facilitates cell identification in fluorescent microscopy; and (4) it labels functionally distinct populations more precisely than other current reporter systems (e.g., combining Ly6a-GFP and CD45 is insufficient to enrich pre-HSCs or myeloid progenitors in Kit+ HCCs [Boisset et al., 2015]).

Figure 7

Schematic Summary of Transcriptome Dynamics of Pre-HSCs to HSCs Maturation by Tracking Runx1-mKO2 and Ly6a-GFP Reporters in Mouse Embryo

At E10.5 of mouse development, Kit+ HCCs emerge from Runx1-mKO2+ HECs through EHT. These include definitive myeloid-restricted (Ly6a-GFP−) as well as pre-HSCs/lymphoid potential cells (Ly6a-GFP+). In contrast, the placenta harbors few HECs, yet it is the first site that at E11.5 contains numerous robustly engrafting Runx1-mKO2/Ly6-GFP DP HSCs, suggesting that it is the major niche for (pre-)HSC maturation and expansion before HSCs migrate to the fetal liver. Single-cell RNA-seq with Waterfall and Reactome pathway/GO analyses reveals the transcriptional landscapes and dynamic changes in biological processes during pre-HSC to HSC maturation. During maturation, mitosis/cell cycle and biosynthesis genes become downregulated while regulators of multi-lineage differentiation become upregulated. In addition, the transcriptomes of maturing cells suggest they have encountered interferon signals.

Bioinformatics analyses of the changes in gene expression patterns of DP cells during development allowed us to identify cell surface and transcription factor genes that are potentially involved in the maturation of HSCs. We show that during pre-HSC maturation, interferon-inducible genes are among the most highly induced genes, supporting a recently proposed model implicating placental interferon signaling in pre-HSC maturation (Kim et al., 2016, Li et al., 2014). Another induced gene is the homing receptor Robo4, and its ligand Slit2, which are likewise expressed in the mid-gestation placenta (Duc-Goiran et al., 1994, Liao et al., 2010, Platt and Hunt, 1998), further suggesting a unique function of this niche.

Our Reactome and GO analyses also revealed a conspicuous downregulation of cell-cycle signature genes that occurs as DP pre-HSCs mature into HSCs and colonize the FL. The acquisition of a quiescent HSC phenotype was previously thought to be primarily a postnatal phenomenon (Manesia et al., 2015). However, our cell cycle and ex vivo cell division studies showed that as E11.5 HSC precursor populations mature into E13.5–E14.5 SLAM HSCs, they become less proliferative and give rise to a significant G0 population, albeit at a still much reduced frequency relative to bona fide quiescent BM LT SLAM HSCs. Our data are, therefore, consistent with recent reports showing that functional FL HSCs are already predominantly in the G0/G1 phase of the cell cycle (Batsivari et al., 2017, Zape et al., 2017).

A long-standing but still unanswered question is why conventional pluripotent stem cell in vitro differentiation methods generally fail to produce HSCs. In vitro differentiation of our dual reporter mouse ESCs produced endothelial cells, HECs, and HCCs, but failed to generate DP cells under myeloid in vitro differentiation conditions. It will be interesting to see if DP cell production emerges and correlates with lymphoid or HSC development when using conditions that can unlock these fates (Sugimura et al., 2017, Vo et al., 2018).

In summary, our validated reporter system will prove a useful tool for dissecting the still incompletely understood processes that regulate definitive hematopoiesis, including specification of HECs and lymphoid/HSC potential, EHT, and maturation of pre-HSCs into HSCs. In addition, the image analysis friendly nature of our reporter system may be leveraged in large-scale screens for agents that enable more efficient production of HSCs or cells with lymphoid potential in vitro for modeling or therapeutic purposes.

Experimental Procedures

All animal procedures were conducted in accordance with the animal care guidelines approved by the Institutional Animal Care and Use Committee at Boston Children's Hospital (protocol no. 15-12-3060R). Staged mouse embryos were obtained using timed pregnancies of B6129SF1/J females with Runx1+23-H2BmKO2;Ly6a-GFP hemizygote males. Tg(Runx1+23-H2BmKO2) was generated by replacing GFP with H2BmKO2 in the construct +23 GFP reporter mouse (Bee et al., 2010) and flanked by insulators. Confocal microscopy, flow cytometry/FACS, erythro-myeloid CFU-C, and LP assays, and transplantation assays were performed as described previously (Chen et al., 2011). Fluorescent microscopy raw images were processed using FIJI rolling ball background subtraction, log transform, and pseudo-coloring functions. Single-cell whole transcriptome amplification was performed using Fluidigm C1 System (Fluidigm). Samples were sequenced using a HiSeq 2500 sequencer (100 bp paired-end reads). Single-cell samples were sequenced to an average depth of 1.5 million reads per cell. Sequence reads were aligned to the mouse transcriptome using HISAT version 0.1.6.

Author Contributions

M.J.C. designed, performed, analyzed, and interpreted the experiments. C.K., P.H., P.S., N.K.M., Y.F., M.N., S.L., and Y.-F. L. performed and assisted with the experiments. E.L.d.R., P.C., Y.T., A.H., and S.Y. mapped, analyzed, and interpreted the single-cell data. L.V., D.K.J., Y.Z., and T.E.N. provided support and advice. M.J.C., L.I.Z., G.Q.D., and T.M.S. designed the study. T.M.S. performed image analyses. M.J.C., E.L.d.R., P.C., T.E.N., L.I.Z., G.Q.D., and T.M.S. wrote the manuscript.