Dissection of progenitor compartments resolves developmental trajectories in B-lymphopoiesis

Despite the fact that B cell development represents one of the most carefully characterized developmental pathways in the hematopoietic system, the identity of the earliest committed B cell progenitor remains obscure. Even though the B220⁺CD43⁺ pre-pro–B or Fraction A (FrA) cells (Hardy et al., 1991; Li et al., 1993) contain a substantial fraction of early B-lineage progenitors, it constitutes a heterogeneous population of cells with varying lineage potentials. Despite the development of more advanced isolation strategies (Sen et al., 1990; Rolink et al., 1994; Li et al., 1996; Tudor et al., 2000), a large portion of the cells in the B220⁺CD19⁻ FrA subpopulations retain T-lineage potential (Martin et al., 2003; Rumfelt et al., 2006; Mansson et al., 2010) as well as the ability to generate myeloid cells (Alberti-Servera et al., 2017). The difficulty in identifying CD19-negative lineage committed B cell progenitors indicated that B-lineage cell fate is associated with CD19 expression (Rumfelt et al., 2006). This would be in line with the fact that CD19 is a direct target gene for the transcription factor (TF) PAX5, which forms a regulatory network with EBF1, FOXO1, and TCF3 to establish stable B-lineage commitment (Nutt et al., 1999; Rolink et al., 1999; Mikkola et al., 2002; Ikawa et al., 2004; Pongubala et al., 2008; Welinder et al., 2011; Mansson et al., 2012; Nechanitzky et al., 2013). However, by using mice carrying an Igll1 (λ5) reporter gene (human CD25 [hCD25]) (Mårtensson et al., 1997), it was possible to prospectively isolate B cell–committed progenitors among CD19-negative cells (Mansson et al., 2008, 2010). These B-lineage–committed population phenotypically belongs to a lineage-negative (Lin⁻) B220⁻SCA1^intKIT^intIL7R⁺FLT3⁺ common lymphoid progenitor (LP [CLP]) compartment (Kondo et al., 1997; Karsunky et al., 2008) originally thought to consist of multipotent cells with potential to differentiate to all lymphoid lineages. Further exploration of the CLP compartment revealed that functionally distinct subpopulations could be identified based on the expression of a Rag1 reporter gene (Igarashi et al., 2002; Mansson et al., 2010) or the surface marker Ly6D (Inlay et al., 2009; Mansson et al., 2010). Despite the fact that Ly6D⁺ LP cells generated mainly B-lineage cells after transplantation (Inlay et al., 2009), single-cell (SC) analysis indicated that a substantial fraction of these progenitors still possessed a T-lineage potential that could be evoked by a strong Notch signal (Mansson et al., 2010). Hence, there exists an unresolved heterogeneity in the CD19⁻ progenitor population, which obscures our current understanding of B cell commitment.

To gain insight into the earliest events in B-lymphopoiesis, we applied a strategy where we combine an antibody library screen with genome-wide expression analyses to identify heterogeneously expressed cell-surface proteins on LPs. SC gene expression analyses allowed us to link the surface markers GDNF Family Receptor Alpha 2 (GFRA2) and bone marrow (BM) stroma cell antigen 1 (BST1) to the combined expression of the B-lineage TFs Ebf1 and Pax5. While up-regulation of GFRA2 on Ly6D⁺ LPs correlated with reduction of natural killer (NK) cell but retained T and B cell potential, BST1-positive cells were largely B cell–restricted in vitro as well as in vivo. The formation of the GFRA2⁺ compartments and activation of B-lineage genes in the progenitors were directly dependent on EBF1. Our data reveal that consecutive expression of Ly6D, GFRA2, and BST1 defines a developmental trajectory linking the CLP compartment with the CD19⁺ pro–B ­­­cells in the mouse BM.

To identify surface markers that can be used to isolate subpopulations of LP cells, we explored the expression of 176 surface antigens on LP populations (CLP/LP) (Lin⁻B220⁻IL7R⁺FLT3⁺SCA1^intKIT^intLy6D^+/−; Karsunky et al., 2008; Inlay et al., 2009) from mouse BM using an arrayed antibody library and FACS analysis (Fig. 1 A and Table S1). Although several antigens, including CD102, CD44, CD47, and CD48, could be detected on all the progenitors (Table S1), CD11C, CD120B, CD5, CD162, and CD4 were expressed on a fraction of the CLP/LPs but not on CD19⁺ cells (Jensen et al., 2016). CD24 and CD157 (BST1) were expressed on subpopulations of CLPs and in a substantial fraction of the CD19⁺ cells, indicating a positive correlation with B-lineage cell fate (Fig. 1 A and Table S1). The antibody screen was complemented by RNA-sequencing (RNA-seq) of Ly6D^−/+ LPs (Inlay et al., 2009) depleted of LPAM1 (α4β7)-expressing cells (Possot et al., 2011) and a CD19⁻ B-lymphoid–restricted progenitor population expressing an Igll1 (λ5) hCD25 reporter gene (Mårtensson et al., 1997; Mansson et al., 2008; Fig. 1 B). This identified several differentially expressed surface markers that could be linked to B-lineage development, including BST1 and GFRA2.

To link the expression of differentially expressed surface markers to that of lineage-associated TFs, we used SC qRT-PCR (SC-qRT-PCR) to analyze gene expression patterns in 338 Lin⁻B220⁻IL7R⁺FLT3⁺SCA1^intKIT^int CLP/LP cells (Fig. 1 C). This identified 37 genes that were differentially expressed between subpopulations, including the surface markers Ly6d, Bst1, Gfra2, Cd24a (CD24), and Selplg (CD162), as well as the TFs Ebf1, Pax5, and Foxo1 (Fig. 1, C and D, Table S2). The expression of Bst1, Gfra2, and Cd24a displayed a positive correlation with Pax5 and Ebf1 although the opposite patterns were observed for Selplg, Itgb7, and Csf3r (CD114; Fig. 1 C). These experiments suggested that the expression of Gfra2 and Bst1 could be linked to B-lineage priming in LP compartments.

