ARTICLE https://doi.org/10.1038/s41467-022-30557-4 OPEN Robust temporal map of human in vitro myelopoiesis using single-cell genomics Clara Alsinet 1,2 ✉, Maria Nascimento Primo 1,2, Valentina Lorenzi 1, Erica Bello 1,2, Iva Kelava1, Carla P. Jones 1, Roser Vilarrasa-Blasi 1, Carmen Sancho-Serra1, Andrew J. Knights 1, Jong-Eun Park 3, Beata S. Wyspianska2,4, Gosia Trynka 1,2, David F. Tough4, Andrew Bassett 1,2, Daniel J. Gaffney 1 ✉, Damiana Alvarez-Errico 5 ✉ & Roser Vento-Tormo 1 ✉ 1234567890():,; 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 defi- nitive 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. 1 Wellcome Sanger Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SA, UK. 2 Open Targets, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SA, UK. 3 Graduate School of Medical Science and Engineering, Korea Advanced Institute of Science and Technology (KAIST), Daejeon 34141, Korea. 4 Immunology Research Unit, Medicines Research Centre, GlaxoSmithKline, Stevenage SG1 2NY, UK. 5 Josep Carreras Leukaemia Research Institute (IJC), Badalona, 08916 Barcelona, Catalonia, Spain. ✉email: ca6@sanger.ac.uk;dg13@sanger.ac.uk;dalvarez@carrerasresearch.org; rv4@sanger.ac.uk NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 M acrophages perform a variety of functions, from tissue The differentiation protocol consists of 3 steps: (i) spin-embryoid homoeostasis to immune surveillance and from the body (EB) formation from day1 to day4, (ii) EB myeloid differ- response to infection to the resolution of entiation from day5 onwards (the latest sample used in this study inflammation1–4. They originate during both development and is from day31), and (iii) macrophage differentiation using non- adulthood and acquire tissue-specific functions5. Despite the adherent cells from EB myeloid differentiation phase and lasting commonalities within mammals, there are important differences for 7 days, in this study we used cells from day31 (day31 to between humans and rodent models6. Establishing and char- day31 + 7). To characterise the robustness of our results, we acterising the current human in vitro models is essential to fully generated two independent scRNAseq data sets. The first data set exploit their research and therapeutic potential7. (referred to as Discovery data set) included scRNAseq data from 3 During development, myeloid cells originate from at least two donors at 20 timepoints (Fig. 1a, b). The second data set (here- waves of progenitors: a first wave involving myeloid-biased pro- after, Validation data set) included scRNAseq and scATACseq genitors from the yolk sac (yolk-sac myeloid progenitors, YSMP) data from 6 donors at 7 and 6 time points, respectively (Fig. 1a, b). and a second wave through definitive hematopoietic stem cells The three donors from the Discovery data set were also used in the (HSC)5,8. YSMP are thought to appear during the first 2 weeks of Validation data set, thus generating biological replicates. development in humans and are responsible for producing pri- After quality control, the Discovery data set contained a total of mordial blood5. HSC are not generated until 3–4 post- 135,000 cells (Fig. 1c and Supplementary Fig. 1). To annotate the conceptional weeks (PCW) in the gonad-aorta-mesonephros. cell types in an unbiased manner, we built logistic regression (LR) HSC and myeloid progenitors derived from YSMP colonise the models trained on publicly available single-cell transcriptomics liver, making this foetal organ the main site of hematopoiesis data sets and projected the data into our in vitro data set (Fig. 1b). until mid-pregnancy9. Later, HSC are restricted to the bone We used multiple human developmental data sets to train our marrow, the only hematopoietic site during adulthood10. In mice, models6,20–23 (Supplementary Data 2), including a full gastrula YSMPs generate a wide range of myeloid cells, including mono- and the main foetal hematopoietic organs: yolk sac, liver and cytes and neutrophils, and are thought to be the main myeloid thymus. Cell type labels were assigned based on the mean LR precursors during development11. In humans, we have limited prediction probability of each cell cluster (Fig. 1c, Supplementary knowledge about the progression and regulatory mechanisms Fig. 2 and Supplementary Data 3). Marker gene expression defining YSMP. analysis further supported the LR cell-type annotation (Fig. 1d). In vitro models of macrophage differentiation hold promise to Cell type label transfer24 from the Discovery data set into the not only answer these biological questions but also to become Validation and scATACseq data sets confirmed the presence of therapeutic tools, particularly immunotherapies, and for high- the main cell populations in all 3 data sets (Fig. 1e, f, throughput screening, including drug testing. Macrophages Supplementary Fig. 1 and Supplementary Data 4). derived from human induced pluripotent stem cells (iPSC) show The majority of cells at the initial EB formation stage (day1–4) tissue-resident phenotypes12 and are an attractive alternative to matched gastrulation cell populations (Fig. 1c and Supplementary adult monocyte-derived macrophage cultures13,14. A current Fig. 3). We found primitive streak-like cells, emergent and protocol, developed by vanWilgenburg et al., is a straightforward, advanced mesoderm, and the initial appearance of hemogenic feeder-free process consisting of 3 steps using constant con- endothelium. Despite using cytokines that induce hematopoietic centrations of 1 to 3 cytokines15. It provides long-term, scalable mesoderm (Fig. 1a), we also observed populations related to other production of macrophage precursors without fluorescence- germ layers (i.e. neural crest and endoderm, Fig. 1c). activated cell sorting (FACS) since the cells of interest con- During EB myeloid differentiation (day5–31), the myeloid and tinuously expand and detach from the culture16. Despite this stromal cell compartments emerged (Fig. 1c and Supplementary being an established in vitro macrophage model, the exact Fig. 3). The myeloid populations included a wide range of cell intermediate populations produced are unclear12. This restricts its types, such as erythrocytes, megakaryocytes, mast cells, neutro- applications and limits our true understanding of the cells phil myeloid progenitors (NMP), monocyte DC precursors obtained. Thus, a thorough analysis of the cell identities and (MDP), monocytes and macrophages. Being so unexpected, we dynamics emerging during in vitro differentiation is imperative to validated the presence of the small population of mast cells by establish their likely in vivo counterparts and fully exploit this FACS (Supplementary Fig. 4). Most in vivo counterparts for these technology17–19. cell types were found in the foetal liver and thymus. However, we Here, we profiled the single-cell transcriptome and open did not find any cluster in the Discovery data set that chromatin of >400k and >70k cells, respectively, during corresponded to the in vivo HSC found in the human developing iPSC–myeloid differentiation with the vanWilgenburg protocol15. liver (Supplementary Fig. 2B). Instead, there was a distinct cluster We uncover a wide range of cell states, their ontogeny and of myeloid progenitors (MP) expressing CD34, low levels of underlying transcription factor (TF) networks, that accurately SPINK2 and PTPRC, but not HOXA genes, which are required to map foetal myelopoiesis in the YS. We demonstrate the versatility generate definitive HSC25,26 (Fig. 1d and Supplementary Fig. 5A). of the current in vitro protocol by modifying the media used and MPs had a high prediction probability for the YSMP-trained using CRISPR/Cas9 edited iPSCs. Altogether, we demonstrate model generated with the embryonic YS data set (Fig. 1g, h and that macrophage differentiation from iPSC is a robust system to Supplementary Fig. 2C), suggesting in vitro myelopoiesis study the early stages of human myelopoiesis, and that macro- recapitulates YS differentiation. In vivo YSMP and macrophages phages obtained are able to acquire definitive tissue-resident LR models captured more than one cell type within the in vitro identities. Our data sets can be visualised and downloaded from data set (Fig. 1g, h). To explore this further, we performed the www.HiPImmuneatlas.org. opposite exercise: we trained models on our in vitro cell types and projected them onto the in vivo YS data set. As expected a subset of cells within the published YSMP cluster showed a high Results prediction probability with the in vitro NMP-trained model and a iPSC-derived cells have human yolk sac myelopoiesis features. subset of in vivo macrophages were captured by the MDP-trained We profiled the full differentiation of iPSC into macrophages model (Supplementary Fig. 5B, C). Using the LR results, we from 6 individuals using single-cell RNA and ATAC sequencing annotated the NMP and MDP cell types within the embryonic YS (scRNAseq and scATACseq) (Fig. 1a and Supplementary Data 1). data set (Supplementary Fig. 5B). 2 NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 ARTICLE To quantitatively characterise the macrophage foetal-like confounders of adult+foetal data set integration. The Hofbauer profile observed, we projected adult and foetal macrophages data cells LR model had a higher mean prediction probability for from the human decidual–placental interface into our data set. iPSC-derived macrophages than any of the adult macrophage This unique tissue setting includes both adult/maternal subtypes identified in the placenta (Supplementary Fig. 5D, E). As monocyte-derived macrophages and foetal/placental YS-derived an exception, day31 + 1 macrophages presented a higher score macrophages (Hofbauer cells)20, thus avoiding the technical with adult macrophages LR models (Supplementary Fig. 5D, E). NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 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)23, yolk sac (YS)6, foetal liver (including skin and kidney) (FLi)22, foetal thymus (FTh)21 and placenta (Pla)20. 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 types6 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 types6. 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. This was likely because day31 + 1 macrophages showed an Later, hemogenic endothelium activated TFs relevant for activated state, upregulating inflammatory cytokines such as hematopoiesis including RUNX1, SPI1, RBPJ, MEF2A and CXCL8, CCL7 or IL1B, (Supplementary Fig. 5F) also found on MEF2C. GATA1 was also activated in this transition but it monocyte-derived macrophages in the decidua (adult/maternal showed the highest activation levels in erythrocytes (Fig. 2d). tissue). Next, we used a hepatocellular carcinoma data set27 and In vitro, myelopoiesis started very early in the EB formation projected several tumour-associated macrophages (TAMs) LR stage, but the wide range of myeloid cell types appeared almost models on our iPSC-derived data set (Supplementary Fig. 5G, H). simultaneously starting at day14, and they all endured until day31 We found the foetal-like FOLR2+ TAMs model showed a higher (Fig. 2a and Supplementary Fig. 6A). The first MP appeared at prediction score among end-stage macrophages (day31 + 7) day9–11 followed by erythrocytes on day11–14, and full while the SPP1+ TAMs LR model markedly captured macro- myelopoiesis was achieved on day16–18. In addition to myeloid phages on the activated state (day31 + 1, Supplementary Fig. 5H). cells, advanced mesoderm also differentiated into an intermediate Overall, this indicated that macrophages produced in the iPSC stage of early fibroblasts (day7), giving rise to fibroblasts by day9. protocol have a strong foetal phenotype, and this could be Trajectory analysis on the sample at day21 reconstructed all relevant for their application as in vitro TAM models. myelopoiesis differentiation steps. Hemogenic endothelial cells, derived from the mesoderm, differentiated into MP, which gave Trajectory analysis and underlying regulatory programmes. rise to both megakaryocytes and NMP (Fig. 2e and Supplemen- Myelopoiesis is shaped by transcriptional programmes, including tary Fig. 6A). NMP gave rise to MDP, which differentiated into TFs, epigenetic regulators and post-transcriptional either monocytes or macrophages (Fig. 2e). The MDP to mechanisms28. We set out to reconstruct the main develop- macrophage trajectory was consistently observed from day21 to mental pathways underlying in vitro myelopoiesis and the tran- day31 (Supplementary Fig. 6B). The differentiation pathway of scriptional networks mediating them. The high-throughput macrophages through MP and bypassing the monocytes is single-cell approach used, the high density of time points col- consistent with the first waves of myelopoiesis emerging in the lected and trajectory analysis using scVelo29, allowed us to YS5. reconstruct the differentiation paths giving rise to the wide range Throughout all stages of myelopoiesis, we found high similarity of cell types observed (Fig. 2a and Supplementary Fig. 6A). In between the regulatory programmes activated in vivo (YS and parallel, we compared the underlying regulatory programmes foetal liver) and in vitro (Fig. 2f, g and Supplementary Data 5). mediating such transitions in vivo and in vitro. To this end, we The transition from hemogenic endothelium to MP was measured TF activities by looking at the expression of consensus characterised by the activation of multiple TFs including RUNX1, TF targets30 (Fig. 2b). SPI1 and GATA1 (Fig. 2f, g). The MP to NMP transition showed In step 1 of in vitro differentiation, iPSC differentiated into the further activation of SPI1 and CEBPA. On the contrary, a large primitive streak, which subsequently gave rise to either endoderm number of endothelial TFs, such as SOX2 and the ETV family, or emergent and advanced mesoderm (Fig. 2c). The split of pluripotency factors, such as POU5F1 and NANOG, or lymphoid primitive streak cells observed (Fig. 2c) suggested there is lineage-promoting factors, such as MEF2C, were inactivated32. heterogeneity within this cell type. Later, advanced mesoderm MEF2C is a TF characteristic of definitive HSC that drives differentiated into hemogenic endothelium, which was the lymphoid fate choice32. The low MEF2C activity in in vitro MPs, precursor of myeloid cells (Fig. 2c). Primitive streak and and in vivo YSMPs, further supported the HSC-independent mesoderm were transient populations that disappeared by features of this system (Fig. 2f). day16, while endoderm and hemogenic endothelium were iPSC-derived differentiation towards monocytes and macrophages detectable until at least day31, the latest time point analysed also shared transcriptional programmes with its in vivo YS and foetal (Fig. 2a and Supplementary Fig. 6A). The transition from liver counterparts. Among the few TFs specifically activated in MDP mesoderm to hemogenic endothelium has also been reported from NMP was RFX5, which regulates MHC-II transcription and is during the gastrulation period31. TF activity analysis showed high responsible for a rare hereditary immunodeficiency33. We also conservation of the TF modules in these earlier stages of observed activation of TFs controlling inflammatory programmes in development (Fig. 2d and Supplementary Data 5). Pluripotency monocytes and macrophages, such as JUN, RELA and NFKB1. In TFs activity (POU5F1, NANOG and SOX2) decreased upon macrophages, we found key similarities related to their myeloid differentiation into mesoderm and endoderm. As expected, in identity such as MAF34 and CEBPB35 and potentially underlying both in vivo and in vitro settings, mesoderm activated SMAD4, tissue-specific programmes, such as the alveolar macrophage HOXA9 or SRF while endoderm activated FOXA2 and HNF4A. programme characterised by PPARG36,37 (Fig. 2f, g). 4 NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 ARTICLE DAY 31 b TF activity comparison analysis a DAY 0 TIME in vitro in vivo NCC Neural Crest Cells iPSCs Human End iPSC Endoderm gastrulation MEc STROMAL MesodermEctoderm Endoderm dataset PrS EmM AdM EFib Fib Fibroblasts Hemogenic Endothelium Erythrocytes Human Hemogenic Endothelium HEn yolk sac Human Megakaryocytes Ery Erythrocytes dataset fetal liver Mast cells MDP NMP Mega Myeloid Progenitors dataset MP Megakaryocytes Mast MYELOID Mast Cells Neutrophil Myeloid Prog. Macrophages Monocyte DC Precursors DoRothEA analysis DoRothEA analysis Mo Monocytes along each trajectory’s cell types in matching cell types Mac Macrophages top variable TFs top variable TFs c PAGA velocity graph Day3 EB formation d Transcription factor activity Gastrulation vs iPSC-derived cells iPSC Trajectories Trajectories Trajectories * * * CEBPB * * E2F2 1.5 * * GATA1 iPSC Induced Pluripotent Stem Cells * * * GATA3 * E2F3 0.5 1 * * ESR2 2 PrS Primitive Streak * * * SP3 * MYCN * * MEF2C EmM End Endoderm * EPAS1 0.5 * * * SMAD3 * * MEF2A 1 EmM Emeregent Mesoderm 0 * * NR5A1 * * * ESR1 0 * * RBPJ AdM AdM Advanced Mesoderm * * * SMAD4 * SOX2 * * TCF4 0 PrS HEn Hemogenic Endothelium * CEBPD −0.5 * * KLF9 * * * HOXA9 −0.5 * POU5F1 * TCF12 −1 * * GATA6 * * ZNF263 HEn * * * RUNX2 * ETV4 End * * MNT JUN * * * * ETV4 * * * * ZBTB7A * * ERG * * * NR2F2 * * NANOG * * IKZF1 e PAGA velocity graph Day21 EB formation * * ATF2 * * E2F5 * * LYL1 * * * * STAT3 * * PAX6 * * SPI1 HEn Hemogenic Endothelium * * * * SOX10 * * TFDP1 * * RUNX1 Mac HEn MP Myeloid Progenitors * * * * AHR * * USF2 * * MAF Mast Mast Cells * * * * SRF * * ESRRA * * TBX21 MDP Mega Megakaryocytes * * * NANOG * * SOX13 * * ATF3 Mo * * NR5A1 * * SP3 Fib NMP Neutrophil Myeloid Progenitors * * CREB3L1 * * TFAP2C MDP Monocyte DC Precursors * * * POU5F1 * AR SOX2 * * GATA6 Mo Monocytes * * * MAX NMP CEBPD * * KDM5B * * * * CREM Mac Macrophages * * KLF4 * * THAP1 KLF6 * * Fib Fibroblasts E2F4 * * FOXA2 TEAD1 Mast * * * * KLF5 * MP * * MYC * AHR Mega * * HNF1A * ESR1 * * * HIF1A * * NR1H3 * * * SP1 * * * * SOX9 HNF4A * * * HHEX * * HNF4G * * GATA6 f * * * SOX9 * * PPARA TEAD4 Transcription factor activity Yolk sac vs iPSC-derived cells n_ AS En _iP S y_ S Er SC SC SC S_ AS SC d_ AS SC Ad iPSC M_ AS Em iPSC S_ S SC M_ AS M A Er _GA Pr _GA HE h_G Ad _G iP iP Pr _G En d_G iP Em _G iP Ad _G iP Trajectories M y S Trajectories M t S M do Ad En HNF1A Pr Pr 1 * * * * * * HIF1A * * * NANOG 1 * * CEBPG g * * * SMAD5 0.5 * * * PPARA 0.5 Transcription factor activity Fetal Liver (+Kidney+Skin) vs iPSC-derived cells * * * HOXB13 0 * * * GATA6 0 THAP11 Trajectories * * * PAX6 * * * * −0.5 1 Trajectories Trajectories * * * * ESR2 * * CEBPD * * * * TCF4 −0.5 STAT1 * * * * −1 * * ETV4 * * * * * TFAP2C 1 * * * * NANOG SP3 * * * * BHLHE40 * * * * TEAD1 −1 * * * −1.5 * * * SOX2 * * * * * BHLHE40 * * * ESR1 * ESR1 * * ZEB2 1 * * * * FOXA2 * * * * ETV4 * * * HNF4G * PAX6 0 * * * ZNF263 * * ATF6 * * * * * * * ETV1 * * * JUN * * * SOX9 * * * * THAP11 E2F1 0 * * * * * * * * * TP53 REL * * * * * * NANOG * * * ZEB1 0 * * * * * * NFYB * * * * POU5F1 HIF1A * * * * SP3 * * * * * * * TEAD1 * * * HOXA9 * * * * * MYC * * MEF2C * * * TCF4 −1 * * * * * E2F2 −1 * * * * LHX2 * * * NR5A1 * * * LYL1 * * TEAD4 −1 * * * * * * TFDP1 * * * * MEIS2 RUNX1 * * * * FLI1 * * * * * * MEIS1 * * * * * * NFYA * * * * SOX2 JUNB * * * STAT5B * * * * * * POU5F1 * * * * * RFX5 * * * MEF2C MAF TBP * * * PBX3 * * * JUNB * * * * * * * * * * * * NR5A1 * * * * GATA6 * * MEF2A * * * * GATA1 * * * * * CEBPB * * * * MAFK SPI1 * * * * ATF6 * * * * * * SP4 * * * * PKNOX1 * * * * * * THAP11 * * * * PPARG * * * HNF4G * * * KDM5B * * * * * * LYL1 * * * CEBPA FOS * * * GATA1 * * * * * * * ETV1 * * * MEIS2 * * * * * * MAF * * * * STAT3 * * * MEF2C * * * * * SP3 * * * STAT5B * * * NFE2L1 * * * * NFKB1 * * * FLI1 * * * * * * JUNB * * * MYC * * * FOXK2 * * * E2F2 * * * * RELA * * NR2F2 * * * * * * STAT1 * * * * NFYA NR5A1 * * * * MEF2C * * E2F3 * * * * * * * * * TP53 * * * * RFX5 * RFX5 * * * * MYC * ETV4 * * * E2F4 TFAP2C * * * * * * * * NR2C2 * * * * * * FOXP1 * * * MEF2B * * * * * * SPI1 * * * * REST * * * GATA1 * * * KLF9 * * * * * FOXL2 * * * * RUNX1 SREBF2 * * BATF * * * * * * SP1 * * * MEF2A * * * * * * NFKB1 * * * HHEX * * * LYL1 CREB3 RELA * * * * POU2F2 MYC * * * * * * * * * * * * * * * * ETS1 KLF1 * * * * * * ATF2 * * * SPI1 E2F2 * * * ATF3 * * * * * * NR2C2 * * * * * * ATF4 * * * E2F4 * * * * USF1 * * * STAT4 * * * * * * CEBPA * * * * * * JUN * * * NR2C2 HNF4A * * * * JUNB * * * * * * * * ELK1 n_ YS SC YS PSC _iP YS SC P_ YS MYB * * * * BACH2 SC SC C * * SRF P_ YS P_ _YS c_ S * * * * * * HE th_ MP P_ iP iPS NM MP_ Ma ac_Y CEBPB * * * KLF5 ESR2 iP iP * * * i NMMP_ * * * * * * M MD DP do * * * FOXA1 N * * * * * * PPARG En M st_ FLi M SC ga FLi M C N Me SC _iP i MP Li RUNX1 MP _FL * * STAT3 iPS * * * * * * P ME PP_F Ma ast_ Me ga_ _i * * * SPI1 _iP i P_ FLi Ma C P_ Li c_ Li C C Mo _FL C M S iPS iPS iPS SC NM P_F Ma c_F MP P_F i C N M SC _iP Li P_ FLi M FL MD DP_ C/ C/ iPS Mo NM iP HS ME PP_ NM P_ M HS En_ M H h iPSC to other myeloid lineages iPSC to Endoderm iPSC to Fibroblasts iPSC to Macrophages * # accessible peaks per cell * * * * * * * * * * 30000 * * 20000 * 20000 * * < 20000 < 20000 < < 10000 10000 10000 10000 0 0 0 0 eS s ce so k es r m do rm M ythr ium yo s as s ls s s do k rm H van Me eak tro ge Me erm do rm Pr lium pr itor on or ag rop es ay s +7 C a e M yte ve Cs k Ea es rm br rm s s er iti SC iti SC En rea _D ge el tre ak cyt ea Fi last st M urs og dM de En ode Em rim iPS t de ilM En de 31 tC ro Ma ocy n l e e S tr es ha Er the c la d on oid the Em Prim iP iP St eD oge nt Str ed od r ly od o Ad en eS iP ec ce so so ob ob ve H van Me ar er itiv s v br e c C M M iti Fi t im t d Ad en im eg ic ph nic Pr P g ge M yel Pr en g nc yt ph oc er va Em o em Ad N em ac M * bonferroni adj. p-value < 0.05 < less than 20 cells eu Finally, we evaluated the chromatin accessibility dynamics of progenitors (‘Emergent Mesoderm’ n = 9759, ‘Advanced Mesoderm’ iPSC–macrophage differentiation using scATACseq. The number of n = 8158). These differences were not associated with the number of accessible cell peaks decreased alongside the trajectories identified, expressed genes (‘Hemogenic endothelium’ n = 2514, ‘Emergent suggesting a more restrictive chromatin landscape as cells differ- Mesoderm’ n = 3463, ‘Advanced Mesoderm’ n = 2725). Another entiate. As an exception, hemogenic endothelium had a higher exception was the macrophages after the macrophage differentiation median of accessible peaks per cell (n = 11,074) than its mesoderm phase (‘Macrophages_Day31plus7’), which had more accessibility NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 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 sac6, gastrulation23 and foetal liver22 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 set23, 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 set6, 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 set22. 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. peaks than the macrophages from the EB myeloid differentiation To assess whether the transient immune activation affected phase (‘Macrophages’ n = 6142 vs ‘Macrophages_Day31plus7’ chromatin structure, we analysed the scATACseq macrophages n = 10,801) despite also having a similar number of expressed genes data. We found that 113 TF motifs were significantly enriched on (‘Macrophages’ n = 2312 vs ‘Macrophages_Day31plus7’ n = 2219). day31 + 7 ATAC peaks, compared to day31, and we obtained TF Finally, NMPs had a very low number of accessible peaks (n = 4071), transcriptional activity scores for 55 of these (Fig. 3e and in line with the low number of genes expressed in these cells Supplementary Data 6). Importantly, the top 11 enriched motifs (n = 1782) (Fig. 2h and Supplementary Fig. 6C). (enrichment score >5) at day31 + 7 corresponded to TFs that had been transcriptionally activated at day31 + 1 but were no longer activated at day31 + 7 (Fig. 3e). Indeed, the global TF motif Modifications in the last stage produce diverse macrophages. enrichment profile at day31 + 7 correlated with TF activities at Non-adherent cells at day31 were collected and plated in a fresh day31 + 1 (Pearson correlation: r = 0.45, p < 0.0003) but not at medium with cytokines for 7 days as step 3 of differentiation day31 + 7 (Pearson correlation: r = 0.11, p = 0.38) (Fig. 3e). In (Fig. 3a). We analysed the changes over time (time points: day31, short, this analysis suggested that activation on day31 + 1 was day31 + 1, day31 + 4, day31 + 7), compared the effect of multiple transient transcriptionally but seemed to induce chromatin cytokines (cytokines: M-CSF, GM-CSF, GM-CSF + IL-34) and changes still present on day31 + 7. Therefore, earlier macro- tested the effect of using foetal bovine serum (FBS) vs. defined phages may represent a more naive cell state, amenable to further medium (media: RPMI + FBS, StemPro34). We annotated the reprogramming in response to polarisation cues, which could cell types using LR from the foetal liver data set22 and they have distinct applications than later macrophages. clustered together across experiments (Fig. 3a, b). Next, we per- Activated macrophages can be classified as M1 or M2 formed a neighbourhood entropy analysis38 and observed that depending on whether they kickstart inflammation or resolve it, time points and media composition induced the highest diversity and the cytokines M-CSF and GM-CSF have been classically used within macrophages (homogeneity ROGUE scores: time points to induce these phenotypes, respectively44. We found that specific 0.53, media composition 0.45 and cytokines 0.7, Fig. 3d–f and TFs, such as MAF, ERG and LYL1, had reduced activity scores in Supplementary Fig. 7A). GM-CSF vs. M-CSF macrophages (Fig. 3f and Supplementary For the time point experiment, we analysed samples at day31, Data 6). In contrast, RFX5 showed increased activation in GM- day31 + 1, day31 + 4 and day31 + 7 using M-CSF standard CSF macrophages33. Despite GM-CSF being largely known to stimulation (Fig. 3a, b). The non-adherent cells collected at day31 promote PPARG activation36, this was not present in GM- from the EB myeloid differentiation phase were already mostly CSF + IL-34 culture conditions (Fig. 3f). Of note, IL-34 is macrophages, alongside the main myeloid populations and a essential for the development of microglia from embryonic small subset of fibroblasts (Fig. 3c). By day31 + 7 there was an myeloid precursors45, and GM-CSF + IL34 induces the microglial enrichment in macrophages, representing a total of 94.1% of cells phenotype on monocytes in vitro46. Both the knockdown and (Fig. 3c), which was consistent with the CD14 + /CD64 + cells pharmacological antagonism of PPARG promotes lipopolysac- analysed by FACS (Supplementary Fig. 7B). We observed broad charide (LPS)‐stimulated transition from the M1 to the M2 transcriptional differences between day31 and day31 + 1 macro- phenotype in primary microglia, with the concomitant upregula- phages, as well as between day31 + 1 and day31 + 4/7, while tion of markers such as CD206, TGFb and IL-447. macrophages from day31 + 4 and day31 + 7 overlapped (Fig. 3b). Furthermore, we stimulated the 4 macrophage subsets with We also analysed the changes in TF activity over time30,39 LPS for 6 h to study their TLR4 stimulation response (Fig. 4a). (Supplementary Data 6). From day31 to day31 + 1, the increased These included the macrophages obtained at the EB myeloid activity of JUN, FOS and NFKB1 confirmed the transient but differentiation phase, as well as those obtained after the 7-day robust immune activation identified earlier (Supplementary differentiation with distinct cytokines (M-CSF, GM-CSF and Fig. 5F). This transcriptional activation decreased during GM-CSF + IL-34). As a readout, we evaluated their single-cell day31 + 4 and day31 + 7 and some TF activities returned to transcriptomics profile and cytokine protein levels in the basal levels (e.g. STAT1, NFKB1 and RELA, Fig. 3d). Despite supernatant. Both analyses showed a clear distinction between being globally similar transcriptomically, we observed few but control and LPS-stimulated conditions (Fig. 4b). Known LPS- relevant differences between day31 and day31 + 7, including induced genes (e.g. IL6 or TNF) were overexpressed in all increased MITF and decreased SRF activities, which regulate macrophage populations (Fig. 4c and Supplementary Data 7). phagocytosis40,41, in addition to decreased SREBF1 and Similarly, TNFa and IL-6 protein levels were also upregulated SREBF2 activity, involved in lipid metabolism and macrophage (Fig. 4d). Regarding LPS specific effects, macrophages from the polarisation42,43. EB myeloid differentiation phase (equivalent to day31) 6 NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 ARTICLE a iPSCs > macrophages in vitro differentiation Computational analysis b LR dataset Cell Type Non-adherent FLi Fib Fibroblasts cells only +GMCSF+IL34 Mega Megakaryocytes FLi +GMCSF FLi Mast MastCells +MCSF In vivo publicly Mac FLi NMP NeutrophilMyeloidProgenitors available dataset FLi MDP MonocyteDCprecursors Popescu et al. Mo FLi Mo Monocytes Fetal Liver FLi Mac Macrophages Mac Differentiation Logistic regression All Time points day 31 31+1 31+4 31+7 Non-adh * * Fib Day31 Macrophage MDP Day31+1 Single cell RNA sequencing Differentiation Day31+4 Phase dataset NMP Day31+7 UMAP_2 31 Mega day Non-adh 31+1 31+7 Mac Differentiation Mast UMAP_2 n=52k cells, 3 donors, 4 timepoints +MCSF in StemPro34 +MCSF in RPMI+FBS UMAP_1 Non-adherent UMAP_1 cells only in red additional experimental conditions Time points Experiment Cytokines Experiment Media Experiment * MCSF-only time point c Time points Experiment 100 Percentage of cells Macrophages 75 UMAP_2 Monocytes UMAP_2 UMAP_2 MDP 50 NMP 3 donors, 3 donors, 1 donors, Mast cells 4 timepoints 2 timepoints 25 3 timepoints Megakaryocytes UMAP_1 UMAP_1 UMAP_1 Fibroblasts Day31 Day31 Day31 0 Day31+1 M-CSF 0 1 4 7 Day31+7 M-CSF Day31+1 RPMI+FBS Day31+4 M-CSF Day31+7 GM-CSF Day31+1 StemPro Days in Macrophage differentiation phase Day31+7 M-CSF Day31+7 GM-CSF+IL-34 Day31+7 RPMI+FBS other other Day31+7 StemPro other Macrophages TF activity * * * CEBPB * * * REL d Macrophages - Time points * * * IRF1 1 * STAT1 Day31 vs +1 vs +4 vs +7 * * * ATF1 0.5 e Correlation of TF motif enrichment (ATAC) with TF activation (RNA) * * SP1 * * * SPI1 0 Macrophages Macrophages * * HIF1A * * * ATF4 JUND JUND ATAC p eaks TF mot if enr ichment * * * TBP −0.5 FOSL2 FOS FOSL2 * * * FOXL2 BACH2 BACH2 FOS 20 20 FOSL1 * * * IRF9 FOSL1 * * * STAT2 JUN BATF JUNNFE2 Day31+7 vs Day31 * * * FOSL1 NFE2 * * * NFE2L2 BATF JUNB * * * FOS 15 JUNB 15 * * * JUN Enrich. * * NFKB1 * * RELA BACH1 >5 * * * ERG BACH1 10 * * * SMAD1 r=0.45 <5 10 r=0.11 * * * TCF7L2 * * * POU2F2 p<0.0003 NFE2L1 NFE2L1 p=0.38 * * * SP3 * * * SMAD3 5 5 * * * FOSL2 UMAP_2 * * * RELB * * * SREBF1 * * * SREBF2 0 0 * * * ELK1 ROGUE score = 0.53 * * * SRF −0.5 0.0 0.5 1.0 −0.50 −0.25 0.00 0.25 * * * MITF UMAP_1 * * * RUNX3 TF transcriptional activity TF transcriptional activity Day31 * * * POU5F1 Day31+1 Day31+7 Day31+1 M-CSF * * * SIX2 Day31+4 M-CSF * * GATA6 * * * ASCL1 Day31+7 M-CSF * * * TWIST1 other Macrophages TF activity f Macrophages - Cytokines * * * CREM 1 * * * FOXK2 M-CSF vs GM-CSF vs GM-CSF+IL-34 * * * KLF6 0.5 * * * ELF1 * * * MAX 0 * * * LYL1 * * * ATF7 −0.5 * * * NR2F1 * * * SMAD4 * * * IRF3 * * * E2F1 * * * NRF1 * * E2F4 * * * NR2C2 * * ATF6 * * * MAF * * * SREBF1 * * * SREBF2 * * * PPARG * * * PPARA * * * KLF4 * * * FOXO1 * * * TWIST1 UMAP_2 * * * ELK1 * * * RARA * * * ERG * * * USF2 ROGUE score = 0.70 * * HNF4G * * * MITF UMAP_1 * * * FOXL2 * * * SRF Day31 * * KLF13 * * * MXI1 Day31+7 M-CSF * * RFX5 Day31+7 GM-CSF * * * POU2F1 Day31+7 GM-CSF+IL-34 * * * POU2F2 other * * * SP3 * * * TFAP2C overexpressed metallothionein and matrix metalloproteinase complex member NLRP3, which is involved in the cleavage and genes after LPS stimulation (e.g MT1E and MMP10), along with activation of IL-1B. Accordingly, IL1B was upregulated both at increased IL-10 protein levels (Fig. 4c, d). Macrophages after the RNA and protein levels in this macrophage subset upon LPS 7 days in GM-CSF showed transcriptomic upregulation of CCL1 stimulation (Fig. 4c, d). and CXCL9 cytokines with LPS (Fig. 4c). Macrophages after Finally, we tested whether using a defined medium (StemPro34 7 days in GM-CSF + IL-34 upregulated the inflammasome serum-free media, SP-SFM) at this step would prevent the NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 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. transient macrophage activation. TF activity analysis showed SP- thymus21 CCR7 DC LR model (Supplementary Fig. 8) and the SFM induced similar activation signals (e.g JUN, FOS and expected marker expression (e.g., CCR7, LAMP3 and lack of NFKB1). However, most TF activities were only partially CD14; Fig. 5b, c and Supplementary Fig. 9D). This result showed recovered by day+7, including PPARG activity which did not the ability of iPS-derived DCs to adopt known activated DC significantly decrease from day+1 level. One additional difference profiles described in vivo52. Finally, cells obtained with the DC was the maintained activity of the SREBF1 and SREBF2 TFs in protocol showed dendrite-like structures in contrast to iPSC- SP-SFM, which link lipid metabolism to macrophages' inflam- derived macrophages (Supplementary Fig. 9E). matory response (Supplementary Fig. 7A and Supplementary As in the presence of M-CSF + IL-3 (macrophage protocol) Data 6). These findings showed how media affects macrophage (Fig. 2f, g), myeloid cells appearing in the presence of GM- metabolism and function, which can expand the use of this CSF + FLT3L (DC protocol) activated similar TFs when com- model42,43. pared to their YS and foetal liver in vivo counterparts (Fig. 5d and Supplementary Data 9). Moreover, in vitro cDC2s activated PU.1 GM-CSF + FLT3L iPSC differentiation produces dendritic (SPI1 gene)53 and KLF454 TF networks relevant for in vivo cDC2 cells. Conventional dendritic cells (cDC) present antigens to T cells identity (Fig. 5d). We also observed increased RFX5 activity, and act as messengers between innate and adaptive immunity48. which regulates MHC II gene expression33. A recent study Protocols to induce DC differentiation in vitro are based on sup- postulated a role for CEBPB in the control of DC maturation and plementing factors, including GM-CSF and FLT3L, that act coop- later stages of DC commitment55. Our results showed reduced eratively on cell precursors to drive cDC generation49,50. To produce CEBPB activity in cDC2 cells compared to monocytes (in vivo iPSC-derived DCs we used GM-CSF + FLT3L in the EB myeloid and in vitro, Fig. 5d), which indicates that the in vitro phenotype differentiation phase (instead of M-CSF + IL3) and GM-CSF + IL4 shares features with a functionally mature DC subset charac- in the last phase of differentiation on non-adherent cells from day31 terised by upregulation of costimulatory and MHC class II (instead of M-CSF) (Fig. 5a). We then performed single-cell tran- molecules. scriptomic analysis and annotated cells using LR classifiers that were We analysed step 3 of the DC protocol in detail, from the non- trained on gastrulation23, YS6 as well as foetal liver and thymus21,22 adherent cells produced during EB myeloid differentiation until data sets (Fig. 5a, Supplementary Fig. 8 and Supplementary Data 8). the end of DC differentiation phase. A mean of 47.5% (standard GM-CSF + FLT3L followed by GM-CSF + IL4 stimulation deviation = 3.67) of the cells produced in the last three time produced DCs with a transcriptomic profile similar to cDC2 points were cDC2 cells (Fig. 5f). Thus, this differentiation was less and a small DC population resembling CCR7 + DCs. These efficient than the macrophage protocol, where macrophages populations were exclusive to this protocol and were not found in represented 94.1% of cells by day31 + 7 (Fig. 3c). In contrast to the Discovery data set (M-CSF + IL-3 followed by M-CSF, what was observed with macrophages, the proportion of cDC2 Fig. 1c). All other cell types in the myeloid and stromal cell remained stable (Figs. 3c and 5f). Regarding the TF activity compartments, and their temporal dynamics (Supplementary analysis, the in vitro activation of cDC2 and macrophages Fig. 9A), remained the same as the macrophage protocol (Figs. 5b induced shared regulatory programmes, including activation of and 1c). Marker gene expression analysis supported the cell type NFKB1, RELA or JUN on day+1 of this phase (Figs. 3c and 5e annotations, and cDC2 expressed bonafide DC markers (e.g., and Supplementary Data 9). Interestingly, the particular profile of HLA-DR, CD1C, CLEC10A)48 (Fig. 5c). The identity of cDC2 was TF networks induced in cDC2 on that day (including JUN, REL, further explored using LR models built on adult tonsil cDC2s and SP1 and HIF1A) did not clearly return to basal levels by two subsets of monocyte-derived DCs (moDCs) from adult day31 + 7 (Fig. 5f and Supplementary Data 9). This was in ascites51. To match the cell types present in the adult data set we contrast to what was observed in the macrophage differentiation created a mixed manifold including the samples from day31 to phase (Fig. 3d and Supplementary Data 6). day31 + 7 from both macrophages and DC protocols. The tonsil To confirm the DC phenotype of the iPSC-derived cells in cDC2s LR model showed the highest prediction probability on our new protocol, we analysed canonical DC markers by FACS the cells identified as cDC2s on the mixed manifold (tonsil cDC2 and functionally interrogated the cells using an antigen model AUC 0.837, Supplementary Fig. 9B, C), thus supporting processing (DQ-OVA) and T cell activation assays. We our annotation. One of the moDCs subsets from ascites, which observed positive cells for DC markers (CD1C, CD209, showed blood cDC2 features, overlapped with our in vitro cDC2 CD11c, HLA-DR, CD86) and low levels of the macrophage cells with a lower AUC (blood cDC2-moDC AUC 0.756 marker CD14, thus validating the DC phenotype56 (Fig. 5g Supplementary Fig. 9B, C). This result demonstrated our in vitro and Supplementary Fig. 9F). Functionally, DQ-OVA antigen cells have cDC2 features, although we cannot fully discard a processing is higher in adult moDCs at earlier time points (15 moDC identity. The small subpopulation annotated as CCR7 + to 30 min), whereas cDC processing capacity peaks at DCs showed a strong and specific prediction using a foetal 60 min57. Accordingly, iPSC-derived cDC2 showed DQ-OVA 8 NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 ARTICLE 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. NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 processing at longer time points only (45 and 60 min) CD4 + T cells purified from peripheral blood mononuclear (Fig. 5h). While this time distinction is not commonly used cells. iPSC-derived DCs showed a stronger ability to induce T to distinguish moDCs from cDC2s, it shows that the antigen cell proliferation (Fig. 5i). These results demonstrated the DCs processing behaviour of our iPSC-derived DCs resembles that generated in our newly developed protocol were transcrip- of cDC2s. Finally, the cells produced at the end of the DC and tionally and functionally similar to in vivo cDC2, and could be macrophage differentiation were co-cultured with human used as a model to study DC function in vitro. 10 NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 ARTICLE 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 set6 and foetal liver, skin and kidney data set22. 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. GWAS immune phenotype genes shape iPSC-derived myeloid suppressor cells have decreased levels of PRKCB, which dampens cells. Dysfunctional myeloid differentiation and signalling can DC differentiation and function in vivo60. lead to immune-related disorders58. To dissect potential myeloid Macrophages generated from iPSCs deficient in PRKCB, LSP1 contributions involved in these pathologies, we selected four or ICAM1 exhibited a mixed anti-inflammatory and anti-fibrotic genes linked to immune-related GWAS hits (ICAM1, LSP1, phenotype. KO macrophages upregulated the suppressors of the PRKCB and ZEB2) based on the existing literature. The ICAM1 NFkB-dependent inflammatory pathway KLF371 and ATF372 and LSP1 loci contain SNPs linked to autoimmune inflammatory (Fig. 6g). They also decreased the expression of M2 profibrotic diseases by GWAS (https://www.ebi.ac.uk/gwas/home) and phenotype genes (e.g. FN1, GRN and SPP1, Supplementary interact with each other59. PRKCB is a protein kinase associated Data 11), as well as decreased activity of M2-promoting with inflammatory diseases and blood cell counts in GWAS and is transcription factors (e.g. MAF34 and PPARG73, Fig. 6g). We involved in myeloid DC differentiation60. Finally, ZEB2 is found observed a connection between PRKCB, LSP1 and ICAM1, as in GWAS for blood phenotypes, and regulates hematopoiesis in ICAM1 was downregulated in PRKCB and LSP1 KOs (Supple- mice61 and DC cell fate decisions62. To study their involvement mentary Data 11). This suggested that the 3 genes are part of a in myeloid differentiation and identity, we generated knock-out regulatory network that fine-tunes macrophage phenotype and (KO) iPSC lines using CRISPR/Cas9 in one of our donor cell lines represses tissue healing both by promoting REL/p65-mediated (kolf_2) (Supplementary Fig. 10). inflammation and controlling the expression of profibrotic M2 KO iPSC lines were differentiated into macrophages and DC genes. Notably, these KO lines generated a macrophage alongside wild-type (WT) isogenic lines. Single-cell transcrip- population with regenerative potential due to the concomitant tomic analyses were performed at day0 (iPSC stage) and day31 suppression of the tissue remodelling programmes (Fig. 6g) and (EB myeloid differentiation phase) (Fig. 6a, b). Though we the NFkB pathway. Our data illustrate how such broad defects in observed all cell populations in all conditions (Fig. 6c, d), some of macrophage polarisation could impair the proper resolution of the KOs affected the cell type proportions (Supplementary inflammation, leading to in vivo autoimmune inflammatory Fig. 11A). As expected, knocking out ZEB2 reduced myeloid disorders74. cells to 5.5% versus 65.6% in WT lines in the macrophage Altogether, we found deletion of genes associated with protocol (12-fold decrease), and 2.1% versus 20.4% in WT lines in autoimmunity mediated early hematopoiesis (e.g. ZEB2) and the DC protocol (10-fold decrease) (Supplementary Data 10). The influenced the inflammatory potential (e.g. ICAM1, LSP1 and overall reduction in cell numbers in the ZEB2 KO vs the other PRKCB) of iPSC-derived myeloid cells, coinciding with the lines suggested that ZEB2 absence either blocked myeloid expected phenotype. Therefore, iPSC-differentiation protocols differentiation or induced apoptosis as described for myeloid can be powerful tools to functionally interrogate GWAS hits and leukaemic cells63. On the contrary, PRKCB KO increased the other genes of interest in vitro. proportion of myeloid cells in the DC differentiation protocol (20.4% in WT versus 78% in PRKCB KO, fourfold increase) Discussion (Supplementary Data 10). The other KOs did not markedly The full characterisation and assessment of the robustness, influence myeloid cell proportions. accuracy and efficiency of in vitro protocols are essential to Some myeloid cell types showed transcriptomic differences utilising them as models for disease as well as leveraging them between KO lines. LSP1 and ICAM1 KOs monocytes from the DC in the search for novel therapeutic targets. Myeloid cells are protocol had a transcriptomic profile with intermediate mono- central to immunity and participate in major inflammatory cytes features. They were characterised by the upregulation of the and autoimmune disorders. Experimental protocols to generate HLA genes and downregulation of the S100 gene family64 macrophages that are easy to replicate and amenable to scaling (Fig. 6e). Interestingly, KO monocytes from the macrophage up are paramount to studying human macrophage ontogeny, protocol did not show this profile (Supplementary Fig. 11B). An genetics and function in health and disease. Here, we profiled increased population of intermediate monocytes has been more than 470k single cells across a commonly used, described in many autoimmune diseases such as active Crohn’s straightforward differentiation process from human iPSC to disease65 and rheumatoid arthritis66. Monocytes from the PRKCB terminally differentiated macrophages. We reconstructed, KO cell line after the DC differentiation protocol had a myeloid- using cell trajectories, the in vitro sequence of events, which derived suppressor profile, including low expression of HLA-DR parallel foetal hematopoiesis prior to the establishment of and CD74 with high levels of CCL2 and MMP967–70 (Fig. 6f). This HSC. Moreover, we show this protocol is a valuable resource in was consistent with the observation that myeloid-derived studying different macrophage cell states, multiple myeloid NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 populations, including erythrocytes, megakaryocytes and mast To quantitatively assess the accuracy of our in vitro models, we cells that spontaneously arise, as well as DC, induced by used machine learning tools. We built LR models trained on adjusting the media composition. Finally, we used this model scRNAseq data from developmental atlases mapping the formation to interrogate the functional effect of genes associated with of the immune system and projected them onto the in vitro data inflammatory and autoimmune disorders and interpreted the sets. The computational framework we have established could be results in relation to their in vivo counterparts. adapted to annotate multiple iPSC differentiation protocols. 12 NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 ARTICLE 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. Following this strategy, we found that the initial phases of (KC)81,84 the LR classifier from foetal liver KC did not capture iPSC–macrophage differentiation faithfully recapitulate YS foetal any in vitro subpopulation, indicating that the strong tissue- hematopoiesis and generate foetal-like FOLR2+ macrophages. The resident signature of KC is not recapitulated using this protocol. lack of an HSC cluster in our data, the activation of master regulator In particular, we did not observe activation of KC-determining RUNX1 in the endothelial-to-hematopoietic transition75, and the TF LXRa, RBPJ or SMAD4, probably due to a lack of liver- lack of expression of HOXA genes in the myeloid progenitors25 all derived signals such as Notch ligand DLL4 essential for their suggest the iPSC-derived macrophages protocol closely recapitulate induction85. All considered, the ability to derive certain macro- YS differentiation prior to the establishment of definitive hemato- phage subtypes in vitro, combined with high-resolution single- poiesis. Thus, we propose this model as a unique system for inter- cell analysis provides the unprecedented possibility to directly rogating the early stages of hematopoietic differentiation in humans, manipulate these cells and interrogate the extent of intrinsic which are largely unexplored. macrophage plasticity, which remains a matter of debate86. Our study shows that iPSC–macrophage differentiation gen- Additionally, a marked enrichment of the macrophage population erates a wide range of myeloid cells and presents a detailed list of was observed during the macrophage differentiation phase, TFs that mediate the generation of distinct myeloid progenitors pointing to the need for single-cell approaches at distinct time in vitro. We also observed that EB stimulation with GM-CSF and points. Overall, the exhaustive analysis of this in vitro protocol FLT3L produces cell populations resembling cDC248. However, will guide the applications of this system to model macrophages. these in vitro DCs may also be related to moDCs, due to the Finally, we leveraged this model to experimentally evaluate presence of monocytes in the culture and the resemblance with an genes linked to immune-related phenotypes by GWAS. Interest- in vivo moDC subset, which shows features of blood cDC2s51. ingly, CRISPR/Cas9-mediated KO of GWAS hits (PRKCB, LSP1 Additionally, a small subpopulation of DCs in the culture reca- and ICAM1) in iPSC-derived macrophages and DCs highlighted pitulates the CCR7+ in vivo activation state52, highlighting the their potential role in physiological and pathological cell states of functional potential of iPSC-derived DCs. Nevertheless, the pro- distinct cell types. Specifically, macrophage differentiation in KO tocol for DC generation had a lower efficiency than the macro- lines altered inflammatory and extracellular matrix genes. phages protocol, which is consistent with the known bias of Fibrosis constitutes a pathological feature of most chronic YSMP towards macrophages6,20–23. inflammatory diseases including the ones featured in our Interestingly, we found erythrocytes in our data following a study87,88, and our results open an avenue for therapeutic inter- decline in RUNX1 and SPI1 TF activity76 and a rise in GATA1 vention in these disorders. In line with this, we show that both the activation77. The in vitro iPSC differentiation of erythroid-lineage foetal-like FOLR2+ and the SPP1+ TAM states observed in liver cells has relevant biomedical implications in the study and cancer27 are recapitulated in this system. This suggests these cells treatment of anaemias78. Monocytes are also generated during could also be a faithful model to unravel the role of macrophage differentiation but, despite the complexity and limitations of subtypes in the tumour microenvironment. trajectory analyses79, our analyses consistently indicate that In conclusion, we have defined a comprehensive map of cells produced macrophages bypass the monocyte stage. It is tempting and molecular programmes that underlie iPSC–macrophage dif- to speculate that monocytes from this protocol can also be dif- ferentiation in a dish. Macrophages play an important role for ferentiated into macrophages and polarised to specific immunity in health and disease, and represent key cellular targets functions80,81. However, the lack of surface proteins uniquely for immunotherapy. Our study shows the potential of deeply expressed in monocytes (vs macrophages) makes this study more characterising differentiation protocols at the single-cell level, and challenging. Future work should evaluate if the origin of mac- demonstrates that these iPSC myeloid differentiation protocols rophages imprints on their function. are valuable models for interrogating the very early stages of Macrophages are a heterogeneous cell type and epigenetic hematopoietic formation that have been largely unexplored so far. states are instrumental in the generation of functional and phe- notypic diversity81–83. Nevertheless, proper models to study this Methods heterogeneity are missing. Here, we demonstrated that exposure Reagent/material details can be found in the Supplementary Information. to M-CSF generates macrophages with a transcriptomic profile similar to their unstimulated counterparts, but with a distinct chromatin accessibility landscape. Our combined transcriptome/ Human induced pluripotent stem cell lines. All iPSC lines used in the study were generated by the HIPSCI project. Details on their generation are available at http:// ATAC analysis shows the potential of macrophages, among non- www.hipsci.org. Briefly, we used kolf_2, yemz_1 and vass_1 in the Discovery and adherent cells in the EB myeloid differentiation phase, to be truly DC data sets, and we added ceik_1, eesb_1 and wegi_1 for the Validation data sets. naïve cells early on, whereas 24 h into the macrophage differ- All cells in the KO data set are derived from kolf_2 as a parental line, which was entiation phase, macrophages are reversibly activated. We also subcloned prior to editing (kolf_2_C1). All HIPSCI samples were collected from consenting research volunteers recruited from the NIHR Cambridge BioResource identified multiple features resembling YS myelopoiesis in our (http://www.cambridgebioresource.org.uk), initially under existing ethics rules for iPSC system. Nonetheless, despite the YS origin of Kupffer cells iPSC derivation (Regional Ethics Committee (REC) reference 09/H0304/77, v.2, 4 NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 January 2013), with later samples collected under a revised consent (REC reference Genomics), aiming for 10000 nuclei per sample. Library preparation was carried 09/H0304/77, v.3, 15 March 2013). out according to the manufacturer’s protocol and sequenced on Illumina NovaSeq 6000, aiming for 20000 fragments per cell using the sequencing formats; read 1: 50 cycles; i7 index: 8 cycles, i5 index:16 cycles; read 2: 50 cycles. In vitro differentiation to macrophages and dendritic cells. We used an adaptation of the van Wilgenbrug et al. protocol15. Feeder-free human iPSC lines were cultured in E8 (StemCell Technologies) on vitronectin-coated plates (Life Single-cell RNA seq computational analysis. Cell Ranger (v3.1.0), mapping to Technologies). For the embryoid body (EB) formation, step1, a single-cell sus- GRCh38 (v3.0.0), was used to filter out empty droplets using default values. Cells pension of hiPSC was plated in 100 µl of EB medium – E8 + SCF (20 ng/ml, were further filtered out for the number of genes (<200) and percentage of mito- Peprotech) + VEGF (50 ng/ml, Peprotech) + BMP-4 (50 ng/ml, Peprotech) + chondrial RNA (>8.5%) using Seurat (https://satijalab.org/seurat/ v3.2.2), R v4.0. ROCK inhibitor (10 µM, Sigma) – at a density of 10,000 cells per well in round All cells identified as doublets using SoupOrCell89, DoubletDetection (https://doi. bottom low-attachment 96 well plates (Corning). After 2 days, we changed half the org/10.5281/zenodo.2678041) and Scrublet90 were discarded. Cell genotype calling media (50 µl) and replaced it with fresh EB media. At day4, EB myeloid differ- was performed using SoupOrCell89. All data sets were normalised using sctrans- entiation started, step 2, when EBs were plated in gelatin-coated 6-well plates form in Seurat91, and UMI counts, mitochondrial RNA and cell cycle variables (Sigma Aldrich) at a density of 8–10 EBs per well in EB-Mac medium – StemPro- were regressed out by cell line. Multiple hiPSC lines were integrated using Seurat’s 34 (Life Technologies) + M-CSF (100 ng/ml, Peprotech) and + IL-3 (25 ng/ml, anchor-based method24. After PCA dimensionality reduction and louvain Peprotech). The EB-Mac medium was changed every 4 to 5 days. At day31, step 3, clustering92, data sets were further analysed as described below. non-adherent cells were collected by centrifugation from the culture media and 1×106 cells were cultured in 10 cm tissue culture plates for 7 days in macrophage Single-cell ATAC seq computational analysis. Cell Ranger ATAC pipeline differentiation medium – RPMI (ThermoFisher) + 10% heat-inactivated FBS + M- (v1.2.0), mapping to GRCh38 (v3.0.0), was used for read filtering and barcode cell CSF (100 ng/ml, ThermoFisher). calling. Peaks were re-called using cellatac, an in-house implementation of Cusano- Alternative macrophage differentiation media were used in the macrophage vich’s approach93 (https://github.com/cellgeni/cellatac)94. Peak and cell filtering were differentiation phase, step 3. For the cytokines experiment (Fig. 3b–d), we used performed using cellatac and Signac (https://satijalab.org/signac/ v1.1.1), as described RPMI + 10% heat-inactivated FBS + GM-CSF (50 ng/ml, Peprotech) and in Supplementary Fig. 1B. Normalisation and dimensionality reduction were per- RPMI + 10% heat-inactivated FBS + GM-CSF (10 ng/ml) + IL-34 (100 ng/ml, formed using term frequency-inverse document frequency (TF-IDF) and Singular Peprotech), for the latter cells were plated at 6×105 cells per well of a 6well plate. Value Decomposition (SVD), respectively. SLM from Seurat was used for clustering. For the media experiment (Fig. 3b and Supplementary Fig. 7A), fully defined TF motif analysis was performed using Signac and JASPAR 202095 motifs database. medium – StemPro-34 + M-CSF (100 ng/ml) - was used. Step 1 is shared between macrophages and DCs, while different cytokines are used in steps 2 and 3. For DC differentiation EBs at day4 are plated with EB-DC Cell-type annotation of RNA and ATAC data sets. Both the Discovery and DC media – StemPro-34 + GM-CSF (50 ng/ml) + FLT3L (100 ng/ml, Peprotech) in data sets were annotated using logistic regression (LR) models built on publicly the same types of plates and density as the macrophage protocol. At day31, step 3, available single-cell transcriptomic data sets. The LR prediction models used at non-adherent cells were collected and plated in 10 cm tissue culture plates in DC each step were built based on a general linear model function and a tenfold cross- differentiation medium – RPMI + 10% heat-inactivated FBS + GM-CSF (50 ng/ validation. Briefly, public raw data (Cell Ranger output when available or processed ml) + IL-4 (100 ng/ml, Peprotech). matrices otherwise) was downloaded and re-processed as described for the Dis- covery and DC data sets. Then, public data sets were split into training (70%) and test (30%) sets, ensuring these proportions were accounted for each cell type. Then LPS stimulation of iPSC-derived cells. Four populations of cells derived from the we generated LR models to classify each cell type for each gene in the training macrophage protocol described above were stimulated with LPS (i.e. non-adherent partition of the in vivo data set using normalised data. A ranked gene list based on cells from the EB myeloid differentiation phase -Mac Day0-, cells at the end of the the area under the curve (AUC) of each gene was produced for each cell type. The Macrophage differentiation phase incubated for 7 days with M-CSF (100 ng/ml) optimal number of genes to build the final LR classifier was chosen by building -Mac Day7 M-CSF-, with GM-CSF (50 ng/ml) -Mac Day7 GM-CSF- and with models on the training set and calculating the AUC of the prediction on the test set. GM-CSF (10 ng/ml) + IL-34 (100 ng/ml) -Mac Day7 GM-CSF IL-34-). For each This was repeated with an increasing number of genes down the ranked list sample, LPS (Sigma Aldrich) was added to the media in two wells for a final described above. The number of genes that produced a model on the training set concentration of 2.5 ng/ml while a third well was kept as control. After 2 h, Bre- with the highest AUC when applied to the test set was then used to build the final feldin A (Sigma Aldrich) was added to reach 5ug/ml to all wells and cells were model on the full in vivo data set. This LR prediction model was then used to incubated for an additional 4 h (total 6 h of stimulation with LPS). Media was classify the cells in the in vitro data set. Finally, the mean prediction probability per collected for control and LPS wells for cytokine analysis using V-PLEX Proin- louvain cell cluster was calculated for all the LR models built, and each cluster was flammatory Panel 1 Human Kit (Mesoscale) according to manufacturing instruc- labelled based on the LR model with the highest mean prediction. As an estimate of tions. In parallel, cells were collected using 10 mg/ml Lidocaine (Sigma Aldrich) the strength of the association between the annotation and the labelled clusters, we 2 mM EDTA (ThermoFisher) solution for 5 min at 37 °C, the two wells of LPS calculated the AUC of each annotated cell type based on the LR probabilities. stimulation per population were combined. Single-cell RNAseq analysis was per- For the validation data set, cell type annotations from the Discovery data set formed aiming at 3000 cells per condition as described in ‘10× Genomics Chro- were projected on the transcriptomic and ATAC validation data sets using Seurat’s mium GEMs sample preparation and sequencing’. anchor-based label transfer approach24. Trajectories analysis. Spliced/unspliced RNA expression matrices were generated 10× Genomics Chromium GEM sample preparation and sequencing. Single-cell using the command line tool from velocyto (http://velocyto.org/velocyto.py/tutorial) on transcriptomic analysis on iPSC-to-macrophage differentiation was performed in 3 Python v3.8. scVelo was used for trajectory analysis based on RNA velocity and PAGA iPSC lines for the Discovery data set and 6 hiPSC lines for the Validation data set. graph abstraction as described (https://scvelo.readthedocs.io/DynamicalModeling/ One 6-well well per line was collected using TrypLE (Life Technologies) at 20 time v0.2.2) on Python v3.8. All analyses were performed on a per sample basis. Progenitor points in the Discovery data set, and 2 6-well wells per line at 7 time points in the and derivative cell type associations, shown in Fig. 2 a, were the most prevalent con- Validation data set, between day0 and day38 (day31 EBs +7 days of the macro- nections observed in the scVelo analysis across all samples analysed. Additionally, the phage differentiation phase). At every collection day, all cells from the wells or proposed progenitors shown were observed earlier or simultaneously to the derivative plates were detached, merged, counted, passed through 40 µM filters and resus- cell type along the time series. If the direction inferred with the scVelo analysis con- pended in DPBS (ThermoFisher) + 0.4% BSA (GE Healthcare). At the day31 tradicted the order of appearance of the potential progenitor and derivative cell types, collection day, we additionally collected sample ‘Day31 non-adh’ which corre- the association was not considered. sponds to the non-adherent fraction of the culture, it was processed as the rest but no detachment step was necessary. Cell suspensions were processed using the Chromium Single Cell 3’ kit (v2 for Discovery, v3 for Validation, 10X Genomics), Transcription factor activity analysis. Transcriptomic changes across trajectories aiming at recovering from 3000 to 10000 cells. Library preparation was carried out and time points were studied based on transcription factor activities using DoR- according to the manufacturer’s instructions. Libraries were sequenced, aiming at a othEA and VIPER analysis30. DoRothEA v1.2.1 (https://saezlab.github.io/ minimum coverage of 50000 raw reads per cell, on the Illumina HiSeq 4000 dorothea) required Seurat v4.0.2 (https://satijalab.org/seurat/). Both in vitro and (Discovery) or Novaseq 6000 (Validation) using the sequencing formats; read 1: 26 in vivo data sets were subset based on connected cell types according to the cycles; i7 index: 8 cycles, i5 index: 0 cycles; read 2: 98 cycles (3’ kit v2) or read 1: 28 trajectory analysis. Normalised data was scaled within each subset, and TF activity cycles; i7 index: 8 cycles, i5 index: 0 cycles; read 2: 91 cycles (3’ kit v3). scores were computed for each cell for 271 TFs with high-confidence target-gene Sample preparation and sequencing for the DC data sets were performed as annotation (A, B and C confidence levels, https://saezlab.github.io/dorothea/). described for the Discovery data set. The Knockout data set samples were processed as Heatmaps for in vitro vs in vivo comparison were produced by selecting the top 50 described for the validation data set but only for 2 time points (i.e. day0 and day31). most variable TFs in each data set, and results were merged and plotted using Single-cell ATAC analysis was performed in a subset of the single-cell pheatmap (https://www.rdocumentation.org/packages/pheatmap/versions/1.0.12). suspensions for 6 of the time points of the Validation data set described above. Significance analysis was performed with FindMarkers function from Seurat Single-nuclei suspensions were obtained and processed according to the (https://satijalab.org/seurat/, v4.0.2) using the LR method, p values were adjusted manufacturer’s instructions using Chromium Single Cell ATAC v1.0 (10X using Bonferroni correction. 14 NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 ARTICLE Marker protein and antigen processing DQ-OVA assay analysis by FACS. Received: 7 July 2021; Accepted: 6 May 2022; Macrophages and dendritic cells were detached from 10 cm plates using 10 mg/ml Lidocaine 2 mM EDTA solution for 5 min at 37 °C, collected in DPBS and spun down at 300 × g for 3 min. Samples were then fixed with BD Cytofix/Cytoperm buffers (ThermoFisher) for 20 min at room temperature and washed with DPBS + 1% FBS. Staining with fluorescent-labelled primary antibodies (see Sup- plementary Information for antibody details) was performed in the dark at room temperature for 30 min. After 2 washes with DPBS + 1% FBS, cells were analysed References by FACS in a BD LSR Fortessa II and data was analysed using FlowJo v10. The 1. Okabe, Y. & Medzhitov, R. Tissue biology perspective on macrophages. Nat. gating strategy is provided in Supplementary Fig. 4B. Immunol. 17, 9–17 (2016). Dendritic cells were collected as described above and incubated with DQ-OVA 2. Hoeksema, M. A. & Glass, C. K. Nature and nurture of tissue-specific (ThermoFisher) in the dark at 4 °C or 37 °C for 15 min, 45 min and 60 min, as macrophage phenotypes. Atherosclerosis 281, 159–167 (2019). indicated. Cells were then washed with ice-cold DPBS + 1% FBS and analysed by 3. Lavin, Y., Mortha, A., Rahman, A. & Merad, M. Regulation of macrophage FACS as above. development and function in peripheral tissues. Nat. Rev. Immunol. 15, 731–744 (2015). 4. Gordon, S. & Plüddemann, A. Tissue macrophages: heterogeneity and T cell activation assay by CFSE staining and FACS. Twenty thousand iPSC- functions. BMC Biol. 15, 53 (2017). derived non-adherent cells from the macrophage or DC differentiation protocol 5. Hoeffel, G. & Ginhoux, F. Fetal monocytes and the origins of tissue-resident were differentiated to DC (in presence of GM-CSF and IL-4) or macrophages (in macrophages. Cell. Immunol. 330, 5–15 (2018). the presence of M-CSF) for 7 days. This was performed as described for the last 6. Bian, Z. et al. Deciphering human macrophage development at single-cell phase of differentiation in the protocols used throughout the study but in round- resolution. Nature 582, 571–576 (2020). bottom 96 well plates. In parallel, peripheral blood mononuclear cells (PBMCs) 7. Kaczanowska, S. et al. Genetically engineered myeloid cells rebalance the core were isolated using Ficoll-Paque PLUS (GE Healthcare) density gradient cen- immune suppression program in metastasis. Cell 184, 2033–2052.e21 (2021). trifugation. Naïve T cells were isolated from the PBMC fraction using EasySep® human naïve CD4 + T cell isolation kits II (StemCell Technologies) according to 8. Frame, J. M., McGrath, K. E. & Palis, J. Erythro-myeloid progenitors: the manufacturer’s instructions. After the 7 days of DC and macrophage differ- ‘definitive’ hematopoiesis in the conceptus prior to the emergence of entiation described above, 100.000 naïve T cells, stained with CFSE (Thermo- hematopoietic stem cells. Blood Cells Mol. Dis. 51, 220–225 (2013). Fisher), were stimulated with anti-CD3 (Tonbo Biosciences, 70-0037-U100, 100 ng/ 9. Ivanovs, A. et al. Highly potent human hematopoietic stem cells first emerge 1million cells) in the presence of iPSC-derived dendritic cells, iPSC-derived mac- in the intraembryonic aorta-gonad-mesonephros region. J. Exp. Med. 208, rophages or anti-CD28/CD3 (StemCell Technologies, 10971, 6ul/1million cells). 2417–2427 (2011). After 4 days of co-culture, T cell proliferation was measured by assessing CFSE 10. Holt, P. G. & Jones, C. A. The development of the immune system during dilution. Anti-CD3 alone was tested independently and had no effect on T cell pregnancy and early life. Allergy 55, 688–697 (2000). proliferation rate. The experiment was performed in two independent experiments 11. Freyer, L., Iturri, L., Biton, A. & Perdiguero, E. G. Overlapping definitive using two PBMC donors each time in one of the iPSC lines (kolf_2). The gating progenitor waves divide and conquer to build a layered hematopoietic system. strategy is provided in Supplementary Fig. 9G. Preprint at bioRxiv https://doi.org/10.1101/2020.12.24.424302 (2020). 12. Lee, C. Z. W., Kozaki, T. & Ginhoux, F. Studying tissue macrophages in vitro: are iPSC-derived cells the answer? Nat. Rev. Immunol. 18, 716–725 (2018). CRISPR-Cas9 KO of human induced pluripotent stem cell lines. KO iPSC lines 13. Alasoo, K. et al. Transcriptional profiling of macrophages derived from were generated by substituting an asymmetrical exon with a Puromycin cassette and monocytes and iPS cells identifies a conserved response to LPS and novel expanding those clones with a frame-shift indel in the remaining allele. A hSpCas9 and two small guide RNA expression vectors along with a homology-directed repair alternative transcription. Sci. Rep. 5, 12524 (2015). template vector were used. The template vector harboured an EF1a-Puromycin 14. Cao, X. et al. Differentiation and functional comparison of monocytes and cassette with two flanking 1.5 kb homology arms designed around the asymmetric macrophages from hiPSCs with peripheral blood derivatives. Stem Cell Rep. exon of interest. For each KO line, 2 × 106 iPSC single cells were transfected using 12, 1282–1297 (2019). the Human Stem Cell Nucleofector® Kit 2 (Lonza) with 4 μg, 3 μg and 2 μg of each 15. van Wilgenburg, B., Browne, C., Vowles, J. & Cowley, S. A. Efficient, long term plasmid, respectively, and plated in 10 cm plates. After 72 h, cells were selected in production of monocyte-derived macrophages from human pluripotent stem 3 μg/mL Puromycin (Invivogen) and colonies were expanded and genotyped. Lines cells under partly-defined and fully-defined conditions. PLoS ONE 8, e71098 confirmed to have the Puromycin cassette and the presence of a frame-shift indel by (2013). Sanger sequencing (Supplementary Information) were selected for the experiments. 16. Gutbier, S. et al. Large-scale production of human iPSC-derived macrophages for drug screening. Int. J. Mol. Sci. 21, 4808 (2020). 17. Camp, J. G., Wollny, D. & Treutlein, B. Single-cell genomics to guide human KO iPSC-derived cell types analysis. Transcriptomic alterations between cell stem cell and tissue engineering. Nat. Methods 15, 661–667 (2018). types arising in WT and KO lines were assessed using differential expression from 18. Balan, S. et al. Large-scale human dendritic cell differentiation revealing Seurat. Genes present in 10% of the cells and with a minimal log fold change of Notch-dependent lineage bifurcation and heterogeneity. Cell Rep. 24, 0.25 were selected for differential expression analysis of each cell type in each KO 1902.e6–1915.e6 (2018). line vs their WT counterpart. Only genes with an FDR < 0.05 were considered as 19. Garcia-Alonso, L. et al. Mapping the temporal and spatial dynamics of the significantly differentially expressed. Differential amount of cells produced was calculated using two-sided t test between ZEB2 KO number of cells collected for the human endometrium in vivo and in vitro. Nat. Genet. 53, 1698–1711 (2021). two clones used vs the number of cells collected for all other lines. 20. Vento-Tormo, R. et al. Single-cell reconstruction of the early maternal-fetal interface in humans. Nature 563, 347–353 (2018). 21. Park, J.-E. et al. A cell atlas of human thymic development defines T cell Reporting summary. Further information on research design is available in the Nature repertoire formation. Science 367, eaay3224 (2020). Research Reporting Summary linked to this article. 22. Popescu, D.-M. et al. Decoding human fetal liver haematopoiesis. Nature 574, 365–371 (2019). 23. Tyser, R. C. V. et al. Single-cell transcriptomic characterization of a gastrulating human embryo. Nature 600, 285–289 (2021). Data availability 24. Stuart, T. et al. Comprehensive integration of single-cell data. Cell 177, The single-cell sequencing fastq files are available in the ArrayExpress database (http:// 1888.e21–1902.e21 (2019). www.ebi.ac.uk/arrayexpress) under accession numbers E-MTAB-11623 for single-cell 25. Dou, D. R. et al. Medial HOXA genes demarcate haematopoietic stem cell fate RNAseq and E-MTAB-11616 for single-cell ATACseq. Processed data sets can be during human development. Nat. Cell Biol. 18, 595–606 (2016). queried and downloaded through the web portal www.HiPImmuneatlas.org. Single-cell 26. Jung, H. S. et al. SOX17 integrates HOXA and arterial programs in hemogenic sequencing data was mapped against GRCh38 http://ftp.ensembl.org/pub/release-100/. endothelium to drive definitive lympho-myeloid hematopoiesis. Cell Rep. 34, Publicly available data sets used include Gastrulation http://www.human-gastrula.net/, 108758 (2021). Foetal liver (+kidney +skin) E-MTAB-7407, Foetal thymus E-MTAB-8581, Placenta E- 27. Sharma, A. et al. Onco-fetal reprogramming of endothelial cells drives MTAB-6701, HCC https://doi.org/10.17632/6wmzcskt6k.1, Yolk sac GSE133345 and immunosuppressive macrophages in hepatocellular carcinoma. Cell 183, adult DCs and Macs GSE115006. Source data are provided with this paper. 377.e21–394.e21 (2020). 28. Álvarez-Errico, D., Vento-Tormo, R., Sieweke, M. & Ballestar, E. Epigenetic control of myeloid cell differentiation, identity and function. Nat. Rev. Code availability Immunol. 15, 7–17 (2015). All relevant codes used for data analysis are available from https://github.com/Ventolab/ 29. Bergen, V., Lange, M., Peidli, S., Wolf, F. A. & Theis, F. J. Generalizing RNA iPSCmyeloid. velocity to transient cell states through dynamical modeling. Nat. Biotechnol. 38, 1408–1414 (2020). NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications 15 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 30. Holland, C. H. et al. Robustness and applicability of transcription factor and 59. Hossain, M., Qadri, S. M., Su, Y. & Liu, L. ICAM-1-mediated leukocyte pathway analysis tools on single-cell RNA-seq data. Genome Biol. 21, 36 (2020). adhesion is critical for the activation of endothelial LSP1. Am. J. Physiol. Cell 31. Marcelo, K. L., Goldie, L. C. & Hirschi, K. K. Regulation of endothelial cell Physiol. 304, C895–C904 (2013). differentiation and specification. Circ. Res. 112, 1272–1287 (2013). 60. Farren, M. R. et al. Tumor-induced STAT3 signaling in myeloid cells impairs 32. Stehling-Sun, S., Dade, J., Nutt, S. L., DeKoter, R. P. & Camargo, F. D. dendritic cell generation by decreasing PKCβII abundance. Sci. Signal. 7, ra16 Regulation of lymphoid versus myeloid fate ‘choice’ by the transcription factor (2014). Mef2c. Nat. Immunol. 10, 289–296 (2009). 61. Li, J. et al. The EMT transcription factor Zeb2 controls adult murine 33. Clausen, B. E. et al. Residual MHC class II expression on mature dendritic hematopoietic differentiation by regulating cytokine signaling. Blood 129, cells and activated B cells in RFX5-deficient mice. Immunity 8, 143–155 460–472 (2017). (1998). 62. Scott, C. L. et al. The transcription factor Zeb2 regulates development of 34. Kang, K. et al. Interferon-γ represses M2 gene expression in human conventional and plasmacytoid DCs by repressing Id2. J. Exp. Med. 213, macrophages by disassembling enhancers bound by the transcription factor 897–911 (2016). MAF. Immunity 47, 235.e4–250.e4 (2017). 63. Fardi, M., Mohammadi, A., Baradaran, B. & Safaee, S. ZEB2 knock-down 35. Heinz, S. et al. Simple combinations of lineage-determining transcription induces apoptosis in human myeloid leukemia HL-60 cells. Curr. Gene Ther. factors prime cis-regulatory elements required for macrophage and B cell 21, 149–159 (2021). identities. Mol. Cell 38, 576–589 (2010). 64. Wong, K. L. et al. Gene expression profiling reveals the defining features of the 36. Schneider, C. et al. Induction of the nuclear receptor PPAR-γ by the cytokine classical, intermediate, and nonclassical human monocyte subsets. Blood 118, GM-CSF is critical for the differentiation of fetal monocytes into alveolar e16–e31 (2011). macrophages. Nat. Immunol. 15, 1026–1037 (2014). 65. Grip, O., Bredberg, A., Lindgren, S. & Henriksson, G. Increased 37. Yu, X. et al. The cytokine TGF-β promotes the development and homeostasis subpopulations of CD16 and CD56 blood monocytes in patients with active of alveolar macrophages. Immunity 47, 903.e4–912.e4 (2017). Crohnʼs disease. Inflamm. Bowel Dis. 13, 566–572 (2007). 38. Liu, B. et al. An entropy-based metric for assessing the purity of single cell 66. Rossol, M., Kraus, S., Pierer, M., Baerwald, C. & Wagner, U. The populations. Nat. Commun. 11, 3155 (2020). CD14brightCD16 monocyte subset is expanded in rheumatoid arthritis and 39. Garcia-Alonso, L., Holland, C. H., Ibrahim, M. M., Turei, D. & Saez- promotes expansion of the Th17 cell population. Arthritis Rheum. 64, Rodriguez, J. Benchmark and integration of resources for the estimation of 671–677 (2012). human transcription factor activities. Genome Res. 29, 1363–1375 (2019). 67. Mengos, A. E., Gastineau, D. A. & Gustafson, M. P. The CD14+HLA-DRlo/ 40. Hume, D. A. The many alternative faces of macrophage activation. Front. neg monocyte: an immunosuppressive phenotype that restrains responses to Immunol. 6, 370 (2015). cancer immunotherapy. Front. Immunol. 10, 1147 (2019). 41. Sullivan, A. L. et al. Serum response factor utilizes distinct promoter- and 68. Alban, T. J. et al. Glioblastoma myeloid-derived suppressor cell subsets express enhancer-based mechanisms to regulate cytoskeletal gene expression in differential macrophage migration inhibitory factor receptor profiles that can macrophages. Mol. Cell. Biol. 31, 861–875 (2011). be targeted to reduce immune suppression. Front. Immunol. 11, 1191 (2020). 42. Im, S.-S. et al. Linking lipid metabolism to the innate immune response in 69. Kumar, V., Patel, S., Tcyganov, E. & Gabrilovich, D. I. The nature of myeloid- macrophages through sterol regulatory element binding protein-1a. Cell derived suppressor cells in the tumor microenvironment. Trends Immunol. 37, Metab. 13, 540–549 (2011). 208–220 (2016). 43. Guo, C. et al. Cholesterol homeostatic regulator SCAP-SREBP2 integrates 70. Fleming, V. et al. Targeting myeloid-derived suppressor cells to bypass tumor- NLRP3 inflammasome activation and cholesterol biosynthetic signaling in induced immunosuppression. Front. Immunol. 9, 398 (2018). macrophages. Immunity 49, 842–856.e7 (2018). 71. Knights, A. J. et al. Krüppel-like factor 3 (KLF3) suppresses NF-κB-driven 44. Hamilton, T. A., Zhao, C., Pavicic, P. G. Jr. & Datta, S. Myeloid colony- inflammation in mice. J. Biol. Chem. 295, 6080–6091 (2020). stimulating factors as regulators of macrophage polarization. Front. Immunol. 72. Zhao, Q. et al. Activating transcription factor 3 involved in Pseudomonas 5, 554 (2014). aeruginosa PAO1-induced macrophage senescence. Mol. Immunol. 133, 45. Wang, Y. et al. IL-34 is a tissue-restricted ligand of CSF1R required for the 122–127 (2021). development of Langerhans cells and microglia. Nat. Immunol. 13, 753–760 73. Ricote, M., Li, A. C., Willson, T. M., Kelly, C. J. & Glass, C. K. The peroxisome (2012). proliferator-activated receptor-gamma is a negative regulator of macrophage 46. Ohgidani, M. et al. Direct induction of ramified microglia-like cells from activation. Nature 391, 79–82 (1998). human monocytes: dynamic microglial dysfunction in Nasu-Hakola disease. 74. Ma, W.-T., Gao, F., Gu, K. & Chen, D.-K. The role of monocytes and Sci. Rep. 4, 4957 (2014). macrophages in autoimmune diseases: a comprehensive review. Front. 47. Ji, J. et al. Antagonizing peroxisome proliferator-activated receptor γ facilitates Immunol. 10, 1140 (2019). M1-to-M2 shift of microglia by enhancing autophagy via the LKB1-AMPK 75. Bruveris, F. F. et al. Human yolk sac-like haematopoiesis generates RUNX1-, signaling pathway. Aging Cell 17, e12774 (2018). GFI1- and/or GFI1B-dependent blood and SOX17-positive endothelium. 48. Cabeza-Cabrerizo, M., Cardoso, A., Minutti, C. M., Pereira da Costa, M. & Development 147, dev193037 (2020). Reis e Sousa, C. Dendritic cells revisited. Annu. Rev. Immunol. 39, 131–166 76. Willcockson, M. A. et al. Runx1 promotes murine erythroid progenitor (2021). proliferation and inhibits differentiation by preventing Pu.1 downregulation. 49. Balan, S. & Dalod, M. In vitro generation of human XCR1(+) dendritic cells Proc. Natl Acad. Sci. USA 116, 17841–17847 (2019). from CD34(+) hematopoietic progenitors. Methods Mol. Biol. 1423, 19–37 77. Ludwig, L. S. et al. Transcriptional states and chromatin accessibility (2016). underlying human erythropoiesis. Cell Rep. 27, 3228.e7–3240.e7 (2019). 50. Anselmi, G. et al. Engineered niches support the development of human 78. Ludwig, L. S. et al. Altered translation of GATA1 in Diamond-Blackfan dendritic cells in humanized mice. Nat. Commun. 11, 2054 (2020). anemia. Nat. Med. 20, 748–753 (2014). 51. Tang-Huau, T.-L. et al. Human in vivo-generated monocyte-derived dendritic 79. Bergen, V., Soldatov, R. A., Kharchenko, P. V. & Theis, F. J. RNA velocity— cells and macrophages cross-present antigens through a vacuolar pathway. current challenges and future perspectives. Mol. Syst. Biol. 17, e10282 (2021). Nat. Commun. 9, 2570 (2018). 80. Lavin, Y. et al. Tissue-resident macrophage enhancer landscapes are shaped by 52. Villar, J. & Segura, E. Decoding the heterogeneity of human dendritic cell the local microenvironment. Cell 159, 1312–1326 (2014). subsets. Trends Immunol. 41, 1062–1071 (2020). 81. Gosselin, D. et al. Environment drives selection and function of enhancers 53. Chopin, M. et al. Transcription factor PU.1 promotes conventional dendritic controlling tissue-specific macrophage identities. Cell 159, 1327–1340 (2014). cell identity and function via induction of transcriptional regulator DC- 82. Saeed, S. et al. Epigenetic programming of monocyte-to-macrophage SCRIPT. Immunity 50, 77.e5–90.e5 (2019). differentiation and trained innate immunity. Science 345, 1251086 (2014). 54. Tussiwand, R. et al. Klf4 expression in conventional dendritic cells is required 83. Gosselin, D. et al. An environment-dependent transcriptional network for T helper 2 cell responses. Immunity 42, 916–928 (2015). specifies human microglia identity. Science 356, eaal3222 (2017). 55. Scholz, F. et al. The transcription factor C/EBPβ orchestrates dendritic cell 84. Schulz, C. et al. A lineage of myeloid cells independent of Myb and maturation and functionality under homeostatic and malignant conditions. hematopoietic stem cells. Science 336, 86–90 (2012). Proc. Natl Acad. Sci. USA 117, 26328–26339 (2020). 85. Sakai, M. et al. Liver-derived signals sequentially reprogram myeloid 56. Guilliams, M. et al. Unsupervised high-dimensional analysis aligns dendritic enhancers to initiate and maintain Kupffer cell identity. Immunity 51, cells across tissues and species. Immunity 45, 669–684 (2016). 655–670.e8 (2019). 57. Chiang, M.-C. et al. Differential uptake and cross-presentation of soluble and 86. Ydens, E. et al. Profiling peripheral nerve macrophages reveals two necrotic cell antigen by human DC subsets. Eur. J. Immunol. 46, 329–339 macrophage subsets with distinct localization, transcriptome and response to (2016). injury. Nat. Neurosci. 23, 676–689 (2020). 58. Bassler, K., Schulte-Schrepping, J., Warnat-Herresthal, S., Aschenbrenner, A. C. 87. Bhattacharjee, O., Ayyangar, U., Kurbet, A. S., Ashok, D. & Raghavan, S. & Schultze, J. L. The myeloid cell compartment—cell by cell. Annu. Rev. Unraveling the ECM-Immune Cell Crosstalk In Skin Diseases. Front. Cell Dev. Immunol. 37, 269–293 (2019). Biol. 7, 68 (2019). 16 NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30557-4 ARTICLE 88. Dorrier, C. E. et al. CNS fibroblasts form a fibrotic scar in response to immune and C.S.-S.; C.P.J. performed and analysed the T cell activation experiments; A.J.K. cell infiltration. Nat. Neurosci. 24, 234–244 (2021). designed and performed the ATAC single-cell sample processing; C.A. analysed the data 89. Heaton, H. et al. Souporcell: robust clustering of single-cell RNA-seq data by with contributions from V.L with the trajectories and knock-out analyses, from J.-E.P. genotype without reference genotypes. Nat. Methods 17, 615–620 (2020). with the thymus data set logistic regression and from B.S.W. and D.F.T. with the den- 90. Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: computational identification dritic cells data set; C.A., R.V.-T. and D.A.-E. interpreted the data and wrote the of cell doublets in single-cell transcriptomic data. Cell Syst. 8, 281.e9–291.e9 manuscript with contributions from V.L., R.V.-B., G.T. and A.B.; D.J.G., D.A.-E. and (2019). R.V.-T. supervised the work. All authors read and approved the manuscript. 91. Hafemeister, C. & Satija, R. Normalization and variance stabilization of single- cell RNA-seq data using regularized negative binomial regression. Genome Biol. 20, 296 (2019). Competing interests The authors declare no competing interests. 92. Blondel, V. D., Guillaume, J.-L., Lambiotte, R. & Lefebvre, E. Fast unfolding of communities in large networks. J. Stat. Mech. 2008, P10008 (2008). 93. Cusanovich, D. A. et al. A single-cell atlas of in vivo mammalian chromatin Additional information accessibility. Cell 174, 1309.e18–1324.e18 (2018). Supplementary information The online version contains supplementary material 94. Garcia-Alonso, L. et al. Single-cell roadmap of human gonadal development. available at https://doi.org/10.1038/s41467-022-30557-4. Preprint at Research Square https://doi.org/10.21203/rs.3.rs-496470/v1 (2021). 95. Fornes, O. et al. JASPAR 2020: update of the open-access database of Correspondence and requests for materials should be addressed to Clara Alsinet, Daniel transcription factor binding profiles. Nucleic Acids Res. 48, D87–D92 (2020). J. Gaffney, Damiana Alvarez-Errico or Roser Vento-Tormo. Peer review information Nature Communications thanks Marc Dalod, Nico Lachmann Acknowledgements and Manu Setty for their contribution to the peer review of this work. We thank the Cellular Genetics wet lab support team, Cellular Genetics IT team, Sanger Sequencing operations and Sanger Cytometry Core facility for their essential help. We Reprints and permission information is available at http://www.nature.com/reprints thank the Gene Editing team for providing iPSC knock-out lines. We would also like to thank Ruxandra Tesloianu and Luz Garcia-Alonso for their help setting up the scATAC- Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in seq computational analysis. We thank Jana Eliasova for her help with figure design and published maps and institutional affiliations. Christina Usher and Aidan Maartens for their edits in the text. This work was mainly funded by the Open Targets consortium (OTAR026 and OTAR032 project) and the Wellcome Sanger core funding (WT206194) with additional support from Open Targets Open Access This article is licensed under a Creative Commons projects OTAR037, OTAR2065, OTAR2071. The authors are grateful to the funders for Attribution 4.0 International License, which permits use, sharing, their support and additional care given to their members during the COVID-19 pan- adaptation, distribution and reproduction in any medium or format, as long as you give demic. D.A.-E. thanks CERCA Programme/Generalitat de Catalunya and the Josep appropriate credit to the original author(s) and the source, provide a link to the Creative Carreras Foundation for institutional support. This study makes use of cell lines and data Commons license, and indicate if changes were made. The images or other third party generated by the HiPSci Consortium, funded by The Wellcome Trust and the MRC material in this article are included in the article’s Creative Commons license, unless (Medical Research Council). For the purpose of Open Access, the author has applied a indicated otherwise in a credit line to the material. If material is not included in the CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission. This publication is part of the Human Cell Atlas—www.humancellatlas. article’s Creative Commons license and your intended use is not permitted by statutory org/publications. regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. Author contributions C.A., D.J.G., D.A.-E. and R.V.-T. conceived and designed the experiments and analyses; C.A. and M.N.P. performed the experiments with contributions from E.B., I.K., R.V.-B., © The Author(s) 2022 NATURE COMMUNICATIONS | (2022)13:2885 | https://doi.org/10.1038/s41467-022-30557-4 | www.nature.com/naturecommunications 17