Single-cell multiomics of human fetal hematopoiesis define a developmental-specific population and a fetal signature

Knowledge of our fetal blood system can facilitate the generation of transplantable hematopoietic stem cells (HSCs) derived from pluripotent stem cells (PSCs) for future regenerative medicine¹ and also increase our understanding of congenital blood disorders and pediatric leukemia.² 

Fetal and adult blood cells differ in cell cycle status, molecular profile, and composition of progenitor cells, and hematopoiesis takes place in different niches depending on the stage of development.^3-6 Moreover, some innate immune cells are formed almost exclusively during fetal life.⁷ The blood system emerges in waves and is initiated by HSC-independent hematopoiesis in the yolk sac, where immature erythroid and myeloid cells, critical for the growing embryo, emerge, followed by erythromyeloid and lymphomyeloid progenitors.^5,7-9 From around day 27 (∼Carnegie Stage [CS] 13), hematopoietic stem and progenitor cells (HSPCs) are formed in the aorta-gonad-mesonephros region and later migrate to the fetal liver (FL), where they mature and expand. The bone marrow (BM), the site of blood production in adults, starts to be colonized at the end of the first trimester.⁵ 

The different waves and niches add to heterogeneity, and the sparsity of human fetal samples makes studies in early development demanding. Much of our knowledge regarding embryonic and fetal hematopoiesis has been translated from studies of the murine system, but accumulating reports of the human counterpart complement the current model.^5,10 Of late, single-cell studies have emerged focusing on embryonic/fetal development, increasing our knowledge of blood formation during human ontogeny.^11-13 However, direct comparisons between human fetal and adult hematopoiesis are just beginning to be explored.^6,14,15 Immunophenotypic markers used to purify and prospectively analyze progenitors may not be equally expressed throughout development, making direct comparisons of purified populations less robust. Furthermore, the analysis of transcriptomic data from different niches and stages also poses challenges, as heterogeneity may be lost upon integration.

Here, we set out to directly compare the transcriptome and immunophenotype of embryonic/fetal and adult hematopoiesis using CITE-seq (cellular indexing of transcriptomes and epitopes by sequencing).¹⁶ To better preserve heterogeneity in cells from different niches and stages of development, a projection approach was used.¹⁷ This allowed for a comparison of progenitor composition within the fetal and adult HSPC compartment. The projection approach revealed a multipotent progenitor cluster, which was hardly detectable in the adult BM. Investigating a publicly available data set,⁶ these cells were specifically found in the FL niche. Furthermore, by comparing molecular signatures across developmental stages, we could define a fetal core signature, with genes specifically upregulated or downregulated in fetal cells compared to adult. By using a publicly available RNA-seq data set from patients with acute lymphoblastic leukemia (ALL),¹⁸ the fetal core signature was found to define some cases of pediatric leukemias, thus enabling the separation of specific subgroups of ALL.

In summary, our comparison between different stages of development highlights important transcriptomic and immunophenotypic differences in human fetal and adult primitive blood cell development.

Human FLs were donated from elective medical terminations of pregnancy after informed consent and with the approval of the Regional Ethics Review Board (Lund)/the Swedish Ethical Review Authority and the Swedish National Board of Health and Welfare. Embryos/fetuses were staged according to CS and postconceptional weeks (pcw) and were all from the first trimester of pregnancy (developmental stages CS16-9 pcw) (supplemental Table 1). The FL was mechanically disrupted and dissociated through a 40 μm filter. The cells were frozen in StemCellBanker (Amsbio) or in fetal bovine serum with 10% dimethyl sulfoxide. Samples were stored at –150°C. Cord blood (CB) (n = 3) and adult BM (20 to 30-year olds, males; n = 2) data sets were from a previous study.¹⁹ 

FL samples were thawed on the day of experiment and stained with a panel of CITE-seq antibodies (supplemental Table 2), a hashtag oligonucleotide (HTO) antibody to facilitate sample multiplexing, a Fc receptor blocking reagent (Miltenyi Biotec), CD45-Alexa700 (HI30, Biolegend), CD34-FITC (581, Biolegend) and lineage markers (LIN: CD19-BV605 (SJ25C1, BD Bioscience), CD3-PE-Cyanine5 (UCHT1, Biolegend), CD2-PE (RPA-2.10, Biolegend), CD16-BV421 (3G8, Biolegend), CD14-PE-Cyanine7 (M5E2, Biolegend), and CD235a-PE-Cyanine5 (GA-R2, BD Bioscience)). These were incubated at 4°C for 30 min. Samples were then washed and resuspended in phosphate-buffered saline with 2% fetal bovine serum (GE Lifesciences) and 7-aminoactinomycin D (7-AAD) (BD Bioscience). Up to 10 000 LIN^-CD45⁺CD34⁺ cells were sorted from each donor using a BD FACSAria IIu. The CB and adult BM data sets from¹⁹ were generated in a similar way from LIN^–CD34⁺ cells.