GFRA2 and BST1 are both glycosylphosphatidylinositol-anchored proteins, but whereas GFRA2 is a receptor for the neurotrophic proteins neurotropic factor (GDNF) and neurturin (Sanicola et al., 1997), BST1 is an ADP-ribosyl cyclase. Although no function of GFRA2 has been reported within the hematopoietic system, BST1 has been shown to be expressed on B cell progenitors (McNagny et al., 1988; Ishihara et al., 1996) and to be important for normal immune response (Itoh et al., 1998) as well as for neutrophil function (Funaro et al., 2004). To investigate their expression on BM progenitor populations, we performed FACS analysis for GFRA2 and BST1 on multipotent as well as B-lineage–restricted progenitors (Fig. 1 E). Although the expression of both markers was low or absent in the multipotent progenitors (LSK), they could be detected on a fraction of the classical CLP population. FrA cells were heterogeneous regarding surface levels of these proteins, whereas CD19⁺ pro–B cells were largely positive for expression of both GFRA2 and BST1. The heterogeneous expression of both GFRA2 and BST1 in progenitor compartments and rather robust expression on CD19⁺ pro–B cells indicated that these markers could be used to define early B cell progenitor populations.

Assessing the expression of GFRA2 in the CLP/LP compartment (Lin⁻B220⁻IL7R⁺FLT3⁺SCA1^intKIT^intLy6D^+/−; Karsunky et al., 2008; Inlay et al., 2009) revealed that it was largely restricted to the Ly6D⁺ population (Fig. 2 A). Expression of BST1 was even more restricted and found mainly on GFRA2⁺ cells (Fig. 2 A). To further explore whether these markers identify unique populations within the LP compartments, we examined the expression of CD19 and KIT on Lin⁻ Ly6D⁺GFRA2⁻BST1⁻, Ly6D⁺GFRA2⁺BST1⁻, and Ly6D⁺GFRA2⁺BST1⁻ cells (Fig. S1 A). All the populations contained CD19⁻ and KIT⁺ putative LP cells, in addition to CD19 expressing B cells. In line with this, CD19⁻KIT⁺ cells expressed lower levels of B220 and higher levels of FLT3 compared with Lin⁻CD19⁺KIT⁺ pro–B cells. These data support the idea that GFRA2 and BST1 can be used to identify early CD19⁻ progenitors in the mouse BM. To unravel the functional properties of the identified populations, we seeded Ly6D and GFRA2 double-negative (CLP^−/−), Ly6D single-positive (CLP^+/−), and Ly6D⁺GFRA2⁺ (CLP^+/+) progenitors on OP9 stroma cells in cultures supplemented with Kit ligand, Fms-like tyrosine kinase 3 ligand, IL-7, IL-15, and IL-2. The progressive development into B (CD19⁺) or NK (NK1.1⁺) cells was then followed by FACS analysis (Fig. 2 B). Ly6D⁺GFRA2⁺ cells underwent a rapid expansion and more efficient generation of CD19⁺ cells compared with the GFRA2⁻ populations although the generation of NK cells from GFRA2⁺ cells was minimal. Expression of CD19 in cultures of CLP^−/− cells was preceded by GFRA2 expression (Fig. 2 C), suggesting that B cell development proceeded through a GFRA2⁺CD19⁻ stage. To resolve the developmental hierarchy among the populations, we followed the formation of GFRA2- and BST1-positive cells after seeding Ly6D⁺GFRA2⁻ CLPs on OP9 stroma cells. Although substantial numbers of GFRA2⁺ cells were detected 2 d after seeding (Figs. 2 D and S1 B), BST1⁺ cells were not present until after 6 d of culture. This suggests that GFRA2 and BST1 can be used to define a developmental hierarchy in the early lymphoid compartment in the BM.

To identify residual T-lineage potential in the progenitor subpopulations, single progenitors were seeded on OP9-Delta1 (OP9DL1) cells, which promote the formation of T-lineage cells from multipotent hematopoietic progenitors (Schmitt and Zúñiga-Pflücker, 2002). The cellular content of the generated colonies was analyzed for expression of early T-lineage (CD19⁻NK1.1⁻Thy1.2^hi, or CD19⁻NK1.1⁻Thy1.2^hiCD25⁺) or B-lymphoid (CD19) markers. Although all the growing Ly6D⁺GFRA2⁻ cells displayed T-lineage potential, the Ly6D⁺GFRA2⁺ population contained cells insensitive to the NOTCH signal (Fig. 2 E). Hence, even though about half of the GFRA2⁺ cells retained some T cell potential, all B-lineage committed cells were found within the GFRA2⁺ population (Fig. 2 E). This suggests a cellular heterogeneity that potentially could be explained by changes in relative sensitivity to NOTCH signaling. To examine this possibility, we incubated Ly6D⁻ (CLP^−/−), Ly6D⁺GFRA2⁻ (CLP^+/−), and Ly6D⁺GFRA2⁺ (CLP^+/+) CLP cells on OP9DL1 cells with increasing concentrations of the γ secretase inhibitor DAPT modulating the strength of the NOTCH signal. Although we observed a profound loss of T cell formation from Ly6D⁺GFRA2⁺ progenitors already at the lowest concentrations of DAPT, formation of T-lineage cells from Ly6D⁺GFRA2⁻ cells were less sensitive to the NOTCH signal inhibitor (Figs. 2 F and S2, A and B). At higher concentrations of DAPT, we observed a significant difference in the ability of Ly6D⁺GFRA2⁻ cells to generate T-lineage progenitors compared with Ly6D⁻GFRA2⁻ cells. Thus, GFRA2 expression marks a functionally distinct population within the classical CLP compartment.

BST1 expression was highly restricted within the CLP compartment. When assessing the content of cultures generated from single BST1⁺ LP cells, we found that 90% of the clones contained B-lineage cells when cultured on OP9DL1 cells (Fig. 2 G), revealing a reduced T-lineage potential. Reanalysis of the CLP SC-qPCR data in Fig. 1 by clustering cells based on expression of Bst1, Gfra2, and Ly6d (Fig. S2 C) suggested that although the expression of Gfra2 is associated with increased expression levels of B-lineage genes including Ebf1, Foxo1, and Pax5, the levels are further increased in Bst1-expressing cells. Hence, the expression of Ly6D, GFRA2, and BST1 resolves a functional hierarchy among LP populations in the early CD19⁻ LP compartment.

