Transcriptomic landscape of circulating mononuclear phagocytes in Langerhans cell histiocytosis at the single-cell level

Langerhans cell histiocytosis (LCH) is a clinically heterogeneous myeloid neoplasm that is most common in children but can occur at any age and may be life-threatening in its most severe forms.¹ This pathology is underpinned by alterations in the MAPK pathway, with 50% to 60% of patients exhibiting the BRAF^V600E mutation in lesions.^2,3 Notably, LCH lesions are highly inflammatory, containing a variety of different cytokines such as tumor necrosis factor and SPP1,^4-6 with high levels of proinflammatory cytokines often also detected in the blood of LCH patients.^7-9 To date, front-line therapy for LCH patients with multisystem or multifocal disease has been standard oncologic chemotherapy (vinblastine and prednisone)¹⁰; however, the identification of BRAF^V600E has paved the way for targeted therapies, with preliminary evidence showing promising responses to BRAF inhibitors in adults and children with LCH.^11-16 

Importantly, recent studies show that BRAF^V600E also serves as a biological marker to trace the origin of pathologic LCH cells and is most commonly detected in circulating CD14⁺ monocytes, CD16⁺ monocytes, and CD1c⁺ conventional dendritic cell (cDC2), as well as bone marrow CD34⁺ progenitor cells.³ These findings indicated the possibility that circulating myeloid cells might be directly involved in disease pathogenesis. In line, in vitro studies have now shown that both CD1c⁺ cDC2 and CD14⁺ classical monocytes can be driven to acquire high expression of langerin (CD207) and CD1a after culture with granulocyte-macrophage colony-stimulating factor, transforming growth factor, and additional notch ligation, suggesting their potential to be candidate precursors in the circulation and to contribute to LCH lesions.^3,17,18 What is currently unclear is how the broader phenotypic and molecular landscape of the blood myeloid compartment in LCH is different in patients and what the developmental and functional consequences of BRAF-mediated MAPK pathway deregulation are on the mononuclear phagocyte compartment.

In this study, we focused on mononuclear myeloid cells in LCH and adopted a multilayered approach combining immune-phenotyping and clinical observations with single-cell transcriptomics, unveiling how they may be affected in newly diagnosed patients and in response to BRAF inhibition.

