Single-cell RNA sequencing reveals developmental heterogeneity among early lymphoid progenitors

Abstract

Single-cell RNA sequencing is a powerful technology for assessing heterogeneity within defined cell populations. Here, we describe the heterogeneity of a B220⁺CD117^intCD19^-NK1.1^- uncommitted hematopoietic progenitor having combined lymphoid and myeloid potential. Phenotypic and functional assays revealed four subpopulations within the progenitor with distinct lineage developmental potentials. Among them, the Ly6D⁺SiglecH^-CD11c^- fraction was lymphoid-restricted exhibiting strong B-cell potential, whereas the Ly6D^-SiglecH^-CD11c^- fraction showed mixed lympho-myeloid potential. Single-cell RNA sequencing of these subsets revealed that the latter population comprised a mixture of cells with distinct lymphoid and myeloid transcriptional signatures and identified a subgroup as the potential precursor of Ly6D⁺SiglecH^-CD11c^- Subsequent functional assays confirmed that B220⁺CD117^intCD19^-NK1.1^- single cells are, with rare exceptions, not bipotent for lymphoid and myeloid lineages. A B-cell priming gradient was observed within the Ly6D⁺SiglecH^-CD11c^- subset and we propose a herein newly identified subgroup as the direct precursor of the first B-cell committed stage. Therefore, the apparent multipotency of B220⁺CD117^intCD19^-NK1.1^- progenitors results from underlying heterogeneity at the single-cell level and highlights the validity of single-cell transcriptomics for resolving cellular heterogeneity and developmental relationships among hematopoietic progenitors.

Keywords: hematopoiesis; heterogeneity; lineage priming; multipotentiality; single‐cell RNA sequencing.

Figure 1. Phenotypic and functional heterogeneity of the EPLM progenitor population (B220⁺CD117^intCD19⁻NK1.1⁻) based on expression of Ly6D, SiglecH, and CD11c

1. A

   Representative FACS plots of EPLM from the bone marrow of WT mice with the addition of Ly6D, SiglecH, and CD11c identifying four subpopulations.

2. B, C

   Percentages (B) and absolute numbers (C) of WT EPLM subpopulations (n = 5).

3. D–F

   In vitro limiting dilution analysis of Ly6D⁺, TN, SiglecH⁺, and CD11c⁺ for B‐cell (D), T‐cell (E), or myeloid (F) potentials.

4. G, H

   B‐cell reconstitution of sub‐lethally irradiated Rag2‐deficient mice with 4 × 10³ Ly6D⁺ (n = 5) or TN (n = 4) cells from WT. (G) Representative FACS plots from spleens of reconstituted mice. (H) Quantification of CD19⁺IgM⁺ splenocytes.

Data information: In (B, C and H), data are presented as mean ± SEM. (D–F) Independent repetitions for each experiment are provided in a table (Appendix Fig S2B). (G) MZB: marginal zone B cells; FB: follicular B cells.

Figure 2. Identification and analysis of Flt3Ltg EPLM subpopulations

1. A

   Representative FACS plots of EPLM from the bone marrow of Flt3Ltg mice with the addition of Ly6D, SiglecH, and CD11c identifying four subpopulations.

2. B, C

   Comparison of EPLM subpopulations from WT (n = 5, circles) and Flt3Ltg (n = 5, squares) mice in percentages (B) or absolute numbers (C).

3. D–F

   B‐cell, T‐cell, and myeloid precursor frequencies of EPLM subpopulations from Flt3Ltg mice obtained by limiting dilution performed as in Fig 1.

Data information: In (B and C), data are presented as mean ± SEM. (C) Two‐tailed unpaired Student's t‐tests with P‐values = 8.6 × 10⁻³ (Ly6D⁺), 0.022 (SiglecH⁺), 1.4 × 10⁻³ (TN), and 2 × 10⁻³ (CD11c⁺). (D–F) Independent repetitions for each experiment are provided in a table (Appendix Fig S2B).

Figure 3. Ly6D⁺ and TN EPLM subpopulations are transcriptionally and developmentally related

1. A

   Gene expression ratio of Ly6D⁺ versus TN cells from Flt3Ltg (vertical axis) plotted against the average expression intensity (horizontal axis), showing 1,008 differentially expressed genes (DEG, stars): 493 highly expressed in Ly6D⁺ (Up, red) and 515 highly expressed in TN (Down, blue) identified by bulk RNA‐seq (n = 4).

2. B

   PCA of 152 Ly6D⁺ and 213 TN single cells from Flt3Ltg using as gene set the DEG identified in (A) and colored according to the cell type. Average gene expression was centered to zero.

3. C, D

   Kinetics of CD19⁺ and Ly6D⁺ EPLM generation in vitro. (C) Ly6D/CD19 FACS plots of the in vitro progeny of Ly6D⁺ (upper row) and TN EPLM (lower row) at days 1–3 after initiation of culture. Cells shown are SiglecH⁻CD11c⁻NK1.1⁻. (D) Kinetics of Ly6D⁺ EPLM and CD19⁺ cell generation in vitro from Ly6D⁺ (top graph) and TN EPLM (bottom graph).

4. E

   Heatmap with pairwise Pearson's transcriptome correlation of Ly6D⁺, TN, and pro‐B averaged populations (bulk RNA‐seq, n = 4).

Data information: (A) Dashed horizontal lines: DEG threshold (abs|log₂(FoldChange)| > 1 and FDR < 0.05). (C, D) A representative experiment is shown (n = 3).

Figure 4. Cell clustering identifies three TN and two Ly6D⁺ subgroups with distinct genetic signatures

1. A

   PCA generated as in Fig 3B showing the subgroups revealed by PAM clustering method (see Materials and Methods). Circles: Ly6D⁺, triangles: TN. G1 Ly6D⁺ (n = 56), G2 Ly6D⁺ (n = 82), G3 TN (n = 85), G4 TN (n = 52), G5 TN (n = 56).

2. B

   Heatmap with pairwise Pearson's transcriptome correlation of Ly6D⁺ and TN subgroups. The number of cells (n) per subgroup is specified in (A).

3. C–G

   Violin plots with genes highly expressed in G1 Ly6D⁺ (C), G4 TN (D), G5 TN (E), (G1, G2) Ly6D⁺ and G3 TN (F), or G4 and G5 TN (G) subgroups compared to the other subgroups. The median expression level is shown with a line when more than 50% of the cells express the indicated gene.

Figure 5. Mixed‐ and opposed‐lineage states at the single‐cell level

1. A

   Same PCA plot as in Fig 4A summarizing the genetic signatures of the Ly6D⁺ and TN subgroups revealed by our in silico analysis.

2. B–D

   Scatter plots showing the expression levels in log₂FPKM of selected B and T (B), neutrophil (Neu) and monocyte/macrophages (Mo/Mc) (C) or myeloid (Mye) and lymphoid (Lym) (D) lineage marker pairs in the Ly6D⁺ and TN subpopulations. Dotted vertical and horizontal lines delimit when the transcript of the indicated gene is detected (> 0). Percentages within the double‐positive area of the plot indicate the fraction of cells co‐expressing both genes to the number of cells expressing only one gene (top: gene on vertical axis; bottom: gene on horizontal axis).

3. E

   B‐cell, myeloid, and bipotent (B/Mye) developmental potential of the indicated single‐cell sorted populations from Flt3Ltg. Three independent experiments were performed for Ly6D⁺ and TN cells and one for the control pro‐B and GMP cells. Shown is mean ± SEM.