Robust temporal map of human in vitro myelopoiesis using single-cell genomics

Abstract

Myeloid cells are central to homeostasis and immunity. Characterising in vitro myelopoiesis protocols is imperative for their use in research, immunotherapies, and understanding human myelopoiesis. Here, we generate a >470K cells molecular map of human induced pluripotent stem cells (iPSC) differentiation into macrophages. Integration with in vivo single-cell atlases shows in vitro differentiation recapitulates features of yolk sac hematopoiesis, before definitive hematopoietic stem cells (HSC) emerge. The diversity of myeloid cells generated, including mast cells and monocytes, suggests that HSC-independent hematopoiesis can produce multiple myeloid lineages. We uncover poorly described myeloid progenitors and conservation between in vivo and in vitro regulatory programs. Additionally, we develop a protocol to produce iPSC-derived dendritic cells (DC) resembling cDC2. Using CRISPR/Cas9 knock-outs, we validate the effects of key transcription factors in macrophage and DC ontogeny. This roadmap of myeloid differentiation is an important resource for investigating human fetal hematopoiesis and new therapeutic opportunities.

Fig. 1. iPSC macrophage differentiation produces a range of foetal myeloid and stromal cells.

a Schematic illustration of the in vitro differentiation protocol from iPSC to macrophages highlighting the time points sampled for scRNAseq and scATACseq profiling. Full culture well/plate was collected at each time point and at day31 the non-adherent fraction of the culture was also processed independently (day31 Non-adh). The protocol was repeated twice to generate the Discovery and Validation data sets. b Computational workflow diagram for cell-type annotation. Briefly, LR models were used to annotate the Discovery scRNAseq data set based on publicly available in vivo data sets of human gastrulation (Gas), yolk sac (YS), foetal liver (including skin and kidney) (FLi), foetal thymus (FTh) and placenta (Pla). Then, annotations were transferred to the scRNAseq and scATACseq Validation data sets. c UMAP projections of the Discovery scRNAseq data labelled by cell type. In vivo data sets supporting the cell type annotation and the area under the curve (AUC) for the best performing LR model are listed. *In vivo data set LR model of the AUC shown. (right) UMAP projections of the Discovery data set labelled by time point. d Dot plot showing canonical markers expression for each of the cell types. Colours depict the mean gene expression and dot size the percentage of expressing cells. e UMAP projections of the scRNAseq Validation data set (n = 62,000) labelled by cell type. f UMAP projections of the scATACseq Validation data set (n = 71,000) labelled by cell type. g Heatmap showing the mean logistic regression models’ predicted probabilities of the YS hematopoietic cell types for each of the cell types in the Discovery scRNAseq data set. h UMAP projections of the scRNAseq Discovery data coloured by the LR models’ predicted probabilities of the YS hematopoietic cell types. iPSC induced pluripotent stem cell, EB embryoid body, Mac macrophage, LR logistic regression, UMAP uniform manifold approximation and projection, AUC area under the curve, ATAC assay for transposase-accessible chromatin, YS yolk sac, YSMP yolk sac myeloid-biased progenitors. Source data are provided in the Source data file.

Fig. 2. Cell population dynamics.

a Diagram illustrating the dynamic emergence of the different cell types over the course of the in vitro differentiation protocol (Discovery data set), weaker links are shown by discontinuous lines. b Schematic representation of the computational workflow used to compare transcription factor (TF) dynamics in vivo and in vitro. Briefly, TF activities were computed at branching points along the in vitro differentiation trajectory (Discovery data set) and were compared to TF activities in matched cell types in the in vivo human yolk sac, gastrulation and foetal liver data sets. c RNA velocity analysis and PAGA graph abstraction of the cells present at day3 (Embryoid body (EB) formation) of the differentiation protocol (Discovery data set) showing the developmental relationships between cell types. d Transcription factor activities computed with DoRothEA for the identified cell types present at day3 of the in vitro differentiation protocol and matched cell types in the in vivo gastrulation data set, relevant TFs are shown in bold, asterisks highlight significantly different activity vs its progenitor, Bonferroni adjusted p < 0.05. e RNA velocity analysis and PAGA graph abstraction of the cells present at day21 (EB myeloid differentiation) of the differentiation protocol (Discovery data set) showing the developmental relationships between the cell types. f Transcription factor activities computed with DoRothEA for the identified cell types present at day21 of the in vitro differentiation protocol and matched cell types in the in vivo yolk sac data set, relevant TF shown in bold, asterisks highlight significantly different activity vs its progenitor, Bonferroni adjusted p < 0.05. g Same analysis as f with matched cell types in the in vivo foetal liver, skin and kidney data set. h Violin plots showing the number of accessible peaks per cell type in the scATACseq Validation data set. Each plot shows a distinct lineage, asterisks highlight two-tailed t test Bonferroni adjusted p < 0.05. iPSC induced pluripotent stem cells, ATAC assay for transposase-accessible chromatin. Source data are provided in the Source data file.

Fig. 3. Evaluation of the macrophage phase.