To study whether the functional differences between GFRA2⁺ and negative progenitors were reflected in unique molecular signatures, we performed RNA-seq analysis from Ly6D⁻, Ly6D⁺GFRA2⁻ and Ly6D⁺GFRA2⁺ CLP/LP cells (Fig. 3 A). Among the genes up-regulated in the GFRA2⁺ cells were Ebf1 and Pax5 encoding two transcriptional regulators essential for early B cell development. To investigate if the protein levels of EBF1 and PAX5 recapitulated that of the mRNA, we used FACS to explore their presence within defined fractions of the CLP/LP compartment (Fig. 3 B and S2 D). Although low amounts of EBF1 could be detected in a minority of the Ly6D⁺GFRA2⁻ progenitors, 87% of the GFRA2⁺ cells expressed EBF1 at levels approaching that of CD19⁺ cells. PAX5 could only be detected at low levels in 22% of the GFRA2⁺ cells, whereas the protein was expressed in 99% of the CD19⁺ cells.

The increased expression of EBF1 in GFRA2⁺ cells prompted us to explore the importance of this TF for the formation of LP populations. Although the fraction of total CLPs generated after transplantation of Ebf1^−/− fetal liver (FL) LSK cells was comparable to Wt cells, the GFRA2⁺ fraction was selectively reduced (Fig. 3 C). Hence, the formation of GFRA2⁺ cells is critically dependent on EBF1. Because it has been reported that Gfra2 is a direct target for EBF1 (Treiber et al., 2010), we investigated whether the loss of GFRA2⁺ cells reflected a developmental block or primarily targeted the expression of the surface marker per se. To this end, we analyzed the association of genes differentially expressed in Ly6D⁺GFRA2⁻ compared with Ly6D⁺GFRA2⁺ CLPs with EBF1 binding sites as defined by chromatin immunoprecipitation sequencing (ChIP-seq) data from in vitro expanded pro–B cells. This suggested that almost two thirds of the genes up-regulated and about half of the genes repressed could be linked to EBF1 binding sites (Fig. 3 D). Analyzing the mRNA levels of the same set of differentially expressed genes in in vitro expanded Wt, Ebf1^−/−, or Pax5^−/− FL pro–B cells revealed that the majority of the up-regulated genes were expressed at high levels in Wt pro–B cells, whereas the expression was decreased in Ebf1^−/− progenitors (Fig. 3 E). Pax5^−/− FL cells expressed the majority of the genes at levels similar to Wt cells. Examining the chromatin accessibility in regions linked to the differentially expressed genes by assay for transposase accessible chromatin by using sequencing (ATAC-seq) analysis from Wt and Ebf1^−/− progenitor cells revealed that several of these sites were devoid of accessibility in the absence of EBF1 (Fig. 3 F). Analyzing the expression of genes expressed at lower levels in GFRA2-positive LPs revealed that a majority of these were overexpressed in the Ebf1^−/− cells (Fig. 3 G). These data strongly support the idea that GFRA2 expression defines a unique EBF1-dependent population of B-lineage primed cells within the classical Ly6D⁺ BLP compartment.

Having established GFRA2 and BST1 expression identified unique progenitor populations in the CLP compartment, we investigated the possibility of using these markers for identification of lineage-restricted progenitors in the B220⁺CD19⁻ FrA/pre-pro–B cell compartment. The protocols for isolation of these progenitors has been gradually improved (Hardy et al., 1991; Li et al., 1996; Rumfelt et al., 2006), and the highest enrichment for committed B cell progenitors has been reported for Lin⁻CD19⁻B220⁺CD24^low/−CD43⁺AA4.1^highKIT^lowIL7R⁺ cells (Rumfelt et al., 2006). Exploring surface expression of GFRA2 and BST1 expression on FrA cells (Fig. 4 A), we observed a similar pattern as for CLP/LP cells, with GFRA2⁻BST1⁻ (FrA^−/−), GFRA2⁺BST1⁻ (FrA^+/−) and GFRA2⁺BST1⁺ (FrA^−/−) populations. When seeded on OP9 stroma cells, both the GFRA2⁺ populations displayed a similar proliferative capacity; however, we detected a significant difference in the fraction of CD19⁺ cells generated from BST1⁺ progenitors after 2 d of culture (Fig. 4 B). In contrast, the GFRA2⁻BST1⁻ cells did not expand as well as the other populations, and the generation of CD19⁺ cells was slower and preceded by formation of GFRA2⁺ cells. At day 12 after seeding, we noted a reduced expression of GFRA2 on a fraction of the CD19⁺ cells generated by GFRA2⁺ progenitors, a finding likely resulting from the heterogeneous expression of this marker in CD19⁺ progenitor compartments (Fig. 1 E). Thus, as for the CLP/LP compartment, expression of GFRA2 could be linked to a high B-lineage potential. To detail this, we sorted the three populations and seeded SCs on OP9 cells. Although most (∼60%) GFRA2⁺ progenitors generated pure B cell colonies, only ∼10% of the GFRA2⁻ cells generated colonies under these conditions, and approximately one third of these was composed of NK cells (Fig. 4 C). When the progenitors were seeded on OP9DL1 stroma cells, the GFRA2⁻ cells generated almost exclusively T-lineage colonies, GFRA2⁺BST1⁻ cells generated similar numbers of B- and T-lineage clones (Fig. 4, D and E), and the BST1⁺ cells generated ∼90% B-lineage clones. Similar cloning frequencies and lineage potentials were observed after seeding single Lin⁻CD19⁻B220⁺GFRA2⁺ cells (Fig. S3, A–C), revealing reduced need for additional markers for identification of B-lineage progenitors within the FrA compartment. Transplantation of GFRA2-positive or -negative B220⁺ progenitors to immune-deficient NOD scid γc^−/− (NSG; CD45.1) hosts revealed that, even though both populations were able to generate B-lineage cells, the GFRA2 negative cells generated more CD19⁻ cells, including a substantial fraction of NK1.1⁺ cells (Fig. 4 F). Similar cellular compositions were observed in blood at 2 and 4 wk after transplantation (Fig. 4 F). In the spleen, we detected NK cells generated from the GFRA2⁻ cells, and in three of the transplanted animals, we detected CD3⁺ cells (Fig. 4 G). The BM of mice transplanted with GFRA2⁺ cells (Fig. 4 H) contained a higher fraction of donor cells, and the finding that 20% of these were CD19⁺IgM⁻ cells indicates that the GFRA2⁺ progenitors were able to remain in an immature stage for at least 4 wk. Hence, the lineage potential of GFRA2- and BST1-positive FrA cells is highly similar to that of the CLP/LP compartment.