We initially recruited 233 consecutive pediatric patients with newly diagnosed LCH (age <18 years) admitted to Beijing Children’s Hospital from May 2018 to December 2019. Of these, 217 patients who had available plasma samples and were tested for cell-free BRAF^V600E (cfBRAF^V600E) mutations were analyzed in this study (Table 1; and supplemental Table 1 available on the Blood Web site). Forty-nine patient blood samples were selected at random for flow cytometric analysis (supplemental Table 2). We also analyzed samples from 22 children with relapsed or refractory LCH bearing the BRAF^V600E mutation, who were treated with the BRAF inhibitor dabrafenib (http://www.chictr.org.cn, identifier: ChiCTR2000032844; supplemental Table 3). For the detailed sample information please refer to supplemental Methods. This study was approved by the Beijing Children’s Hospital Institutional Review Board (IEC-C-008-A08-2018-75 and IEC-C-028-A10-2016-K-31) and was conducted in accordance with The Declaration of Helsinki. Informed consents were obtained from guardians of the patients.

Healthy control and LCH skin samples were digested with 1.2 U/mL dispase II (Roche, Indianapolis, IN) and 0.25% trypsin-EDTA (Gibco) containing 0.25 mg/mL DNase I (DN25; Sigma), to obtain single-cell suspensions, as previously published.¹⁹ Blood samples were treated with lysis buffer (Becton Dickinson) to remove erythrocytes and were then labeled with antibodies. For the detailed information of antibody labeling for flow cytometry, please refer to supplemental Methods.

Single-cell capture and downstream library construction following a modified single-cell tagged reverse transcription protocol were achieved as previously reported.^19-22 Sequencing was performed on the Illumina Hiseq X Ten platform (Novogene). scRNA-seq data analysis is described in the supplemental Methods.

Survival rates were analyzed by the Kaplan-Meier method, and subgroups were compared with the log-rank test. Tests were performed using SPSS 25.0 software (SPSS Inc, Chicago, IL). Wilcoxon-rank-sum test in R version 3.6.2 was used for comparisons of cell proportions based on flow cytometry analysis, cfBRAF^V600E levels between any 2 groups, and gene expression levels between 2 cell clusters. Kruskal-Wallis test was used for significance verification. We considered P < .05 statistically significant.

See supplemental Methods for additional description of methods.

We enrolled 217 newly diagnosed pediatric LCH patients during the period of May 2018 to December 2019 (Table 1; supplemental Table 1). The median age at disease diagnosis was 3.1 years (range, 0.07-13.60), and the ratio of males to females was 1.44:1. Clinical classification placed 119 patients (54.8%) in the single system (SS) low-risk group, 65 (30%) in the multisystem (MS) low-risk group, and 33 (15.2%) in the MS high-risk group. Forty-nine blood samples were randomly selected, and after confirming that they were broadly representative of the entire newly diagnosed cohort, they were subjected to detailed flow cytometric analysis (supplemental Table 2). In addition, 14 blood and 4 skin samples were used for scRNAseq (Figure 1A). The BRAF^V600E status was determined in lesions of 174 LCH patients, and 95 were found to harbor the BRAF^V600E mutation, among which 52 patients also had cfBRAF^V600E detectable in their plasma at diagnosis (supplemental Figure 1A-B). Of the 79 patients who tested negative for BRAF^V600E in tissue biopsy, 77 also did not present with cfBRAF^V600E in plasma (supplemental Figure 1A-B). Here, the plasma cfBRAF^V600E was detected in just 87.5% of high-risk patients with tissue mutation, similar to our recent result,²³ whereas the frequency was 100% in previous studies.^24,25 This difference may result from the varying cohort size and ethnic diversity.

The presence of cfBRAF^V600E in LCH patient plasma is an adverse prognostic indicator and has been associated with more severe clinical manifestations of the disease,^23,24,26 but it is unclear how this might be linked to changes in the circulating mononuclear cell compartment. We allocated each patient within the clinical groupings to different cfBRAF^V600E groups: cfBRAF^V600E negative (cfBRAF^V600E [−]; cfBRAF^V600E ddPCR levels <0.1%, n = 154), cfBRAF^V600E low mutation (0.1% to 3.1%, n = 47), and cfBRAF^V600E high mutation (≥3.1%, n = 16; supplemental Figure 1C-D). When we compared the mutation load and clinical risk groups, a significant association between the level of cfBRAF^V600E and clinical classification of LCH in newly diagnosed children was detected (Figure 1B).

BRAF^V600E mutation in LCH is commonly detected in circulating monocytes and dendritic cells (DCs),³ and aberrant activation of the MAPK pathway can affect myeloid cell differentiation in mice^27,28; however, the influence of BRAF^V600E mutation dysregulated MAPK activation on the proportion of circulating mononuclear phagocytes in LCH patients remains unknown. First, the proportions of CD14⁺ monocytes, CD11c⁺BDCA2⁻CD14⁻ cells (containing conventional DCs [cDCs] and CD16⁺ monocytes) and BDCA2⁺CD11c^–CD14⁻ pDCs in peripheral blood (PB) of newly diagnosed LCH patients (n = 49) were analyzed in the different clinical groups and healthy controls by flow cytometry (Figure 1C). This revealed higher frequencies of CD14⁺ monocytes in all clinical groups of LCH patients compared with controls but no significant difference in frequency among the subgroups of LCH (Figure 1D). On the contrary, the proportions of CD11c⁺BDCA2⁻CD14⁻ cells (containing cDCs and CD16⁺ monocytes) were significantly lower in LCH patients than controls, again with no significant difference between LCH subgroups. Intriguingly, although pDC frequencies were lower in MS patients than controls, a further significant decrease in MS high-risk patients compared with any of the other groups was detected (Figure 1D). The same patterns were present when patients were categorized according to cfBRAF^V600E level (Figure 1E). Overall, pediatric LCH is in general related to significantly higher relative frequencies of CD14⁺ monocytes, whereas lower pDC frequencies are in addition significantly associated with disease severity and high mutational load.

Immune-phenotyping analysis provided preliminary evidence that monocyte and DC subsets in blood are affected in proportions in LCH. Therefore, we next investigated the molecular alterations occurring in these cells by using fluorescence-activated cell sorting for their enrichment and further subjecting them to a modified version of single-cell tagged reverse transcription-seq.^19-22 To this end, 14 patients were selected from the newly diagnosed cohort that together represented all clinical subgroups, MS high risk (n = 6), MS low risk (n = 4), and SS low risk (n = 4), and exhibited a range of cfBRAF^V600E levels (Figure 2A). A total of 1988 CD45⁺Lin⁻ (CD3, CD19, and CD56) HLA-DR⁺ cells were isolated from the LCH peripheral blood samples (Figure 2B; supplemental Figure 2A-B), as well as 739 CD45⁺Lin⁻HLA-DR⁺ cells from the blood of 12 healthy age-matched children as controls. Finally, unsupervised clustering analysis was applied to these data, and 11 cell clusters were identified and annotated according to their expression of a curated list of known feature genes (Figure 2C; supplemental Figure 2C-H).

We identified CD14⁺ classical monocytes (CD14⁺ Mo), CD16⁺⁺ nonclassical monocytes (CD16⁺⁺ Mo), and CD14⁺CD16⁺ intermediate monocytes (iMo) based on the expression of CD14 and FCGR3A (CD16). We also identified a minor cDC1 cluster expressing CLEC9A, XCR1, and CADM1 and a major cDC2 cluster expressing CD1c and FCER1A (Figure 2C; supplemental Figure 2E). Furthermore, the recently described proinflammatory DC3 subset (exhibiting higher levels of LMNA, CD14, CD163, and CDKN1A) was also detected within cDC2 cluster,^29,30 yet being not fully discriminated from cDC2, whose proportion was significantly higher in cDCs from LCH patients than that from healthy controls at the transcriptomic level (supplemental Figure 2G). On the other hand, pDCs highly expressed CLEC4C and NRP1 (Figure 2C). According to the expression of canonical marker genes, we also identified B cells, plasma cells, and immature progenitors (Figure 2B-C).

Having identified the link between plasma level of cfBRAF^V600E and clinical severity in our cohort, we went on to investigate the transcriptomes of the CD1c⁺ cDCs (cDC2/DC3) and 3 monocyte subpopulation clusters (CD14⁺, CD16⁺⁺, and iMo) in each cell subset individually by comparing each mutational load group (n = 4 for high- and low-cfBRAF^V600E, and n = 6 for cfBRAF^V600E negative) with the cells from healthy control blood samples serving as reference. Each of the predominant monocyte and DC clusters were present in blood from every patient and control tested, although at variable proportions (supplemental Figure 2I). Examining the top 5 most highly DEGs of these 4 clusters in each mutational group, we found that CD16⁺ Mo and cDC2/DC3 of high-mutation patients highly expressed C1QA, C1QC, or C1QB, which encode the polypeptide of serum complement subcomponent C1q (Figure 2D-E). In CD14⁺ Mo of high-mutation patients, we observed high expression of CLU, which is involved in neurodegenerative diseases and cancerous processes.³¹ Recently, the examination of brain biopsy specimens from LCH neurodegenerative disease patients reported the perivascular infiltration by cells with monocytic phenotype (CD14⁺CD33⁺CD163⁺P2RY12⁻) carrying the BRAF^V600E mutation.³² Future studies will be needed to assess whether brain-infiltrating monocytes in LCH correspond to the CLU-expressing monocytes detected here in the high-mutation group.

Subsequently, we focused on identifying those DEGs from the different cfBRAF^V600E mutational groups that were involved in the RAS-MAPK-ERK pathway. Overall, all cell populations from each of the mutational groups expressed significantly higher levels of genes related to the RAS-MAPK-ERK pathway than the corresponding cell populations from healthy controls (Figure 3A; supplemental Figure 3A). Using single-sample gene set enrichment analysis (ssGSEA), we observed higher RAS-MAPK-ERK related signaling activities in CD14⁺ Mo, CD16⁺⁺ Mo, iMo, and cDC2/DC3 from the high-mutation group than those from the other 2 mutation groups (Figure 3B). Thus, in our cohort of newly diagnosed LCH patients, almost all individuals tested exhibited diverse aberrations in expression levels of RAS-MAPK-ERK pathway-related genes across the circulating myeloid compartment with the DEGs varying by cfBRAF^V600E mutational group and cell subset. Interestingly, the expression of genes involved in RAS and MAPK-related pathways was enriched in natural killer/T cells of LCH patients compared with healthy counterpart but not in B cells or plasma cells of LCH patients (supplemental Figure 3B-C).

We then focused on the high-mutation patients and used cytoscape³³ to visualize the molecular interaction network between the RAS-MAPK-ERK related genes and transcription factors in the identified DEGs in each cell population. We found evidence that the different cell types exhibited differential regulation of the networks around RAS-MAPK-ERK: for example, in CD14⁺ Mo, SOS1, a gene involved in controlling the activation of the RAS/MAPK signaling pathway, was related to CEBPD, but the later one was linked with JAK3 in CD16⁺⁺ Mo (Figure 3C). Together these observations provide evidence for widespread and differential MAPK-associated pathway activation and inflammation in circulating mononuclear phagocytes before their migration to the tissues.

Next, effects of targeted treatment of the RAS-MAPK-ERK pathway using a specific BRAF enzyme inhibitor dabrafenib were assessed at the clinical and transcriptional level. A cohort of 22 children with LCH who were positive for the BRAF^V600E mutation had refractory disease or progressive organ dysfunction and received dabrafenib treatment was investigated. After 1 month of treatment with dabrafenib, 19 patients showed disease regression, with 3 other patients in a stable condition, accompanied by improvement of laboratory tests and imaging examinations (Figure 4A; supplemental Figure 4A-B). After 3 months of targeted therapy, just 1 patient showed evidence of disease progression (Figure 4A). One- and two-year progression-free survival rates were 63.9 ± 12.2% and 47.9 ± 16.6%, respectively (supplemental Figure 4C). Plasma analyses of cfBRAF^V600E revealed a significant decrease of cfBRAF^V600E load after 1 month of treatment, consistent with clinical responses (Figure 4B).

To investigate the effects of dabrafenib treatment on the circulating myeloid compartment, we performed scRNA-seq on PB CD45⁺Lin⁻HLA-DR⁺ cells from 5 children within the targeted therapy cohort before treatment and after 1 month of therapy (Figure 1A; supplemental Figure 4D). A total of 1261 cells passed quality control and were used in downstream analysis (supplemental Figure 4E). Using the same curated list of feature genes as above (Figure 2C; supplemental Figure 2E,G), we similarly annotated the 10 clusters (Figure 4C; supplemental Figure 4F). Focusing on these subsets, significant downregulation of a series of genes related to the RAS-MAPK-ERK pathway following treatment was uncovered, including FBP1 and PSMA2 specifically in CD14⁺ Mo, and KRAS, RASA2, and CD81 in CD16⁺⁺ Mo (Figure 4D). CD14⁺ Mo showed posttreatment downregulation of CD63, involved in tumor progression,³⁴ along with the antiapoptosis gene BCL2A1,³⁵ whereas expression of the cancer suppression gene EGR1³⁶ was elevated (Figure 4D). In CD16⁺⁺ Mo, the expression level of C1QB, a characteristic gene previously identified in patients with high levels of BRAF^V600E mutation, decreased significantly (Figure 4D). Conducting the same analysis across the other cell subsets revealed several interesting phenomena: in cDC2/DC3, the expression of AREG, which is related to activation of the MAPK pathway,³⁷ and IFITM1, which is implicated in metastasis in colorectal cancer,³⁸ were significantly downregulated after 1 month of treatment (supplemental Figure 4G). Although certain changes appeared to be cell type specific, the expression of IGLL5, encoding one of the immunoglobulin λ-like polypeptides and considered a prognostic marker in several cancers,^39,40 was generally upregulated across mononuclear phagocytes clusters after treatment (Figure 4D; supplemental Figure 4G).

Using GSEA,⁴¹ we found that expression of genes involved in multiple signaling pathways was specifically enriched in pretherapy CD14⁺ Mo compared with posttherapy, including pathways related to inflammation (interferon response, interleukin-1 family and complement), cellular metabolism (cellular response to hypoxia and glycolysis), and RAS-MAPK signaling (MAPK6/MAPK4 and regulation of RAS by GTPase-activating proteins (GAPs); Figure 4E; supplemental Figure 4H-I; supplemental Table 13). Similarly, expression of genes involved in inflammation and MAPK-related pathways was enriched in the cDC2/DC3, CD16⁺⁺ Mo, and pDCs in pretherapy samples compared with posttherapy samples; for instance, oncogenic-MAPK signaling, reported to be intricately linked to tumorigenesis, was enriched in pretherapy CD16⁺⁺ Mo, whereas the genes related to cellular of responses to stress including activation of cell cycle stagnation, apoptosis or senescence were enriched in cDC2/DC3 in pretherapy sample (Figure 4E).⁴² 

Taken together, this analysis points toward inhibition of the MAPK pathway, prevention of inflammation, and regulation of cell metabolism as key mechanisms underlying the efficacy of targeted therapy with BRAF inhibitor.

Skin is one of the most frequently affected organs in LCH,⁴³ with a variety of cutaneous manifestations, including widespread multiple rashes. As expected, immunohistochemistry analysis revealed numerous CD207⁺CD1a⁺ cells distributed in the epidermis and dermis (supplemental Figure 5A). To understand the heterogeneity among different clusters expressing CD207 and CD1a, we sorted CD45⁺ cells from LCH skin lesions and CD207⁺CD1a⁺ LCs or CD45⁺ cells from healthy epidermis as controls using flow cytometry (supplemental Figure 5B-C) and then performed scRNA-seq on these populations. After quality control and removing nonhematopoietic cells (supplemental Figure 5D-F), our analysis revealed 7 clusters in LCH skin and 3 clusters in normal skin, respectively, which were then annotated based on expression of known feature genes (Figure 5A-C; supplemental Figure 5G).

We found that the major population of CD207⁺CD1a⁺ cells in LCH skin contained 4 heterogeneous clusters, which were also readily discriminated when the dataset was integrated with that from normal skin (supplemental Figure 5H). The LC_proliferating cluster (LC_Prolif.) uniquely expressed proliferation-related genes, including MKI67 and CENPF (Figure 5A-C). Notably, SPP1, which is associated with T-cell activation and recruitment and reported specifically upregulated in LCH cells,⁵ and CCR6,⁴⁴ expressed by pathologic dendritic cells and may contribute to their accumulation in nonlymphoid organs such as skin and bone, had the lowest expression in LC_control (Figure 5B-C). Four heterogeneous LC clusters in lesional skin had partially overlapping patterns of gene expression; for example, they highly expressed CD58, whereas T cells expressed CD2, both genes involved in T cell activation (supplemental Figure 5I).⁴⁵ We also identified other cell populations, including macrophages (CD14 and CD68), neutrophils (FCGR3B and VNN2), pDCs (CLEC4C and NRP1), and T cells (CD3E and CD2; Figure 5B-C; supplemental Figure 5G).

When focusing on the expression of RAS-MAPK-ERK cascade related genes, we found that MAPK14, MAPK9, RRAS, and JAK2 were significantly upregulated in LC_1 compared with LC_control, LC_0, and LC_2 (Figure 5D). Consistently, expression of gene sets associated with the RAS-MAPK-ERK related pathway was significantly enriched in LC_1 in comparison with LC_control, LC_0 or LC_2 (Figure 5E), assuming that different cues may contribute to heterogeneities across LC subclusters, such as differentiation hierarchy,⁴⁶ mutation status,⁴⁷ and lesion microenvironment. In addition, LC_1 showed enriched expression of genes involved in various inflammation-related and collagen degradation–related pathways compared with LC_control (supplemental Figure 5I). These results, especially in LC_1 cluster, indicate the inflammatory myeloid neoplastic nature of LCH cells.

Having performed correlative analysis, we found that the correlation coefficient between LC_1 from LCH skin and cDC2/DC3 from LCH PB was higher than that of other circulating cell clusters at single-cell resolution, supported by integrated analysis with published LCH lesion 10× dataset (Figure 5F; supplemental Figure 5J).

Blood mononuclear phagocytes in healthy humans have been extensively described by RNA sequencing and mass spectrometry,^29,48 but in rare diseases, such as pediatric LCH, relatively little is known of their altered biology. In this study, we adopted a multilayered approach combining not only traditional immune-phenotyping and clinical observations, but also single-cell transcriptomics to propose that the oncogenic BRAF^V600E mutation may lead to widespread MAPK activation in circulating mononuclear phagocytes. Importantly, this set of experiments also revealed the possible mechanisms underlying the therapeutic effects of BRAF inhibitor.

We uncovered widespread alterations to expression of genes related to RAS-MAPK-ERK across all patient groups and myeloid mononuclear phagocyte subsets, with the highest RAS-MAPK-ERK pathway activity observed in patients with high mutational load. This provides an additional important perspective to LCH immunopathogenesis, suggesting that the oncogenic BRAF^V600E mutation may lead to activation of RAS-MAPK-ERK pathway not only in LCH lesions, but also in circulating mononuclear phagocytes, the presumed precursors of tissue LCH cells. In addition, activation of MAPK pathway was also detected in cells from lineages that do not commonly harbor the MAPK activating mutation,^3,25 such as T cells and natural killer cells, highlighting the possible contribution of cell-extrinsic MAPK activation. Future studies aiming to assess simultaneously both the transcriptomic profiles and the BRAF^V600E mutation status in individual cells will be needed to understand the extent to which the MAPK activation is mediated by the cell-intrinsic or -extrinsic effects, such as hypercytokinemia.

We next asked how mononuclear phagocytes would respond to the BRAF inhibitor dabrafenib. Previous reports on the mechanism underlying BRAF inhibitors are mainly in the studies of melanoma. In addition to its molecularly targeted inhibition, BRAF inhibitor has immunomodulatory effects.⁴⁹ On one hand, some of the immunosuppressive effect signals related to the BRAF-mutant tumor microenvironment, such as interleukin-1 and interleukin-6, decreased in posttherapy mononuclear phagocytes.^50,51 On the other hand, the activities of multiple proinflammatory pathways including interferon and chemokine signaling significantly decreased to levels of healthy individuals in posttherapy mononuclear phagocytes.⁴⁹ Our study suggests that dabrafenib may adopt a similar inflammatory regulatory mechanism in LCH. Of note, metabolic rewiring in posttherapy compared with pretherapy samples was evident in all mononuclear phagocytes subsets investigated here. Specifically, the signaling of glycolysis and pentose phosphate pathway, which are essential to the proliferation and survival of cells bearing BRAF^V600E,^52,53 are significantly downregulated by dabrafenib, suggesting a previously unknown mechanism of this targeted therapy. Moreover, the hypoxia status of circulating cells seems to be improved in the meanwhile. These findings imply that dabrafenib has the potential to modulate metabolic status in addition to its expected role in immunomodulation. Whether and to what extent metabolic intervention mediates the effect of BRAF inhibitor and its role in the treatment of LCH remains to be further investigated.

In line with previous reports in adults and children,^3,54 higher frequencies of CD14⁺ Mo compared with healthy controls were also observed in LCH patients. In contrast, pDCs were markedly and significantly less frequent in all LCH subgroups, with the lowest frequencies detected in MS high-risk and high-mutation patients, highlighting pDC frequency as a promising cellular biomarker to be further developed for assessing disease severity and for the evaluation of treatment efficacy by the current standard of care or novel treatment options. Constitutively active MAPK pathway induced by BRAF^V600E may lead to a myeloid differentiation bias, mirrored by the higher proportions of CD14⁺ monocytes and lower proportions of pDC in the HLA-DR⁺ population.^54-56 Another explanation could be their influx to the inflamed tissue, because Halbritter et al⁴⁶ and we identified this cell cluster in LCH lesions.

One caveat of this study is that the relatively low number of sequenced cells may complicate some comparisons, although most comparisons are of biological significance because flow cytometry enrichment of concerned cell populations was carried out, and most clusters in analysis contained sufficient cells. Integration of multisourced single-cell omics derived from LCH patients in the future would undoubtedly further improve the understanding of the pathogenic and therapeutic mechanisms of the disease.

The authors thank all patients and their families who gave them permission to include samples and clinical data in this study.

This study was supported by grants from the National Key Research and Development Program of China, Stem Cell and Translational Research (2017YFA0103401, 2016YFA0100601, and 2020YFA0112402), the National Natural Science Foundation of China (31425012, 81890991, 31871173, 81800102, 82070202, and 81600077), the Program for Guangdong Introducing Innovative and Entrepreneurial Teams (2017ZT07S347), the Key Research and Development Program of Guangdong Province (2019B020234002), and the China Postdoctoral Science Foundation (2018M643373). F.G. is supported by Singapore Immunology Network (SIgN) core funding. F.G. is a European Molecular Biology Organization (EMBO) YIP awardee and is supported by the Singapore National Research Foundation Senior Investigatorship (NRFI) NRF2016NRF-NRFI001-02. C.Z.W.L. is supported by the A*STAR Graduate Scholarship.

Contribution: B.L., R.Z., F.G., and Y.L. designed the study; H.S. and Z. Bian performed the sample preparation and the fluorescence-activated cell sorting with help from Y.Z., Y.N., and Z. Bai; T.H. and H.S. performed the single-cell RNA sequencing with help from J.Z., C.L., and X.L.; L.C., Z.L., N.L., H.M., and D.W. collected and prepared the samples; H.H. and T.H. performed the bioinformatics analysis with help from J.H., Y.G., Z.L., E.K., H.S., Z. Bian, and C.Z.W.L.; and H.S., H.H, L.C., E.K., Z. Bian, T.H., C.Z.W.L., Y.L., R.Z., F.G., and B.L. wrote the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interest.

Correspondence: Yu Lan, Key Laboratory for Regenerative Medicine of Ministry of Education, Institute of Hematology, School of Medicine, Jinan University, Guangzhou, Guangdong, 510632, China; e-mail: [email protected]; Florent Ginhoux, Singapore Immunology Network (SIgN), A(∗)STAR, 8A Biomedical Grove, Immunos Building, Singapore 138648, Singapore; e-mail: [email protected]; Rui Zhang, Laboratory of Hematologic Diseases, Beijing Pediatric Research Institute, Beijing Children's Hospital, Capital Medical University, National Center for Children’s Health, Beijing, 10045, China; e-mail: [email protected]; and Bing Liu, State Key Laboratory of Experimental Hematology, Institute of Hematology, Fifth Medical Center of Chinese PLA General Hospital, Beijing, 100071, China; e-mail: [email protected].