After sorting, the samples were loaded on to the 10× genomics 3’ version 2 platform (10× Genomics) and single-cell RNA-seq was performed according to the manufacturers’ instructions with minor changes according to a previous study¹⁶ to allow for CITE-seq. After reverse transcription and cDNA amplification, the resulting libraries were sequenced on a NOVAseq (Illumina). After sequencing, the binary base call files were processed using Cellranger mkfastq to produce FASTQ files, and Cellranger count was used for alignment to the human reference genome (Hg38), filtering, barcode counting, and unique molecular identifier counting. The resulting matrix files were further analyzed in sequential steps using Seurat²⁰ and Scarf.¹⁷ 

The output from Cellranger was processed using Seurat (V4).²⁰ Samples were filtered based on unique molecular identifier, gene, and mitochondrial counts to remove low-quality cells (supplemental Table 1). The cells were demultiplexed using the hashtag antibody, and doublets were removed. Each donor was treated separately. RNA reads were log-normalized and antibody derived tag (ADT) and HTO data were CLR (centered-log ratio) normalized using the NormalizeData function. CS22 FL samples were integrated using the IntegrateData function, and cell cycle effects were regressed out using the ScaleData function. The output data were used to make a unified uniform manifold approximation and projection (UMAP). The cells were clustered using FindClusters (resolution 0.7), which resulted in 11 clusters. These clusters were associated with a hematopoietic progenitor cell state based on the results from the FindAllMarkers function.

The scaled and normalized RNA and ADT data of FL CS22 were used to perform weighted nearest neighbor analysis (WNN).²⁰ The FindMultiModalNeighbors function of Seurat (V4) was used with the first 30 PCs of the RNA principal component analysis (PCA) and the first 18 PCs of the ADT PCA with default settings, after which RunUMAP was performed using WNN and default settings. The cells were clustered using FindClusters, with the smart local moving algorithm (resolution of 0.5), resulting in 10 clusters. The metadata and the UMAP coordinates were exported and used in Python for visualization.

The CLR normalized ADT values from Seurat were loaded into Python, transformed by their antilog (base e), and then multiplied by 1000 to enable gating in FlowJo. The values were exported using the write_FCS function and loaded into FlowJo V10 (BD), where conventional gating was performed using internal negative controls to set the gates. The gated populations were exported as csv, and using Python, the ADT values from cells in the exported gates were matched to the original data to identify the cell identities.

The ADTs were used to immunophenotypically define the following populations (gated from LIN^-CD34⁺CD45⁺ cells): HSC^CD49F: LIN^-CD34⁺CD38^-/lowCD45RA^–CD90⁺CD49F⁺; HSC^CD90: LIN^–CD34⁺CD38^–/lowCD45RA^-CD90⁺; MPP: LIN^–CD34⁺CD38^–/lowCD45RA^–CD90^–; LMPP: LIN^–CD34⁺CD38^–/lowCD45RA⁺CD90^–; CMP^CD123: LIN^–CD34⁺CD38⁺CD10^–CD45RA^–CD123⁺; CMP^CD135: LIN^–CD34⁺CD38⁺CD10^–CD45RA^–CD135⁺; GMP^CD123: LIN^–CD34⁺CD38⁺CD10^–CD45RA⁺CD123⁺; GMP^CD135: LIN^–CD34⁺CD38⁺CD10^–CD45RA⁺CD135⁺; MEP^CD123: LIN^-CD34⁺CD38⁺CD10^-CD45RA^–CD123^–; MEP^CD135: LIN^–CD34⁺CD38⁺CD10^–CD45RA^–CD135^–; CLP^CD10: LIN^–CD34⁺CD38⁺CD45RA⁺CD10⁺; CLP^IL7R: LIN^–CD34⁺CD38⁺CD45RA⁺IL7R⁺; CLP^CD7: LIN^–CD34⁺CD38⁺CD45RA⁺CD7⁺.

The projection of cells on the FL CS22 reference UMAP was done in Scarf.¹⁷ In brief, the UMAP coordinates, HVGs, and cluster information of the FL CS22 analysis in Seurat were loaded into Scarf, and by using the run_mapping function (with k = 11), the nearest neighbors of each cell projected were calculated. This was performed for FL CS16, FL CS22, FL 9 pcw, CB, and adult BM, the latter 2 from Sommarin et al,¹⁹ for immunophenotypically defined populations from Ranzoni et al,¹² and with data from Roy et al.⁶ For the data from Ranzoni et al,¹² CORAL was set to “True” in the run_mapping function to account for the difference in technology used to produce the data. With the adult BM as a reference, the same settings as for the FL CS22 were used, projecting FL CS16, FL CS22, FL 9pcw, CB, and adult BM cells. The projected cells were classified using the get_target_classes function (threshold = 0.4). The mapping score was calculated using get_mapping_score, and the size of each cell in the reference map was set proportional to the mapping score.

The FL CS22 and adult BM data sets were integrated with Harmony (version 0.1.0).²¹ The FL and BM Seurat objects were merged (merge function). The data were log-normalized, and highly variable genes were calculated with the FindVariableFeatures function. Cell cycle effects were regressed out using the ScaleData function. PCA was calculated using the RunPCA function, after which the RunHarmony function was used. The Harmony reduction was then used in FindNeighbors (dimensions 1:15), FindClusters (resolution = 0.8), and RunUMAP (dimensions 1:15, 15 neighbors, minimum distance = 0.005, spread = 6). The metadata and the UMAP coordinates were exported and used in Python for visualization.