As the functional relevance of GFRA2 in early B cell development is largely unexplored, we analyzed the presence of LP populations as well as different B-lineage compartments in the BM from GFRA2-deficient mice. We detected a comparable number of cells as well as highly similar proportions of CLPs and B cell populations in the absence of GFRA2 (Fig. S3, D–F), suggesting that this protein is not essential for B cell development in the BM.

GFRA2 signals via the tyrosine kinase receptor RET (Cik et al., 2000) has been reported to be of importance for normal hematopoiesis in the fetus (Fonseca-Pereira et al., 2014). To explore whether GFRA2 was expressed in the FL, we dissected FL from Wt mice and explored the expression of surface markers by FACS. Similar to what we observed in the BM, CD19⁺ cells in the FL expressed high levels of both GFRA2 and BST1 (Fig. S3 G). GFRA2 expression and to some extent BST1 was also detected on FL FrA (Fig. S3 I) as well as CLPs (Fig. S3 H). Despite robust expression of GFRA2 on FL cells, we noted no differences in the presence of LSK cells, CLPs, total CD19, or B1 progenitors in Wt compared with Gfra2^−/− FL (Fig. S3, J–L). Hence, although expression of GFRA2 in the LP compartments in the BM is linked to B-lineage progenitors, we were unable to detect evidence for a critical function of this protein in B cell development.

To identify potential molecular differences between GFRA2- and BST1-expressing cells in the FrA or CLP compartments, we analyzed the expression of lineage-associated genes in the defined populations by SC-qPCR (Fig. 5 A and Table S3). This revealed an association of GFRA2⁻ CLP/LP and GFRA2⁻ FrA cells with a relative absence of B-lineage priming. In contrast, B-lineage genes, including Pax5, Pou2af1, and Igll1, were expressed in both GFRA2-positive CLP/LP and FrA cells. Component analysis (C; Fig. 5 B) revealed that the GFRA2⁻ FrA cells clustered with GFRA2⁻ CLPs whereas all the GFRA2⁺ progenitor populations displayed a high degree of similarity in gene expression. These data suggest that GFRA2 expression defines similar populations in both the classical CLP compartment and FrA, arguing against these cells representing distinct progenitor populations.

To further investigate the molecular heterogeneity in B-lineage progenitor subpopulations we performed SC RNA sequencing of classic CLPs (Karsunky et al., 2008), FrA cells as well as CD19⁺IgM⁻ pro–B cells. Analyzing the data by using dimension reduction (DDRTree; Fig. 5 C) suggested a clear discrimination between the CLPs (red) and pro–B cells (blue) with FrA cells (green) and some CLPs bridging the two populations. In line with the observation that FrA is rather heterogeneous (Rumfelt et al., 2006), we detected cells from this population clustering with both the more mature B cells and the CLP compartments. The impact of cell cycle on the clustering of cells shows that the C3 dimension is largely dependent on cell cycle status and that a large portion of the examined cells from all populations was in either S or G2/M phase (Fig. 5 D). Expression levels of Ebf1 (Fig. 5 E), Pax5, or Igll1 (Fig. S4) were generally increased in the C1 and C2 dimensions of the DDRTree plot, suggesting progressive differentiation. Although we detected mRNA encoding Ly6D in several populations, the expression of Gfra2 or Bst1 was more restricted to a few cells of more differentiated phenotype. These data suggest that, even though FrA as well as the Ly6D-expressing BLP compartment is enriched for B cell–primed progenitors, both populations display a high molecular heterogeneity.

Having established that GFRA2 and BST1 expression identified committed B-lineage progenitors in both the classical CLP and FrA populations, we wanted to examine whether these markers could be used to create a unified model for early B cell development in the BM. To identify committed and primed progenitors in the lineage negative compartment, we used reporter transgenic mice carrying Rag1-GFP and Igll1-hCD25 reporter constructs (Mansson et al., 2010; Fig. S5 A). Even though the Lin⁻IL7R⁺KIT^Int CLP population was highly enriched for reporter-expressing cells, a substantial fraction of the Rag-expressing cells could be detected in the IL7R⁺KIT⁻ fraction, indicating that B-lineage progenitors might be identified outside of the classical CLP (Kondo et al., 1997) or FrA (Rumfelt et al., 2006) compartment.

As would be predicted from the analysis of reporter mice, we detected expression of Ly6D, GFRA2, and BST1 in both the KIT^Int and the KIT⁻ populations (Figs. 6 A and S5 B). qPCR analysis suggested that the expression of GFRA2 was linked to increased expression of Pax5 mRNA (Fig. S5 C), and exploring Ig heavy chain (IgH) DJ recombination in SCs (Rumfelt et al., 2006) identified fewer cells with remaining germline alleles in the GFRA2⁺ compartment of the Ly6D⁺ cells (Fig. S5 D).

Investigation of the myeloid potential of SCs in vitro suggested that high expression of IL7R was sufficient to identify progenitors with reduced myeloid potential (Fig. 6 B). Hence, even though a few myeloid clones with more than 100 cells were generated from IL7R⁺Ly6D⁻ cells, the inclusion of additional markers would be of limited value. Exploring B-lineage potential of SCs suggested that expression of IL7R and intermediate expression of KIT was associated with higher cloning frequencies (Fig. 6 C). However, we detected cells with potential to generate B-lineage clones in both IL7R⁺KIT^int and IL7R⁺KIT⁻GFRA2⁺ cells. Hence, GFRA2 expression identifies cells with B-lineage potential outside the IL7R⁺KIT^Int compartment harboring both the CLP (Kondo et al., 1997) and FrA (Rumfelt et al., 2006) populations. Cultivation of the defined progenitors under conditions stimulating T cell development (Fig. 6 D) generated data highly similar to what we observed for purified CLP and FrA cells (Figs. 2 and 4). GFRA2 expression was associated with formation of both B and T cell clones, whereas the expression of BST1 was linked to a clear dominance of B cell clones (Fig. 6 D). Although the expression of KIT was associated with increased cloning frequency, this had no impact on lineage potential. Exploring the ability of the populations to generate NK cells, we noted NK cell production by IL7R⁺Ly6D⁻ as well as the IL7R⁻KIT⁻B220⁺ progenitors (Fig. 6 E). Although Ly6D expression generally associated with a reduction in the ability to generate NK cells (Inlay et al., 2009), we detected a virtual loss of NK cell potential among all the GFRA2-expressing populations. Within the IL7R⁺ compartment, the expression of KIT could be linked to higher cloning frequency, although this had most limited impact on the lineage potential in the cells. Hence, SC in vitro analysis suggests that the Ly6D⁺ population can be divided into GFRA2-negative cells (BLP1), GFRA2⁺BST1⁻ cells, (BLP2), and GFRA2⁺BST1⁺ progenitors (BLP3).