a Schematic illustration of the in vitro differentiation protocol and cell-type annotation analysis, steps before the macrophage differentiation phase are hidden. Additional experimental conditions (alternative cytokines, top, and media experiment, bottom) are highlighted in red, standard protocol conditions are in black. b (right) UMAP projections of the macrophage phase labelled by cell type and time points. All experiments are pooled. (right) UMAP projections highlighting the cells included in each experiment. c Stacked area plot of the cell type percentages in each time point. Only samples from the time points experiment (M-CSF) were included. d (left) UMAP projection highlighting macrophages from the time points experiment (M-CSF only) and coloured by time point. (right) Heatmap of the transcription factor activity scores calculated using DoRothEA across time points, relevant TFs are in bold, asterisks highlight significantly different activity vs day31, Bonferroni adjusted p < 0.05. e TF motif enrichment values in macrophage ATAC open peaks at day31 + 7 vs day31 plotted against TF transcriptional activity score at day31 + 1 (left) or day31 + 7 (right). Pearson correlation’s r and exact p values are shown. f (left) UMAP projection highlighting macrophages from the cytokines experiment and coloured by time point and cytokine cocktail used. (right) Heatmap of the TF activity scores across time points and cytokines, relevant TFs are in bold, asterisks highlight significantly different activity vs day31, Bonferroni adjusted p < 0.05. iPSC induced pluripotent stem cells, TF transcription factor, ATAC assay for transposase-accessible chromatin. Source data are provided in the Source data file.

Fig. 4. LPS stimulation of distinct subtypes of macrophages.

a Schematic illustration of the in vitro differentiation protocol and cell-type annotation analysis, steps before the macrophage differentiation phase are hidden. Additional experimental conditions (LPS stimulation at distinct time points) are highlighted in red. b UMAP projections of scRNAseq analysis of LPS stimulated samples and matched controls for 4 populations of macrophages. c (top) Dot plot of genes overexpressed >3 log fold-change in either of the 4 LPS stimulated samples vs their control, (bottom) Dot plot of genes significantly overexpressed in 1 of the 4 populations analysed. d Cytokines protein expression levels in supernatants after LPS stimulation (n = 2) and controls (n = 1) for the 4 experimental conditions analysed. Media was collected from samples processed for scRNAseq (shown in b, c). Source data are provided in the Source data file.

Fig. 5. Modification of differentiation cytokines produces dendritic cells.

a (left) Schematic illustration of the in vitro differentiation protocol from iPSC to dendritic cells highlighting the time points when samples were collected for scRNAseq profiling. (right) Computational workflow diagram for cell-type annotation. b UMAP projections of the scRNAseq data labelled by cell type (left) and time point (right). In vivo data sets supporting the cell type annotation and the area under the curve (AUC) for the best performing LR model are listed. *In vivo data set LR model of the AUC shown. c Dot plot of canonical marker genes expression for each cell type. d Transcription factor (TF) activities computed with DoRothEA for the identified cell types present at day21 of the in vitro differentiation protocol and matched cell types in the in vivo yolk sac data set and foetal liver, skin and kidney data set. Relevant TFs are in bold. e TF activities computed with DoRothEA for cDC2 identified at the last 4 time points of the differentiation protocol (step 3). Relevant TFs are in bold, asterisks highlight significantly different activity vs day31, Bonferroni adjusted p < 0.05. f Stacked area plot showing the proportions of the major cell types from day31 + 1 to day31 + 7. g Flow cytometry histograms showing the protein levels of cDC2 marker genes and CD14 as a negative marker in non-adherent cells at the end of the DC differentiation phase (day31 + 7), matched unstained controls are shown in grey. h Flow cytometry histograms for BODIPY™ FL DQ-ovalbumin processing by non-adherent cells at the end of the DC differentiation phase (day31 + 7) incubated for 15, 45 and 60 min at 37 °C. In grey are matched samples kept at 4 °C as a negative control. i Flow cytometry histograms for a T cell activation assay, peaks with lower CFSE signal than unstimulated T cell negative control (in grey) correspond to proliferative/activated T cells by the presence of anti-CD28, iPSC-derived macrophages or iPSC-derived DCs. Plots shown are representative of two donors and two independent experiments. Source data are provided in the Source data file.

Fig. 6. Effect on macrophage differentiation of ICAM1, LSP1, PRKCB and ZEB2 KO.

a (left) Schematic illustration of the in vitro differentiation protocols from iPSC to macrophages (MAC, top) or dendritic cells (DC, bottom) used to evaluate the effects of ICAM1, LSP1, PRKCB or ZEB2 knock-outs (KO). Samples were collected at day0 and day31 of the protocols and profiled with scRNAseq. (right) Computational workflow diagram for cell type annotation. Briefly, cell type annotations were transferred from scRNAseq data of the macrophages (Discovery data set) and DC protocols described in the previous sections. b UMAP projections of scRNAseq data from both KO protocols labelled by time point. c UMAP projections of scRNAseq data generated from the iPSC-to-macrophages KO protocol (one UMAP per KO and wild type) coloured by cell type. d UMAP projections of scRNAseq data generated from the iPSC-to-DC KO protocol (one UMAP per KO plus wild type) coloured by cell type. e Dot plot showing the average expression of intermediate monocyte–associated genes in the monocytes produced by each KO and the wild type in the iPSC-to-DC protocol. f Dot plot showing the average expression of genes associated to myeloid-derived suppressor cells in the monocytes produced by each KO and the wild type in the iPSC-to-DC protocol. g (left) Dot plot showing the average expression of M2-associated genes in the macrophages produced by each KO and the wild type in the iPSC-to-macrophages protocol. (right) Transcription factor activities computed with DoRothEA for macrophages produced by each KO and the wild type.