To analyze differentially expressed genes (DEGs) in the CS22 MPP-I cluster compared with HSC and MPP-II, respectively, the FindMarkers function of Seurat was used. DEGs (adjusted [adj] P value < .001; log2 fold change (FC) > |0.5|) and cells were hierarchically clustered using the clustermap function of Seaborn (Euclidean distances and Ward linkage).

CellRadar (using data from HemaExplorer²²) (available at https://karlssong.github.io/cellradar/) was used to investigate lineage affiliation and to create a radar plot using the min-max scaled median value of marker genes in each cluster (adj P value < .001; log2 FC > 0.25). Furthermore, the CellCycleScoring function in Seurat was used for cell cycle phase classification for the CS22 FL cells.

To perform DEG analysis, each cluster of the FL CS22 and the adult BM cells projected to a FL cluster were pseudobulked into 2 replicates per donor sample using the make_bulk function of Scarf.¹⁷ These pseudobulked populations were then exported as csv and loaded into DESeq2²³ using DESeqDataSetFromMatrix (design = ∼sample, with “sample” referring to FL CS22 or adult BM), after which DEG analysis was performed (DESeq function with default settings) between FL CS22 and projected adult BM for each of the CS22 clusters. Clusters DC-I and Ly-III were excluded from the analysis because of their low cell count (<60 cells). The resulting files were then analyzed for upregulated (fetal signature) and downregulated (adult signature) genes using the adjusted P value < .05 and log2 FC > |1| to make gene lists for each cluster. For validation purposes, DEG analysis was also performed in DESeq2 using pseudobulked samples without pseudoreplicates (n = 2 per developmental stage).

The DEseq2 normalized counts for FL CS22 and projected adult BM clusters (including pseudoreplicates; n = 4 per developmental stage) were exported for Gene Set Enrichment Analysis (GSEA 4.3.2)²⁴ to analyze hallmark gene sets (h.all.v7.5.1).²⁵ Gene sets with a false discovery rate q-value of <0.05 were considered significantly enriched/depleted (nr of permutations: 1000, permutation type: gene set). The normalized values for all the CS22 FL clusters (with pseudoreplicates, excluding DC-I and Ly-III) were additionally exported for GSEA, now focused on defining differences between MPP-I/HSC and the MPP-I/MPP-II clusters.

Pseudobulked data of the clusters and the ADT-defined immunophenotypic populations (pseudobulking done as above) were loaded into DESeq2 (DESeqDataSetFromMatrix, design = ∼1). After running the DESeq function, the data were transformed using the vst function, and the top 500 most variable genes were used in the prcomp function of R to perform PCA.

The significant genes (adj P value < .05; log2 FC > |1|) for each cluster of the DESeq2 differential testing (described above) were visualized using the Venn tool (http://bioinformatics.psb.ugent.be/webtools/Venn/). Upregulated and downregulated genes were displayed separately. Clusters DC-I and Ly-III were excluded. Genes found to be differentially expressed in all primitive, lymphoid, or myeloid clusters were included in a gene set named ‘fetal core signature.' DEGs were also generated without pseudoreplicates (adj P value < .05; log2 FC > |1|) as well as using FindMarkers in Seurat (adj P value < .05; log2 FC > |0.5|), and these were analyzed as above. The fetal core signature as well as the gene sets generated without pseudoreplicates and with Seurat were further validated using a publicly available bulk RNA-seq data set.¹⁵ The fragments per kilobase of transcript per million mapped reads (FPKM) values were log2-transformed (log2(FPKM +1)), and PCA was performed using the gene sets described above and for validation purposes with the top 500 variable genes and the top variable genes excluding genes in the fetal core signature.

The single-cell FL and adult BM data sets were combined by using the ZarrMerge function of Scarf. The genes specific for FL (upregulated in primitive, Ly, or My) were used to subset the combined FL-BM data set. Using the clustermap function of seaborn, the cells and genes were clustered using hierarchical clustering (Euclidean distances and Ward linkage). In addition, all HOX genes expressed in >10 cells were used to create a heatmap using the plot_marker_heatmap function of scarf. Furthermore, using Seurat, the data from Roy et al⁶ were log-normalized and scaled, after which the AddModuleScore function was used to calculate scores for the fetal and adult signature genes for each cell.

The fetal gene signature was investigated in an ETV6::RUNX1 iPSC model using publicly available data from a previous study.¹⁵ FPKM values were log2-transformed (log2(FPKM +1)), scaled, and centered using the scale function in R and a heatmap created from the identified universal upregulated fetal core genes (14 out of 15 genes were identified). The z-scored mean expression of the same genes was also calculated for each sample and visualized in a box plot.

To investigate the fetal gene signature in childhood ALL, a publicly available bulk RNA-seq data set of almost 2000 patients was used.¹⁸ The following cytogenetic abnormalities were selected: KMT2A::AFF1, KMT2A::MLLT1, KMT2A::MLLT3, BCR::ABL1, ETV6::RUNX1, and high or low hyperdiploidy. The HTSeq files of these were loaded into DEseq2 using the DESeqDataSetFromHTSeqCount function (design = ∼RNA_lib, to account for batch effects from different sequencing methods). Next, the DESeq function was run, followed by variance stabilizing transformation by the vst function, and batch effects were removed with removeBatchEffect from the Limma package. The fetal- and adult-specific gene lists were translated into hg37 Ensembl IDs using the getBM function of Biomart (557 Ensembl IDs were detected from the 629 fetal core genes). The genes were used for PCA on all leukemia subtypes above, pooled as well as individually, using the prcomp function with default settings. To validate whether the PC1 variation seen was due to a random effect, a number of Ensembl IDs matching the number of fetal core Ensembl IDs found in the data set were randomly selected after removing genes for which the average expression across the samples was lower than for the lowest-expressed fetal core gene, and a comparison PCA was performed (5 iterations). In addition, for leukemia subtypes where adult counterparts were represented, the fetal core genes were z-scored, and the mean values for upregulated and downregulated fetal core genes were plotted for each age group. For KMT2A::AFF1, the top 25 genes and the lowest 25 genes from PC1 were used to generate a heatmap, and a second heatmap was generated looking specifically at expression of different HOX genes.

To generate statistics for the contribution of the fetal core signature to the PCA, a bootstrapping approach was used. Here, random sampling of all genes was performed with 10 000 iterations. The variance explained by PC1 of each iteration was used to generate a histogram. The variance described by PC1 of the fetal core signature was compared with the distribution to generate a P value. Statistical analysis for the z-scored fetal core was calculated with ANOVA function of R (aov), followed by the TukeyHSD function to generate the adjusted P values.

The emerging blood system is a complex mixture of cells from different niches and waves acting simultaneously, generating a heterogeneity that necessitates the use of single-cell assays. In addition, immunophenotypic markers classically used to purify progenitors in CB and adult BM may not necessarily be conserved throughout development. To address these issues, we analyzed FL HSPCs using CITE-seq, a high-throughput single-cell RNA-seq method wherein information on immunophenotype is sampled simultaneously with the transcriptome, through oligo-barcode–tagged antibodies, ADTs.¹⁶ FL samples were obtained from medical terminations of pregnancy in the first trimester, ranging in age from CS16 to 9 pcw (supplemental Table 1). HSPCs were selected based on the expression of CD45 and CD34 and for lack of mature lineage markers (LIN^-) (Figure 1A; supplemental Figure 1A). The analysis focused on CS22, where 6732 cells from 2 embryos were captured for further analysis after quality control (supplemental Table 1). Eleven clusters were identified and assigned to a hematopoietic state based on the expression of lineage-associated genes and transcription factors (Figure 1B-C; supplemental Figure 1B). Accordingly, the UMAP showed an HSC cluster at the top, expressing genes like MLLT3 (protein AF9) and BST2.²⁶ Subsequent to HSCs, cells with a multipotent progenitor (MPP)-like signature were located, followed by a distinct separation into a megakaryocyte-erythroid progenitor axis (MEP), a granulocyte-monocyte progenitor axis (GMP), a dendritic cell-monocyte cluster (DC-Mono), and a lymphoid progenitor axis (Ly-I/Ly-II). The lymphoid axis showed a gradual maturation from Ly-I to Ly-II, and the Ly-II cluster was molecularly more primed toward the B lineage (Figure 1C). The Ly-III was T lineage primed but contained less than 60 cells, as did the DC-I (dendritic cell precursor) cluster. A cluster located below the MPP clusters had cell cycle–associated genes among the DEGs and was named cycling (Cyc) (Figure 1C).

To further investigate relatedness between the clusters, samples were pseudobulked (2 replicates per donor), and PCA was performed. As expected, the more primitive populations clustered together (HSC, MPP-I, and MPP-II), and the lymphoid and myeloid clusters were clearly separated, with the Ly-II cluster positioned further away from the primitive clusters than Ly-I (Figure 1D).

To validate the hematopoietic state assigned to each cluster, a publicly available data set with immunophenotypically defined fetal HSPCs was projected on the CS22 UMAP.¹² Immunophenotypic HSCs, MPPs, and CLPs mapped to associated clusters as expected. However, the purified CMPs and MEPs mapped to both MEP and GMP clusters, whereas the purified GMPs projected to the Ly-I as well as the DC-Mono clusters (supplemental Figure 1C).

Thus, CITE-seq captured the cellular heterogeneity in early blood development, and different progenitor populations could be defined based on the expression of lineage-associated genes. The data set can be viewed in an interactive format on the following platform: https://singlecell.broadinstitute.org/single_cell/study/SCP2022/fetal-liver-hematopoiesis#study-visualize.

We next sought to further investigate the heterogeneity observed in some of the immunophenotypically defined progenitors using the ADTs (Figure 2A). First, the transcriptome and the ADT data were integrated. The WNN²⁰ UMAP generated showed high similarity with the transcriptionally defined cell states (supplemental Figure 2A). Next, the abundance of the individual ADTs was investigated within each cluster (supplemental Figure 2B). The HSC cluster expressed CD90 and CD49F and lacked expression of CD38, all in agreement with markers used in conventional purifications of the HSC population (Figure 2B).^27,28 Furthermore, CD201 (EPCR), a marker shown to define fetal HSCs,²⁹ was also distinctly expressed in the molecular HSC cluster. The comparatively immature lymphoid cluster (Ly-I) expressed CD10 and CD45RA, but was low in CD38 expression, indicative of lympho-myeloid primed progenitors (LMPPs),³⁰ whereas the relatively more mature lymphoid cluster (Ly-II), in addition to CD10 and CD45RA, expressed CD38, corresponding to a common lymphoid progenitor (CLP) phenotype.³⁰ This cluster also had a distinct expression of the interleukin 7 receptor (IL7R), known to mark CLPs in the murine system.³¹ The MEP cluster was positive for CD71, as expected (Figure 2B).^14,32 However, CD7, CD123 (IL3R), and CD135 (FLT3), the latter 2 markers traditionally used to define GMPs and common myeloid progenitors (CMPs),^30,33 showed a broader expression pattern (Figure 2B).

Next, the ADT information was used to perform conventional immunophenotypic gating (supplemental Figure 2C; supplemental Table 3). Molecularly defined HSC, MPP, LMPP, and MEP clusters were well captured by the conventional adult immunophenotype (Figure 2C; supplemental Figure 2D).

Conventional CLPs are typically defined based on the expression of CD10^30,34 and using this gating strategy, most of the captured cells were indeed part of a lymphoid transcriptional cluster. However, almost 60% of the cells belonged to the lymphomyeloid Ly-I cell state, likely because of issues discriminating the cells based on CD38 expression. Therefore, other gating strategies were investigated to identify a transcriptionally purer lymphoid cell state. In CB, CD7, when combined with a lack of CD38 expression, has been shown to enrich for lymphoid progenitors.^35,36 However, in the FL, CD7 surface expression measured with the ADT antibody was virtually absent in the CD38-negative fraction (Figure 2B).

Using CD7 in combination with CD38 and CD45RA, we defined a heterogeneous CD7⁺ CLP-like population with cells from both lymphoid and myeloid clusters, in agreement with earlier observations.^35-37 Using IL7R³¹ in combination with CD45RA and CD38, almost 80% of cells captured belonged to the Ly-II cell state, and the remaining cells were defined as Ly-I. Thus, IL7R was by far the best surface marker to capture the Ly-II cell state during early development (Figure 2C).

Immunophenotypic CMPs and GMPs showed high heterogeneity as measured by the number of molecular clusters and their relative proportions captured within each immunophenotype (Figure 2C; supplemental Figure 2D). Thus, in immunophenotype-gated GMPs, 16% to 21% of the cells molecularly belonged to the Ly-II lymphoid cell state (Figure 2C; supplemental Figure 2D).^30,33 A PCA of pseudobulked immunophenotypically defined populations showed GMPs located in the vicinity of lymphoid populations, in agreement with the previous analysis (Figure 2D; supplemental Figure 2E). When combining the immunophenotypic populations with the molecularly defined clusters from Figure 1D, heterogeneity was seen for the myeloid GMP and CMP populations, whereas for most other populations, the immunophenotypic population correlated with the corresponding molecular cluster (Figure 2E; supplemental Figure 2F).

Overall, our analysis showed that immunophenotypic gating of HSCs and LMPPs captured populations that were molecularly fairly homogenous, whereas CMPs and GMPs were heterogenous.

Next, we set out to directly compare the fetal cells with our adult BM data set of HSPCs, analyzed with CITE-seq in a similar way¹⁹ (Figure 3A; supplemental Figure 3A). Looking at immunophenotypically defined populations using PCA showed that the populations were mainly separated by developmental stage (PC1; 51% of variance) (supplemental Figure 3B). To preserve heterogeneity within the sample, a nearest neighbor projection approach was next used.¹⁷ A reference map was built, and the investigated cells were mapped onto the reference, so only molecular differences within the reference sample were considered. Using FL CS22 HSPCs as a reference, 37% of the cells in adult BM mapped to the fetal HSC cluster, whereas only 7% mapped to this cluster in the FL CS22 sample. In contrast, only 6% of the adult cells mapped to the lymphomyeloid Ly-I cluster, compared with 15% in CS22 (Figure 3B-E). Furthermore, the MPP-I cluster, which represented about 17% of the CS22 FL cells, was hardly detectable in the BM sample (3% of total cells). To narrow down the time window in development when these changes in cell states occurred, CB-derived cells were analyzed, generating similar results as those from adult BM (Figure 3B-D).¹⁹ Investigating 2 more fetal time points, CS16 and 9 pcw (supplemental Table 1), 6% and 13%, respectively, of the fetal cells mapped to MPP-I, identifying MPP-I as a fetal–specific cell state (Figure 3B-D).

Intriguingly, when using adult BM as a reference, almost no FL cells mapped to the adult HSC clusters and few to the more primitive MPP clusters, whereas the adult MPP-III was enriched in the CS22 and 9 pcw developmental stages. Furthermore, adult MEP-I was enriched in all stages of the FL (supplemental Figure 3C-E). The fetal-associated MPP-I mainly mapped to MPP-III on the adult BM map (supplemental Figure 3F). Thus, transcriptionally defined adult HSCs were virtually absent in first-trimester FLs.

We next asked if the differences in cellular states observed could be detected with an integration approach. Integration of fetal CS22 and adult BM samples with Harmony²¹ showed that transcriptionally defined HSCs (Harmony cluster 0) clearly increased with developmental stage and a lymphoid cluster (Harmony cluster 8) decreased with age. Harmony clusters 1 and 2 were dominated by fetal cells and contained a large fraction of the fetal MPP-I cluster (supplemental Figure 4A-D).

To better characterize the fetal-specific MPP-I, the cells were compared with the HSC and MPP-II clusters. Among DEGs, an upregulation of MYC and heat shock–related genes was detected, and an erythromyeloid lineage affiliation was observed with CellRadar (Figure 3F-G). Looking at cell cycle scores, the MPP-I had an increased number of cycling cells compared with the HSC but less than the MEP cluster, thus matching a trajectory from HSC to MPP-I to MEP. Of note, in this analysis, the cluster named Cyc was mainly enriched for G2/M (supplemental Figure 4E). GSEA on pseudobulked clusters^24,25 showed upregulation of MYC targets, MTORC1 signaling, and unfolded protein response and oxidative phosphorylation gene sets in MPP-I compared with the HSC and MPP-II clusters, indicative of a more active metabolic state (supplemental Figure 4F).

To confirm the fetal specificity of the MPP-I cell state, a published data set that included samples from different developmental stages was used.⁶ Mapping these data onto our CS22 data showed enrichment of MPP-I in the early FL sample, whereas the fetal, pediatric, and adult BM were almost completely depleted of the MPP-I cell state (Figure 3H).

Taken together, these data show that the MPP-I cluster harbors transient, fetal-specific, erythromyeloid–primed multipotent progenitors.

To further understand how fetal and adult cell types differ, the projected adult cells were compared with their fetal counterparts using DESeq2 (Figure 4A). The clusters were pseudobulked for the analysis, and as with the immunophenotypic samples, PC1 separated the clusters based on developmental stage (supplemental Figure 5A-B). DEGs were analyzed for each cluster (Figure 4B) and GSEA^24,25 identified a general enrichment of pathways involved in proliferation in fetal clusters, whereas gene sets involved in inflammation and immune system processes were downregulated, in agreement with an earlier study.⁶ The fetal-specific MPP-I cells had only about 130 adult counterparts, making the analysis less robust, but showed upregulation of MYC targets in the fetal MPP-I, as was also detected in the comparison with the fetal HSC and MPP-II clusters (Figure 4C; supplemental Figure 4F).

To identify common upregulated (fetal signature) and downregulated (adult signature) genes, Venn diagrams were used. A gene set of 179 upregulated and 450 downregulated genes was identified in the primitive, lymphoid (Ly), and myeloid (My) groups (Figure 4D; supplemental Figure 5C; supplemental Tables 4-5). The gene signature retrieved was validated on a data set with fetal, induced (i) PSC-derived and adult BM HSPCs,¹⁵ and could separate samples by developmental stage, showing the biological significance of the derived fetal core signature (supplemental Figure 5D). Performing the analysis on pseudobulked samples without the pseudoreplicates generated fewer but highly overlapping genes with the fetal core signature. In contrast, no clear separation of the samples by developmental stage was seen using DEGs identified with a single-cell method (Seurat) (supplemental Figure 5E-H).

Among the genes upregulated in the fetal core signature were LIN28B, IGF2BP1, and IGF2BP3, all part of the LIN28B-let7 axis and known to be involved in self-renewal of fetal HSCs in the murine system³⁸; Delta-like gene 1 (DLK1), a gene shown to be a negative regulator of HSC formation in the mouse embryo;³⁹ and CHD7, an epigenetic remodeler known to interact with RUNX1 and inhibit differentiation⁴⁰ (Figure 4D-E). Furthermore, CD7 was more generally expressed among fetal progenitors, as was also seen at the protein level with the ADT marker (Figure 2B). Among the downregulated genes, the HLA complex of both class I and class II dominated. We also noticed that homeobox (HOX) gene clusters were differentially expressed in fetal and adult HSPCs, where HOXB genes were in general more expressed in fetal and HOXA were more expressed in adult primitive clusters (Figure 4F).

Applying our core to another data set with samples from different developmental stages,⁶ the fetal signature was seen to gradually decrease throughout development, whereas the adult signature increased (Figure 4G). Thus, by directly comparing fetal and adult counterparts from our projection analysis, a fetal core signature with genes specifically upregulated or downregulated in fetal compared with adult counterparts was identified.

By investigating neonatal blood spots and through twin studies, translocations giving rise to ALL in children have been backtracked to birth and an in utero origin.² One question is whether a remnant of the fetal signature can be detected in the pediatric leukemic cells. Looking at the upregulated universal fetal core genes in a publicly available RNA-seq data set of a preleukemic iPSC model of ETV6::RUNX1 (TEL::AML1),¹⁵ fetal genes were to a large extent expressed in the iPS-derived normal (known to recapitulate fetal development) and ETV6::RUNX1-expressing hematopoietic progenitors (Figure 5A).

Next, using a publicly available RNA-seq data set with almost 2000 samples from patients with B-ALL, the fetal core signature was investigated in B-ALLs of known in utero origin.¹⁸ B-ALLs with different driver mutations of interest were selected, pooled, and batch-corrected (KMT2A::AFF1, KMT2A::MLLT1, KMT2A::MLLT3, BCR::ABL1, ETV6::RUNX1, and high or low hyperdiploidy, in total, 645 B-ALL samples; supplemental Figure 6A). The selected samples were pooled, and a PCA plot was generated using the fetal core signature. Intriguingly, the PC1 seemed to separate samples based on age, thus indicating that the fetal core may discriminate between pediatric and adult B-ALL (Figure 5B-C). Investigating the subtypes of B-ALL individually, an association with ETV6::RUNX1 and high hyperdiploid ALL, 2 very frequent cytogenetic abnormalities in children, was seen (P value: 1 × 10^-4 and 3 × 10^-4, respectively),² whereas PC1 variance in low hyperdiploid, BCR::ABL1, and KMT2A::AFF1 did not significantly differ from random variance (Figure 5D; supplemental Figure 6B-C and data not shown).

To investigate if a fetal core program was retained specifically in pediatric B-ALL, a link between our fetal core and age was next examined in leukemia subtypes where both pediatric and adult patients were represented. KMT2A::AFF1 (MLL::AF4) almost exclusively gives rise to B-ALL and is found in many different age groups, including infants, for whom a clear in utero origin has been demonstrated.⁴¹ Furthermore, in infants, few additional mutations are needed in contrast to older children.⁴² The fetal core program was specifically maintained in KMT2A::AFF1 infants compared with older patients. In the other 2 types of leukemia investigated, high hyperdiploid and BCR::ABL1, the latter a subtype seen in about 25% of adult B-ALLs but only 3% to 4% of pediatric cases,⁴³ an association with the downregulated gene set was seen (Figure 5E; supplemental Figure 6D).

To look further into the difference between infant, pediatric, and adult KMT2A::AFF1 patients, a heatmap of DEGs was made (supplemental Figure 6E). HOXA3, HOXA7, and HOXA10 were less expressed in some infants. Looking specifically at the HOX genes expressed in our fetal and adult data sets (Figure 4E), one cluster with mostly infants had lower expression of the different HOXA genes investigated but higher expression of HOXB3 and HOXB4 (Figure 5F). Indeed, a correlation was observed between HOXA expression, prognosis, and age in an earlier study.⁴⁴ 

Thus, in some subtypes of pediatric leukemia associated with an in utero origin, part of a fetal core signature seems to remain in the leukemic cells. The data have the potential to explain differences in pediatric and adult leukemia carrying similar driver mutations, which is potentially important for prognosis and therapy response.

The hematopoietic system is a paradigmatic hierarchical organization of tissues, and careful immunophenotypic isolation and characterization of progenitors have formed the foundation of the hematopoietic hierarchy of today in both mice and human.^45,46 Recent developments in single-cell omics have allowed for further purification and analysis of heterogeneity, and with these findings, the hematopoietic hierarchy can be regarded as a continuum rather than a step-wise hierarchical organization.^47,48 Until recently, these studies have mainly investigated adult hematopoiesis, but more reports looking into human blood development are rapidly emerging.^{6,11-13,49,50} Most of these studies focused on molecular heterogeneity, but recently CITE-seq was used to link the immunophenotype of second-trimester HSCs to engraftment potential.⁵⁰ Within this study, we interrogated FL HSPCs using single-cell RNA-seq, relying on CITE-seq to also capture the immunophenotype along with the transcriptome. By focusing on early primitive hematopoiesis, the whole HSPC population could be interrogated and the relative proportions of all CD34-positive populations compared across developmental time points, thereby also capturing the immunophenotype enabling a unique view of the classically defined progenitors. Although most of the conventional markers of HSPCs (CD90, CD45RA, CD71, and IL7R) showed good correlation with the transcriptionally defined clusters, markers of GMPs and CMPs, that is, CD123 and CD135 (FLT3), were more generally expressed. This was even more apparent when using the ADTs to perform immunophenotypic gating, which showed substantial molecular heterogeneity within GMPs and CMPs. This heterogeneity in defining myeloid progenitors was not specific to fetal development but was also seen in adults.¹⁹ Looking at lymphoid surface markers, CD10 mainly captured the lymphomyeloid cluster Ly-I, whereas our data revealed that IL7R in combination with CD38 and CD45RA could capture the Ly-II cell state with high purity in the embryo. This was the opposite of what was seen in adults, where the population captured by CD10 comprised 70% Ly-II cells vs 43% when using IL7R (data not shown and a previous study¹⁹). Thus, both fetal and adult myeloid immunophenotypically defined progenitors captured molecularly heterogenous populations, whereas HSCs, LMPPs, and MEPs defined fairly pure molecular populations.

An advantage with high-throughput single-cell RNA-seq methods is that there is no need to preselect populations based on immunophenotype. However, difficulties remain when cells from different stages and niches are compared. Earlier studies have merged different stages onto the same analysis, relying on batch correction to enable sample integration.^6,11,49 However, when merging cell populations from different stages into the same analysis, the developmental differences may take over and heterogeneity within the sample itself may be lost. Another problem with integrating data is that it can be difficult to determine if integration is incomplete or if there is a real difference between the samples. Our integrated UMAP is a good example of this, where without the projection approach, the minimal contribution of adult cells in certain clusters could easily be misinterpreted as incomplete integration (supplemental Figure 4). The projection approach used herein allows for signals responsible for lineage determination to be kept.¹⁷ By mapping the query cells onto the reference map, only molecular differences that separate the reference cells will be investigated. From this analysis, cell state differences over development could be captured, and an increase in molecular HSCs was observed with age, whereas the relative fraction of lymphoid progenitors (Ly-I/Ly-II) decreased, in agreement with an earlier study.⁶ In our approach, 2 different reference maps were used: FL CS22 and adult BM. The projection analysis of FL CS22 on the adult BM reference revealed substantial molecular differences between HSCs in FL and adult BM. The molecular FL HSCs mainly projected to adult BM MPP-I and MPP-III, but not to the adult HSC molecular clusters. The discrepancy between fetal and adult HSC states is important to understand as it may facilitate the development of functional HSCs from PSCs for future regenerative medicine.¹ 

Furthermore, the MPP-I cluster was found to be almost exclusively a fetal cell state, suggesting that this population could constitute a fetal–specific progenitor population arising from an HSC-independent wave. The MPP-I cells displayed an erythromyeloid gene signature while maintaining a primitive gene program. No suitable statistical method was found to verify the differences in fetal and adult cell states, as the cells were from 2 donors. That said, using another data set,⁶ our MPP-I was shown to be virtually absent from both pediatric and adult samples, as well as fetal BM, but not from donor-matched FL. For future characterization, studies such as ours should be performed on enriched yolk sac and early FL progenitors. A combination of immunophenotypic markers presented in our study (supplemental Figure 2B) may allow for prospective isolation and further characterization of the progenitor population and reveal its origin. Such studies may also provide markers and molecular signatures that can separate early from definitive hematopoietic waves.

By using the projection analysis, fetal and adult counterparts could be transcriptionally compared, and a core set of genes linked to fetal or adult identity could be identified. Multiple methods exist to perform DEG analysis on single-cell data sets, and in general, pseudobulking approaches are found to be more reliable.^51,52 Pseudoreplicates were made to account for intradonor variability as well as to increase statistical power in the analysis as donors were limited. The derived core was validated and compared by other methods and shown to be highly biologically relevant, and our data highlights substantial ontogeny-dependent transcriptional differences.

Childhood B-ALL has in many cases been shown to have a fetal origin,² and the leukemic cells could potentially retain expression of fetal-specific genes or lack expression of adult-specific genes.⁵³ Thus, our fetal core signature was used to investigate childhood leukemia from publicly available data of B-ALL samples with different cytogenetic abnormalities.¹⁸ The fetal core signature correlated with the age of the patient in some leukemia subtypes, indicating that a fetal gene expression program may remain in these patients. For KMT2A::AFF1, it was clear that HOXA genes were downregulated in some infants compared with those in adults, in agreement with an earlier study where a correlation was observed between lower expression of HOXA genes, poor prognosis, and young age.⁴⁴ Of note, in our fetal-adult data set HOXA genes were less expressed in fetal HSCs compared with adult (Figure 4F).

Our study links immunophenotypes with transcriptomes at the single-cell level, providing a unique map of human fetal blood development. Comparison to adult hematopoiesis gives insight into how the hematopoietic progenitor compartment changes with development. A fetal core signature with universally upregulated or downregulated genes depending on developmental stage could be identified, a gene set that could separate specific subtypes of B-ALL samples based on age. The importance of the fetal core genes in the formation of the preleukemic clone in utero and in leukemia progression and development in general remains to be explored.

The authors acknowledge the personnel of the Lund Stem Cell Center FACS core facility, Anna Fossum and Zhi Ma, and Center for Translational Genomics (CTG) at Lund University for their assistance with the CITE-seq experiments, especially Linda Geironson Ulfsson and Eva Erlandsson, and Ariana Calderón for help with the graphical abstract. The authors acknowledge support from the StemTherapy Program, the Crafoord Foundation (C.B.), the Swedish Childhood Cancer Foundation (C.B. and G.K.), the Ragnar Söderberg Foundation (C.B. and G.K.), the Knut and Alice Wallenberg Foundation (G.K.), the Swedish Research Council (C.B. and G.K.), and the Swedish Cancer Society (G.K.).

Contribution: M.N.E.S., G.K., and C.B. conceptualized the study; M.N.E.S. performed the methodology and conducted the formal analysis; M.N.E.S., P.D., and R.O. curated the data; M.N.E.S., S.P., and R.O. conducted the investigation; M.N.E.S. and C.B. visualized the data and wrote the original draft of the manuscript; M.N.E.S., S.P., G.K., and C.B. reviewed and edited the manuscript; C.B. and G.K. acquired funds; and S.S., G.K., and C.B. supervised the study.

Conflict-of-interest disclosure: P.D. and G.K. are board members and have equity in Nygen Analytics AB, all unrelated to this work. The remaining authors declare no competing financial interests.

Correspondence: Göran Karlsson, Division of Molecular Hematology, Lund Stem Cell Center, Lund University, BMC B12, 221 84 Lund, Sweden; e-mail: goran.karlsson@med.lu.se; and Charlotta Böiers, Division of Molecular Hematology, Lund Stem Cell Center, Lund University, BMC B12, 221 84 Lund, Sweden; e-mail: charlotta.boiers@med.lu.se.