To resolve the in vivo potential of the progenitors, we transplanted 2,000 Ly6D⁻ KIT^int/-IL7R⁺, BLP1, BLP2, and BLP3 cells as well as IL7R⁻KIT⁻B220⁺ progenitors to CD45.1⁺ NSG mice. We examined the presence of CD45.2⁺ cells in the spleen and BM at 4 wk after transplantation (Fig. 6, F and G). Although the Ly6D⁻KIT^int/-IL7R⁺ as well as the BLP1-2 cells all generated comparable levels of CD45.2⁺ cells in the spleen (Fig. 6 G), BLP3 cells were less efficient. B220⁺IL7R⁻ cells were largely unable to generate CD45.2⁺ cells in the spleen. All the reconstituting populations generated mainly CD19⁺ B-lineage cells, and none of the transplanted LP populations generated a macroscopically identifiable thymus at 4 wk after transplantation. However, we detected some formation of CD3⁺ cells in the spleen and BM from the Ly6D⁻ KIT^int/-IL7R⁺ progenitors and BLP1 cells. The generation of NK and myeloid cells in the spleen was low from all the tested populations. The BM generally contained lower levels of CD45.2⁺ cells (Fig. 6 F). Although all the populations generated CD19⁺ cells, the Ly6D⁻ KIT^int/-IL7R⁺ and BLP1 as well as the IL7R⁻B220⁺ cells generated a substantial fraction of CD3⁺ cells. These cells displayed low expression of Pax5 and high mRNA levels encoding the TFs of Gata3 and Tcf7 (Fig. S5 E), verifying that they represent T-lineage cells. Hence, expression of GFRA2 on LP cells is associated with functional commitment to B-lymphoid development.

Although seminal work in the 1990s (Sen et al., 1990; Hardy et al., 1991; Li et al., 1993, 1996; Rolink et al., 1994; Kondo et al., 1997; Tudor et al., 2000) created a foundation for our understanding of B-lineage development, the identity of the earliest committed B cell progenitor has been rather elusive. Despite the fact that expression of B220 enriches for B cell progenitors and appears to proceed the expression of GFRA2 in in vitro differentiation assays (Fig. S1 B), it has been reported to have a limited value as a marker for early committed cells (Martin et al., 2003; Rumfelt et al., 2006; Mansson et al., 2010; Alberti-Servera et al., 2017). Furthermore, even though surface expression of Ly6D identifies cells largely B-lineage restricted in vivo (Inlay et al., 2009), a small population of T-lineage cells was detected in the thymus after transplantation (Inlay et al., 2009). The high-resolution analysis presented in this study unravels both functional and molecular heterogeneity within the Ly6D⁺ BLP population.

The data presented in this study show that the expression of B-lineage–restricted genes is largely confined to the GFRA2⁺ BLP2-3 compartments and that GFRA2-negative progenitors retain the ability to adopt T-lineage fate when exposed to a strong NOTCH signal or when transplanted to lymphopenic mice. Even though we were able to detect CD3⁺ cells generated from Ly6D⁺ progenitors, we failed to identify a macroscopically detectable thymus 4 wk after transplantation. This is likely a consequence of the limited ability of BM-derived LPs to sustain long-term generation of T-lineage cells (Allman et al., 2003) in combination with the use of adult NSG mice as hosts for transplantation experiments.

GFRA2 expression is tightly linked to the expression of EBF1, and a substantial part of the genes associated selectively with the BLP2-3 stage appears to be direct targets for EBF1. This is well in line with the idea that EBF1 is essential for the activation of B-lineage restricted transcriptional programs (Zandi et al., 2008; Treiber et al., 2010) in early progenitor cells. EBF1 expression has also been linked to lineage restriction because ectopic expression of EBF1 has been reported to result in reduced development of T- and NK-lineage cells in vivo (Zhang et al., 2003), and inactivation of the Ebf1 gene in B-lineage cells results in lineage plasticity (Mikkola et al., 2002; Nechanitzky et al., 2013; Ungerbäck et al., 2015). This is likely explained by the ability of EBF1 to suppress expression of Gata3 (Banerjee et al., 2013) and Id2 (Thal et al., 2009), which are essential for development of T and NK cells, respectively. Interestingly, the expression of Bst1 has been associated with functional PAX5 (Pridans et al., 2008; Revilla-I-Domingo et al., 2012) in line with that, and expression of this marker is associated with commitment to the B cell lineage. Although the generation of GFRA2⁺ cells is dependent on EBF1, the expression of Ly6D would appear to be under the control of alternative regulatory circuits (Fig. 3 C; Tsapogas et al., 2011). One possibility could be that this reflects the action of E-proteins because of the importance of these TFs for the formation of Ly6D⁺ progenitors (Inlay et al., 2009; Welinder et al., 2011). The same regulatory network may be involved in the regulation of Ig recombination because, even though we detected an increased frequency of recombination events in the GFRA2 expressing cells, Rag1 expression as well as recombination events are detected outside this compartment.

Blood cell generation is often depicted in a linear model in which multipotent progenitors gradually lose functional potentials as they develop along branches representing the different blood lineages. However, emerging evidence suggests that the developmental trajectories may not follow a strict tree-like structure but rather resemble a network of paths where the goal can be reached via different routes (Paul et al., 2015; Perié et al., 2015). This concept is rather well established in fetal development where progenitor populations with combined B-macrophage or T-macrophage potential (Katsura, 2002) as well as B/T and B/NK bipotent progenitors are found in the FL (Possot et al., 2011; Pereira de Sousa et al., 2012). Our analysis as well as previous work (Martin et al., 2003; Rumfelt et al., 2006; Inlay et al., 2009; Mansson et al., 2010) suggest that B-lineage–committed cells can be found both in the CLP/LP compartment and in FrA. These data could be interpreted as if B cell development proceeds along two parallel developmental trajectories. However, both functional and molecular analysis suggested that GFRA2-expressing cells from the B220⁻ (CLP) and B220⁺ (FrA) populations represented a highly similar or even identical developmental stage. Even though it should not be excluded that development may proceed directly from a GFRA2⁺ stage into the CD19⁺ compartment, our data support that the BST1⁺ population indeed represents a more mature compartment of CD19⁻ cells. Our inability to identify committed B cells outside of the GFRA2⁺ progenitor compartments strongly suggests that the developmental trajectory defined by this marker represents the major pathway for B-lineage development from a functional CLP into a CD19⁺ pro–B cell. Even though the expression of KIT did not impact the lineage potential of the progenitor cells, it should be noted that the frequency of functional progenitor cells was reduced in the KIT-negative populations. Based on our findings, we would propose that the pro–B cell (BLP; Inlay et al., 2009) compartment can be divided into Ly6D⁺GFRA2⁻ (BLP1), GFRA2⁺BST1⁻ (BLP2), and GFRA2⁺BST1⁺ (BLP3) cells (Fig. 7). We believe this subdivision of early progenitor cells will increase our understanding of molecular events regulating normal as well as malignant B cell development.

Wt, Gfra2^−/− (Rossi et al., 1999), Ebf1^−/− (Lin and Grosschedl, 1995), an hCD25 reporter under the regulatory elements of the Lambda5 (λ5, Igll1) promoter (Mårtensson et al., 1997), and Rag1-GFP- (Igarashi et al., 2002) mice were on a C57BL/6 background and backcrossed for more than seven generations. BM samples were harvested from Wt at 6–11 wk and from Gfra2^−/− mice at 8 wk or embryonic day 18 FL. Comparative analyses were based on littermate controls. Animal procedures were performed with consent from the local ethics committees at Lund University and Linköping University.

For isolation of BM LP cells, bones were crushed in a mortar, and the cell suspension passed through a 50-µm filter. For progenitor cell isolation by FACS sorting, BM cells were either sorted unenriched, enriched for KIT⁺ cells by magnet-activated cell sorting column for using anti-CD117 immunomagnetic beads (Miltenyi Biotec), or lineage depleted by using sheep anti–rat IgG beads (Dynabeads; Invitrogen), recognizing purified lineage antibodies against CD4 (GK1.5; BioLegend), CD8A (53–6.7; BioLegend), CD11B (M1/70; BioLegend), Ter119 (TER-119; BD Biosciences), GR1 (RB6-8C5; BioLegend), and CD19 (6D5; BioLegend) before antibody staining. Cells were thereafter stained with antibodies (full information of antibody combinations, clones, and conjugates is shown in Table S4). Cell sorting was performed on BD FACSAriaII or FACSAriaIII cell sorters (BD Biosciences) and analyses on a BD FACSLSRII, FACSAriaIII, FACSCantoII, or FACSFortessa.

The Lyoplate analysis was performed on 5 million unfractionated BM cells per sample according to the manufacturer’s instructions (BD Biosciences). In addition to the Lyoplate antibodies, cells were costained with surface markers identifying CD19⁺ and CLP/LP cells (see Table S4 for full antibody information).

For evaluation of B, NK, and T cell potential, SCs were deposited directly into 96-well plates containing preplated (2,000 cells/well) stroma cells. T-cultures (on OP9-DL1 stroma layers) were supplemented with 10 ng/ml KIT ligand, 10 ng/ml Fms-like tyrosine kinase 3 ligand (FLT3L), and 10 ng/ml IL-7. B/NK cultures (on OP9 stroma layers) were supplemented with the same cytokines as the OP9-DL1 cultures, with the addition of 20 ng/ml IL-15 and 40 ng/ml IL-2. All cytokines were acquired from PeproTech. Cultures were substituted with fresh cytokines after 7 d. OptiMEM supplemented with 10% heat-inactivated FCS, 50 µg/ml gentamicin, and 50 μM β-mercaptoethanol was used for maintaining the OP9/OP9-DL1 stromal cell lines as well as for the co-cultures. Co-cultures were evaluated by FACS with CD19 (ID3), NK1.1 (PK136), CD11c (N418), CD45 (30-F11), CD11b (M1/70), GR1 (RB6-8C5), CD25 (7D4), and CD90.2/Thy1.2 (53–2.1) antibodies. Cells were considered B cells if CD45⁺CD19⁺, NK cells if CD45⁺NK1.1⁺, and T cells if CD45⁺CD90.2⁺CD19⁻NK1.1⁻CD11c⁻CD11b⁻GR1⁻ or CD45⁺CD90.2⁺CD25⁺CD19⁻NK1.1⁻CD11c⁻CD11b⁻GR1⁻. Co-cultures were analyzed at days 12–16 with a BD FACSLSRII analyzer (BD Biosciences).

Myeloid potential of SCs in vitro was determined by cultivation of cells in OptiMEM supplemented with 10% heat-inactivated FCS, 50 µg/ml gentamicin, and 50 μM β-mercaptoethanol. The medium was complemented with 5 ng/ml IL-3, 5 ng/ml GM-CSF, 5 ng/ml M-CSF, 10 ng/ml G-CSF, 25 ng/ml Kit-ligand, and 25 ng/ml FLT3L (PeproTech). SC dispersion was based on limiting dilution (170 sorted cells diluted into 3.4 ml of medium) to generate on average one cell per 20 µl per well.

SCs, with 10- and 20-cell controls, were sorted into 96-well plates with 4 µl lysis buffer (10% NP-40, 10 mM deoxynucleotide triphosphate, 0.1 M 1,4-dithiothreitol, RNAsout [Invitrogen], and nuclease-free water) and immediately frozen at −80°C. Upon thawing, Taq/SSIII enzyme reaction mix (CellsDirect One-Step RT-qPCR kit; Invitrogen) and 96 TaqMan assays at a final dilution of 0.025× were added to each well, excluding four wells where no RT control of Taq enzyme (Invitrogen) reaction mix was added. Preamplification cycling was preformed according to previously described protocols (Pina et al., 2012) with minor alterations; the first cycle, 50°C at 15 min, was extended to 60 min. The preamplified SCs were loaded on to a 48.48 or 96.96 Dynamic Array Chip for Gene Expression (Fluidigm) with 48 or 96 TaqMan assays, respectively. The unfractionated CLP Fluidigm data were analyzed by using the SCExV webtool (Lang et al., 2015). Gene expression was normalized by using the median expression and positive and negative control wells as well as cells showing the same expression patterns as those that were removed. Random forest clustering was applied as described in (Shi, 2006), and the resulting grouping was reordered to follow the expression changes of PAX5, EBF1, and CD114. Statistics on differentially expressed genes were calculated using the MAST R package according to the hurdle model, and differentials were called at a cut-off of P < 0.01 (Finak et al., 2015). The subfractionated CLP and FrA Fluidigm data were read into a Rscexv object, and all cells showing XRNA expression >25 CT (cycle of threshold) as well as the multicell controls were dropped. The data were normalized to the median expression per cell and median absolute deviation z-scored. Unsupervised clustering was performed by using the random forest clustering implemented in Rscexv. Subsequently, the expression of the B cell marker genes (Pax5, Rag1, Ebf1, and Igll1) was used to manually group the FACS sort and random forest clustered samples into the B cell–specific groups “not B cell committed,” “preparing,” and “B cell committed.” Statistics on differentially expressed genes were calculated using the MAST R package according to the hurdle model, and differentials were called at a cut-off of P < 0.05 (Finak et al., 2015).

Single Lin⁻CD19⁻Kit^int/negIL7R⁺Ly6D⁺ cells were index sorted, based on the mentioned markers and additionally GFRA2, BST1, and B220, directly into lysis buffer in 96-well plates. D_HJ_H rearrangement status was investigated as previously described (Rumfelt et al., 2006), with the exceptions that the Advantage 2 Polymerase Mix and PCR buffer (Takara Bio) were used, and PCR products were visualized on 2% agarose gels (E-Gel96; Invitrogen) by using an E-Base integrated power supply device (Thermo Fisher Scientific). The validity of these PCR products was not verified by sequence analysis. SC-qPCR data can be found in GEO accession no. GSE92540.

Total RNA for RNA-seq was isolated from sorted CLP/LPs by using RNAeasy Micro Kit (Qiagen) according to manufacturer’s recommendations. Libraries were constructed by using NuGEN’s Ovation Ultralow Library systems (NuGEN Technologies) and were subsequently subjected to 50–76 cycles of sequencing on a HiSeq 2000 or NextSeq500 (Illumina). For analysis of LPAM1⁻Ly6D⁻ CLPs, LPAM1⁻Ly6D⁺ CLPs, and hCD25 (Igll1)⁺ LPs, reads were mapped by using HISAT (Kim et al., 2015) against the mm10 version of the mouse genome and counted against the Gencode M4 annotation. Differentially expressed genes were identified by using DESeq2 (Love et al., 2014) for R using a false discovery rate (FDR) <10% and reduced to surface marker proteins only by using a surface protein atlas (Bausch-Fluck et al., 2015). Replicate samples were averaged before generating heatmaps. For analysis of Ly6D⁻ GFRA2⁻, Ly6D⁺ GFRA2⁻, and Ly6D⁺ GFRA2⁺ CLPs and hCD25⁺ LPs, the reads were aligned to mouse reference genome (mm10/GRCm38) using TopHat (Kim et al., 2013). If not indicated, further analyses were performed using the HOMER platform (Heinz et al., 2010). After conversion of Bam-files to Tag directories with the command makeTagdirectory.pl, samples were analyzed by using normalization to 10M mapped reads by the analyzeRepeats.pl command with the options –count exons -condeseGenes. For analysis of statistical significance among differently expressed genes between the three different CLP populations, the data were analyzed using analyzeRepeats.pl with the –noadj option followed by the getDiffExpression.pl command by using EdgeR (Robinson et al., 2010) analyzing -AvsA with parameters -fdr 0.10 -log2fold 1. Heatmaps were created by hierarchical clustering in Cluster3 (de Hoon et al., 2004) followed by visualization in Java Treeview (Saldanha, 2004).

RNA-seq data from Wt, EBF^−/−, and Pax5^−/− pro–B cells, ATAC-seq data from Wt and EBF^−/− pro–B cells, and pro–B EBF1 ChIP-seq data were mapped against reference genome (mm10/GRCm38) by using TopHat (Kim et al., 2013) and Bowtie 2 (Langmead and Salzberg, 2012), respectively. Further analyses were performed using HOMER (Heinz et al., 2010) with an initial conversion of Bam files to tag directories by using makeTagdirectory.pl. Expression matrices of RNA-seq data normalized to 10M mapped reads were created by the analyzeRepeats.pl command with the options –count exons -condeseGenes. ATAC-seq peaks were identified by using findPeaks.pl in HOMER with the settings -style histone -size 75 -minDist 75 -minTagThreshold 6 -L 8 -F 8. ChIP-seq peaks were identified by using findPeaks.pl in HOMER with the setting -style factor. Files for fishbone heatmaps were created by using annotatePeaks.pl -size 5000 -hist 25 -ghist, and files for histograms were created using annotatPeaks.pl with the options -size 5000 -hist 25 in HOMER.

Heatmaps were finalized by hierarchical clustering in Cluster3 (de Hoon et al., 2004) followed by visualization in Java Treeview (Saldanha, 2004). Gene expression heatmaps were clustered by using log transformation center genes on means, and average linkage (uncentered correlation), hierarchical clustering of genes and arrays, and fishbone heatmaps by clustering of genes only by using average linkage (uncentered correlation). Publicly available data retrieved from GEO and used in this report was pro–B EBF1 ChIP (accession no. GSE69227; Ungerbäck et al., 2015). RNA-seq and ATAC-seq data from Wt and EBF1-deficient cells as well as CLP subfractions were deposited in GEO (accession no. GSE92434 and GSE92540).

SC RNA-seq data were generated on the 10X emulsion platform (10X Genomics) according to the manufacturer’s instructions. The SC RNA-seq from the Chromium platform was processed using CellRanger v1.3 and default settings, and subsequent normalization, quality control, filtering, differential gene expression analysis, and cell-cycle assignment were performed in R by using Seurat v2.2 (Butler et al., 2018). Dimension reduction was performed using DDRTree over three components. SC RNA-seq data was deposited in GEO (accession no. GSE92540).

Subfractions of CLP and FrA cells were plated at a concentration of 50 cells per well in 12-well plates containing preplated (27,000 cells per well) OP9 stroma cells. Cultures were kept in OptiMEM supplemented with 10% heat-inactivated FCS, 50 µg/ml gentamicin, and 50 μM β-mercaptoethanol and were supplemented with 10 ng/ml KIT ligand, 10 ng/ml FLT3L, 10 ng/ml IL-7, 20 ng/ml IL-15, and 40 ng/ml IL-2. All cytokines were acquired from PeproTech. Cultures were substituted with fresh cytokines after 7 d. Cell cultures were analyzed by FACS at days 2, 4, 6, 8, and 12 and stained against NK1.1 (PK136; BioLegend), CD19 (1D3; eBioscience), CD45.2 (104; BioLegend), BST1 (BP3; BioLegend), B220 (RA3-6B2; BioLegend), and GFRA2 (R&D Systems). Cell expansion in cultures was estimated by FACS through the addition of a standardized number of beads (Trucount Control; BD Biosciences) to all FACS samples.

Bones (femur, tibiae, crista iliac) from 10–12-wk-old CD45.2 mice were taken, and hematopoietic BM cells were extracted by crushing the bones in mortar. Red cells were removed by adding ice-cold Erythrocyte Lysis Buffer (1.5 M NH₄Cl, 100 mM NaHCO₃, and 10 mM EDTA, pH 7.4). Lineage-positive cells were removed by staining with purified antibodies CD19(6D5), CD8a(53–6.7), CD4(GK1.5), Gr1(RB6-8C5), Mac1(M1/70), and Ter119(TER-119), followed by negative selection by using Dynabeads Sheep Anti–Rat IgG (11035; Invitrogen). Remaining cells were next incubated with biotinylated mGFRα-2 (BAF429; R&D Systems) and further stained with lineage-PECy5 Gr1(RB6-8C5), Mac1(M1/70), Ter119(TER-119), CD11c(N418), NK1.1(PK136), CD3 ε(145−2c11), CD19-PECy7 (ID3; eBioscience), B220-Alexa Flour 700 (RA3-6B2), and Streptavidin-Brilliant Violet (BV) 605. Lin⁻CD19⁻B220⁺GFRA2⁺ and Lin⁻CD19⁻B220⁺GFRA2⁻ cells were sorted on BD FACSAriaIII, and 2,000 cells per 150 µl 1% PBS/FCS were transplanted by intravenous injection via tail vein to 9–11-wk-old preconditioned (2.25 Gray) NOD.Cg-Prkdc^scid Il2rg^tm1WjI/SzJ (NSG) mice. Peripheral blood was analyzed 2 wk after injection, and at 4 wk peripheral blood, spleen, and BM were taken for flow cytometric analysis by using IgM-APCCy7 (RMM-1) and CD4-APC (GK1.5) or CD49b-APC (DX5), CD19-PECy7 (ID3; eBiosciences), CD3ε-BV605 (17A2), CD8-PacificBlue (53–6.7), NK1.1-FITC (PK136), CD45.1-PE (A20; eBiosciences), and CD45.2-Alexa Fluor 700 (104). All antibodies are from BioLegend, unless otherwise indicated.

For analysis of EBF1-deficient cells, 2,000 sorted LSK cells from Wt or Ebf1^−/− FL were transplanted into sublethally irradiated NSG mice by tail vein injection. The BM samples of transplanted mice were subjected to FACS analysis 6 wk after transplantation to identify LP populations.

Subfractions of CLP cells were plated at a concentration of 20 cells per well in 24-well plates containing preplated (13,000 cells per well) OP9DL1 stroma cells. Cultures were kept in OptiMEM supplemented with 10% heat-inactivated FCS, 50 µg/ml gentamicin, and 50 μM β-mercaptoethanol and were supplemented with 10 ng/ml KIT ligand, 10 ng/ml FLT3L, and 10 ng/ml IL-7 (PeproTech). Cultures were substituted with fresh DAPT and cytokines after 7 d. Cultures were analyzed by FACS at day 14 and stained against CD19 (ID3), NK1.1 (PK136), CD11c (N418), CD45 (30-F11), CD11b (M1/70), GR1 (RB6-8C5), and CD90.2/Thy1.2 (53–2.1).

WT mice were individually lineage-depleted by using sheep anti–rat IgG beads (Dynabeads; Invitrogen) recognizing purified lineage antibodies against CD4 (GK1.5; BioLegend), CD8A (53–6.7; BioLegend), CD11B (M1/70; BioLegend), Ter119 (TER-119; BD Biosciences), and GR1 (RB6-8C5; BioLegend) and subsequently stained against surface markers (Table S4). Fixation, permeabilization, and staining of TFs were done by using the Transcription Factor Buffer Set (BD Pharmingen) according to the manufacturer’s instructions.

Table S1 includes detailed data from the antibody screening process. Tables S2 and S3 show differentially expressed genes in the SC-qPCR experiments. Table S4 shows the FACS antibodies used to identify LP populations. Tables S5 and S6 show statistical analysis of SC in vitro culture data. Fig. S1 shows back gating of defined populations and FACS data from kinetics experiments. Fig. S2 shows SC PCR analysis of CLP cells and data from DMSO control experiments for the Notch inhibition experiments in Fig. 2 as well as intracellular FACS analysis of EBF1 and PAX5. Fig. S3 shows in vitro analysis of B220⁺ progenitor populations and FACS data from defined progenitors in BM and FL in Wt and Gfra2^−/− mice. Fig. S4 shows DDRTree projections highlighting gene expression patterns obtained from SC-RNA-seq analysis, and Fig. S5 shows FACS plots from reporter transgenic mice, gating strategies, qPCR data from sorted populations, and SC IgH DJ analysis.

We thank Liselotte Lenner, Linda Bergström, Gerd Sten, Maria Malmberg, Eva Erlandsson, and Maria Björklund for expert technical assistance.

This work was supported by grants from the Swedish Cancer Society, the Swedish Childhood Cancer Foundation, the Swedish Research Council, including Center grants to Stemtherapy and Biocare, and the Knut and Alice Wallenbergs Foundation, and by donations from Henry Hallberg and Lund as well as Linköping Universities.

The authors declare no competing financial interests.

Author contributions: C.T. Jensen, J. Åhsberg, M.N.E. Sommarin, K. Okuyama, R. Somasundaram, T. Strid, and J. Ungerback designed, conducted, and analyzed the experiments. S. Soneji, S. Lang, and T. Strid performed bio-informatics analysis of data. M.S. Airaksinen and J. Kupari contributed to the experiments with GFRA2^−/− mice. D. Bryder aided in the SC-RNA-seq experiments. S. Soneji, D. Bryder, M. Sigvardsson, and G. Karlsson designed experiments and analyzed data. All authors contributed to the writing of the manuscript.