B1 cells are self-renewing innate-like B lymphocytes that provide the first line of defense against pathogens. B1 cells primarily reside in the peritoneal cavity and are known to originate from various fetal tissues, yet their developmental pathways and the mechanisms underlying maintenance of B1 cells throughout adulthood remain unclear.
We performed high-throughput single-cell analysis of the transcriptomes and B-cell receptor repertoires of peritoneal B cells of neonates, young adults, and elderly mice. Gene expression analysis of 31,718 peritoneal B cells showed that the neonate peritoneal cavity contained many B1 progenitors, and neonate B cell specific clustering revealed two trajectories of peritoneal B1 cell development, including pre-BCR dependent and pre-BCR independent pathways. We also detected profound age-related changes in B1 cell transcriptomes: clear difference in senescence genetic program was evident in differentially aged B1 cells, and we found an example that a B1 subset only present in the oldest mice was marked by expression of the fatty-acid receptor CD36. We also performed antibody gene sequencing of 15,967 peritoneal B cells from the three age groups and discovered that B1 cell aging was associated with clonal expansion and two B1 cell clones expanded in the aged mice had the same CDR-H3 sequence (AGDYDGYWYFDV) as a pathogenically linked cell type from a recent study of an atherosclerosis mouse model.
Beyond offering an unprecedent data resource to explore the cell-to-cell variation in B cells, our study has revealed that B1 precursor subsets are present in the neonate peritoneal cavity and dissected the developmental pathway of the precursor cells. Besides, this study has found the expression of CD36 on the B1 cells in the aged mice. And the single-cell B-cell receptor sequencing reveals B1 cell aging is associated with clonal expansion.
B lymphocytes contribute to both humoral and cellular immunity, serving as antibody-producing and antigen-presenting cells . B cells are generally divided into two major subsets: B1 and B2 cells. B2 cells generate the well-described germinal center and T-cell-dependent immune responses against protein antigens, whereas B1 cells produce a rapid T-cell-independent antibody response—primarily targeting lipids and polysaccharides [2, 3]. In mouse, B1 cells primarily reside in body cavities including the peritoneal cavity—a fluid-filled space located between the wall of the abdomen and organs [4, 5]. B1 cells secrete most of the circulating natural antibodies, including those which recognize endogenous (i.e., self) antigens detectable on apoptotic cells and oxidized lipids [6,7,8], natural anti-influenza virus IgM antibodies , and gut mucosa protective IgA antibodies . B1 cells also secrete immune regulatory cytokines including interleukin-10  and granulocyte-monocyte colony stimulating factor . Both regulatory cytokines and natural antibodies produced by B1 cells are essential for maintaining physiological and immune homeostasis [3, 4]. Furthermore, peritoneal B1 cells have been shown to exert phagocytic and microbicidal abilities .
Recent studies have revealed that B1 cell precursors are distinct from B2 cell precursors in terms of their distinct transcriptional signature and divergence in their networks of transcription factors governing specification and fate commitment [14, 15]. B1 cell development occurs in two fetal waves and one adult wave [15, 16], whereas B2 cells are generated and continuously replenished via hematopoiesis in adults. Precursors of B1 cells have been detected in diverse tissues, including the yolk sac , para-aortic splanchnopleure , aorta-gonad-mesonephros , placenta , fetal and neonatal liver , fetal omentum [22,23,24], spleen , and bone marrow . While B1 cell progenitors are present in adult bone marrow, several lines of evidence indicate that their contribution to adult B1 cell pool is quite limited under physiological condition [26,27,28,29]. A recent lineage tracing experiment with the hematopoietic stem cell restricted Pdzk1ip1-CreER line in tamoxifen-treated adult mice revealed that the contribution of adult hematopoiesis to B1 cell subset is below 10% over a 40-week observation period .
Given the various origins of B1 cell precursors and multiple development waves, it follows that the postnatal B1 cell pool must be a heterogeneous population, and likely one that differs substantially with aging. We do know that B1 cells are maintained by self-renewal in adults [27, 30,31,32] and their functions are known to be affected by aging [33,34,35,36,37,38,39]. Of particular note, a report revealed that passive transfer of IgG-depleted serum from elderly (18–24 months) mice had no effect on immune-deficient CB17-SCID mice which infected with streptococcus pneumoniae, whereas IgG-depleted serum from young adult (3 months) mice was protective . This result clearly demonstrates that the protection against pneumococcal infection afforded by natural IgM antibodies in vivo diminishes with age. Another well-known age-related phenomenon with B1 cells is their general expansion in senescent mice . However, much remains unclear about the mechanisms underlying many of these differences observed in aged mice.
We also know that the B-cell receptor (BCR) repertoire of B1 cells is affected by aging: dynamic changes in the B1 cell BCR repertoire during aging were recently profiled by sequencing transcripts of the complementarity determining region 3 of the heavy chain (CDR-H3) . This elegant study revealed that the BCR repertoire of B1 cells is continuously shaped by stringent selection throughout the lifetime of a mouse, and becomes progressively restricted with age. However, this study was conducted using bulk sequencing approaches, which required pre-sorting of B cell subsets and which limited their capacity to characterize cell-to-cell variations in the BCR repertoire of aging B1 cells.
Here, using single-cell RNA-sequencing (scRNA-seq) of peritoneal CD19+ cells from neonates, young adults, and elderly mice, we revealed that a large majority of peritoneal B cell subsets were B1 cells in adults. We found that neonate peritoneal B1 cell progenitors develop along two distinct trajectories towards B1 cells, including a pre-BCR dependent pathway and a pre-BCR independent pathway. Aging-related changes were obvious in the transcriptomes of B1 cells from the variously aged mice, and we discovered that a subset only present in the elderly mice could be easily sorted based on the surface marker CD36. Subsequently, single-cell B-cell receptor sequencing (scBCR-seq) of peritoneal CD19+ cells from the variously aged mice experimentally tracked how the BCR repertoire of peritoneal B1 cells became progressively restrictive with age. Moreover, our scBCR-seq data also support profiling of B1 cell clonal expansion, and we detected the accumulation of cells with an atherosclerosis-associated CDR-H3 sequence.
A large majority of B cells in the mouse peritoneal cavity are B1 cells
In order to investigate the B1 cell development in the peritoneal cavity, we purified individual CD19+ cells in the peritoneal cavity by magnetic-activated cell sorting from neonates, young adults, and elderly mice for scRNA-seq and scBCR-seq analysis (Fig. 1a). After sequencing quality control (Additional file 1: Fig. S1a−d) and depletion of the contaminated non-B cells (Additional file 1: Fig. S1e−g), we analyzed 31,718 B cells derived from mice of the three age groups (Fig. 1b). Dimensional reduction analysis by uniform manifold approximation and projection (UMAP) revealed that almost all cells expressed the known B cell regulators including Cd19, Ighm, Tcf3, Ebf1, Pax5, and Foxo1 (Fig. 1c), confirming that the examined cells were indeed B cells. Fcer2a (which encodes CD23) and Ighd were used as markers to distinguish B2 from B1 cells . Both Fcer2a and Ighd were highly expressed in the left bottom subset of the UMAP projection (Fig. 1d), while the B1 cell genes Bhlhe41 [41, 42] and Zbtb32 [42, 43] were highly expressed in most of the other cells (Fig. 1d).
The mutually exclusive expression of the mature B cell gene Ly86  and B cell precursor gene Vpreb3  suggested the presence of both mature B cells and their precursors in examined CD19+ cells (Fig. 1e). As mentioned above, the mature peritoneal B cells contain B1 cells (Bhlhe41+Fcer2a−) and B2 cells (Bhlhe41−Fcer2a+). Thus, the peritoneal CD19+ cells were roughly grouped into 3 clusters: B1 cells, B2 cells, and B cell precursors (Fig. 1f; Additional file 2: Table S1). Interestingly, all B cell precursors were sourced from neonates and accounted for 63% of the total peritoneal B cells of neonatal mice (Fig. 1g). The mature B1 cells accounted for 21% in neonatal mice, while their ratio increased to 87% in young adults and 88% in elderly mice (Fig. 1g).
In addition to the reported B cell regulators Sox4  and Cd93 , we discovered potential new genes that were expressed in peritoneal B cell precursors, for example, Tifa, Atp1b1, Chchd10, and Arl5c (Fig. 1h). Our scRNA-seq analysis also revealed additional B1 cell genes including Ctla4, Anxa2, Zbtb20, Sox5, Plac8, and Ass1, as well as B2 cell genes such as Icosl, Ms4a4c, Neurl3, Cr2, and Ccr7 (Fig. 1h). These scRNA-seq results for neonates, young adults, and elderly mice show that a large majority of B cells in the mouse peritoneal cavity are B1 cells and that B cell precursors are present in the neonate peritoneal cavity.
The neonate peritoneal cavity contains a high proportion of B1 cell progenitors
To analyze peritoneal B cell precursors in early life, we re-clustered the neonatal samples (Fig. 2a). A distinct aggregation pattern of B cells was observed in the re-clustered UMAP projection, wherein all of the neonatal cells were grouped into nine clusters (Fig. 2a and Additional file 1: Fig. S2a and Additional file 2: Table S1). The cells in clusters 1, 2, 3, and 4 (denoted as C1-C4) expressed no detectable Ms4a1 (which encodes CD20) but expressed Il7r and Cd93 (Fig. 2b). C1 and C4 cells expressed the known DNA recombination genes Rag1 and Rag2 (Fig. 2b). C1-C2 cells expressed a pre-BCR surrogate light chain gene Vpreb1 (Fig. 2b).
The expression of the early B cell genes suggested that C1-C4 cells were early B cell progenitors. To experimentally validate neonate peritoneal B cell progenitors, we performed flow cytometry analyses. Consistent with scRNA-seq data, C1-C4 cells were identified as CD19+CD20−IL-7Rα+CD93+ cells (Fig. 2c). These B cell progenitors occurred at a high frequency in the first few days after birth (until day 11), after which their frequency dramatically decreased (Fig. 2d). We also noted that B cell progenitors production ceased at about weaning age (Fig. 2d).
The lack of MHCII molecules and CD138 expression on fetal B1 progenitors distinguished them from bone marrow B2 progenitors [48,49,50]. C1-C4 cells expressed almost no detectable H2-Aa (encoding the MHCII molecules subunit), Ciita (MHCII molecules regulator), or Sdc1 (which encodes CD138) (Additional file 1: Fig. S2b), and very few C1-C4 cells expressed MHCII molecules or CD138 on their surface (Fig. 2e), suggesting a B1 cell progenitor identity. Moreover, scoring the single cells with the signatures of C1-C9 clusters revealed that C1, C3, C4 and most of C2 cells expressed B1 cell progenitor signature (Fig. 2f and Additional file 1: Fig. S2c). Pursuing this, we sorted C1-C4 progenitors from neonates and adoptively transferred them into sub-lethally irradiated Rag2−/− mice. Engrafted control adult bone marrow CD19+CD20−IL-7Rα+CD93+ progenitors mainly produced B2 cells, whereas neonate peritoneal CD19+CD20−IL-7Rα+CD93+ C1-C4 progenitors generated B1 cells (Fig. 2g).
C1-C2 cells expressed typical pro-B cell genes (Igll1, Vpreb1, Vpreb2, Dntt, and Kit), whereas C3-C4 cells had a pre-B cell signature (Cxcr4, Klf2, and Rhoh) (Fig. 2h). Comparing C1 with C2 cells, we observed C1 cells had higher expression levels of surrogate light chain genes (Vpreb1, Vpreb3) and Ig related genes (Cd79a, Cd79b, Mzb1) (Fig. 2i; Additional file 3: Table S2). C2 cells had a prominent proliferative signature marked by the expression of cell cycle regulators (Mki67, H2afz, Ran) (Fig. 2i; Additional file 4: Table S3). Comparing cells of the C3 with C4 clusters, C3 cells were clearly marked by the expression of cell proliferation regulators (Mki67, H2afz, Ran), whereas C4 cells had relatively higher expression of the anti-proliferation gene Btg1 and the Ig-associated gene (Cd79a, Cd79b, Mzb1, and Vpreb3) (Fig. 2j; Additional file 5: Table S4; Additional file 6: Table S5). Interestingly, we discovered many potential cluster-specific genes (Additional file 1: Fig. S2d, e): for example, Fcrl6, Speer2, and Selenom (in C1 cells); Bok, Dnph1, Impa2 (in C2 cells); Rgs16, Serinc5, and Il2ra (in C4 cells).
Successive developmental stages of bone marrow B2 cells can be distinguished based on the correlated expression of various cell surface markers . We investigated similar surface markers to experimentally identify C1, C2, C3, and C4 cells. Within the neonate peritoneal CD19+CD20−IL-7Rα+CD93+ cell population, the C1-C2 cells were C-kit+CD43+, and could be further divided into Ki67− C1 cells and Ki67+ C2 cells. C3 cells were Ly6d+CD25−Ki67+, and C4 cells were Ly6d+CD25+Ki67− (Additional file 1: Fig. S2f, g and 3d).
A previous Bcl11aeGFP reporter mice study showed that the transcription factor Bcl11a is highly expressed in early lymphoid progenitors . Using this mouse reporter line, we found a subset of Lin−CD19−C-kitlowSca-1low cells which expressed a high level of Bcl11a in the neonate peritoneal cavity (Fig. 2k and Additional file 1: Fig. S2h). We isolated these cells from the Bcl11aeGFP reporter mice for adoptive cell transfer assays. In contrast to engrafted control adult bone marrow Lin−CD19−C-kitlowSca-1low Bcl11a-eGFP+ cells mainly generating B2 cells, engrafted neonate peritoneal Lin−CD19−C-kitlowSca-1lowBcl11a-eGFP+ cells produced B1 cells (Fig. 2l), demonstrating that early lymphoid progenitors with B1 cell potential are present in the peritoneal cavity of neonate mice.
Single-cell analysis reveals two trajectories of B1 cell development in neonate peritoneal cavity
C5 and C6 cells expressed the known B cell precursor genes Vpreb3 and Cd93 (Fig. 2b, 3a), whereas C7-C9 cells expressed the mature B cell gene Ly86 and multiple MHCII molecules related genes (H2-Aa, H2-Ob, H2-Eb1) (Fig. 3a and Additional file 1: Fig. S3b). Interestingly, these MHCII molecules related genes were also highly expressed in C6 cells, but were very low in C5 cells (Fig. 3a and Additional file 1: Fig. S3b), enabling discernment of C6 cells from C5 cells. Finally, we noted that C7 cells expressed relatively high levels of B2 cell marker genes (Fcer2a and Ighd), while C8-C9 cells expressed B1 cell marker genes (Bhlhe41 and Cd9) (Fig. 3a and Additional file 1: Fig. S3c); and C8 and C9 cells were well demarcated by their expression of proliferation regulators (Tuba1b, H2afa, Hmgb2) (Fig. 3a).
We next performed flow cytometry and adoptive cell transfer analyses with C5-C9 cells. Consistent with our scRNA-seq data, C5 cells were identified as CD19+CD20+CD93+MHCII− cells and C6 cells were CD19+CD20+CD93+MHCII+ (Fig. 3b). C7 cells were CD19+CD20+CD93−CD23+, C8 cells were CD19+CD20+CD93−CD23−Ki67−, and C9 cells were CD19+CD20+CD93−CD23−Ki67+ (Fig. 3c and Additional file 1: Fig. S3d). Notably, transplantation assays demonstrated that C5, C6, C8, and C9 cells produced B1 cells, but not B2 cells in vivo (Fig. 3d).
We performed pseudotime analysis of our neonate B cell data using Palantir to help understand the relationship(s) between developmental stages and to ascertain a peritoneal B1 cell development trajectory . According to the obtained pseudotime values, C1 cells appear to represent the beginning stage for B1 cell development among CD19+ cells (Fig. 3e, f). Analysis of the V(D)J recombination status of both Ig heavy and light chains showed that few cells in C1-C2 had completed recombination (Fig. 3g), confirming that C1-C2 cells represented the earliest B1 cell progenitors among the detected clusters. It was notable that only a few of the cells from the C1-C2 clusters expressed Ig light chain constant-region transcripts (Igkc, Iglc1, Iglc2, Iglc3) (Fig. 3h), findings suggesting that the majority of C1-C2 cells cannot initiate light chain rearrangement .
C3-C4 cells initially expressed Igkc, and some C4 cells expressed Iglc1, Iglc2, and Iglc3 (Fig. 3h). C5-C9 cells appeared to accomplish for V(D)J combination of both Ig heavy and light chains, and express light chain constant-region transcripts (Fig. 3g, h). Interestingly, the pseudotime analysis suggested that C4 cells may differentiate directly from C1 cells (Fig. 3f). We noted that 5.6% (4/72) of C1 cells completed V(D)J combination of both Ig heavy and kappa chains and also expressed Igkc (Fig. 3g, h), suggesting that Ig heavy and light chain recombination may occur simultaneously in these cells. If so, this could result in functional BCR formation, which may bypass the pre-BCR dependent expansion and selection stages . Taken together, our Palantir analysis supports the speculation that neonate peritoneal B1 cell development apparently comprises two trajectories: a pre-BCR dependent developmental trajectory (C1-C2 to C3-C4) as well as a pre-BCR independent developmental trajectory (C1 to C4).
B1 cells exhibit profound age-related changes in their transcriptomes
In addition to our focused study of the neonate-specific subgroup of peritoneal CD19+ cells, we also re-clustered mature B1 cells (Bhlhe41+Fcer2a−Ly86+) presenting in mice of all the sampled age groups (Fig. 4a). This generated five clusters of B1 cells which we denoted as B1 cell clusters 1–5 (B1-C1 to B1-C5) (Fig. 4a; Additional file 2: Table S1). Perhaps the most obvious trend was that these clusters were differentially enriched in mice of the different age groups (Fig. 4b). That is, the large majority of the total cells from B1-C1 and B1-C2 were from the neonate mice (both the day 5 and 6 groups) (Fig. 4b). In contrast, the vast majority of the B1-C5 cells were found in the most advanced age group of mice (Fig. 4b and Additional file 1: Fig. S4a). The B1-C3 and B1-C4 cells were variously present in young adults and in the elderly mice (Fig. 4b).
Many previous studies have shown that young B1 cells exert different functions than aged B1 cells in a variety of disease contexts, including streptococcus pneumonia infection [39, 56], atherosclerosis [57, 58], and insulin resistance . We therefore used our B1 cells clusters to explore gene expression signatures characteristic of differentially aged B1 cells. Comparing B1-C5 cells against the other non-proliferating adult B1 cells, we detected 1,468 upregulated DEGs in C5 cells, including 71 cellular senescence-related genes by a Reactome  pathway analysis (Fig. 4c and Additional file 1: Fig. S4b; Additional file 7: Table S6; Additional file 8: Table S7). Notably, the known cellular senescence marker genes such as Anx5a , Cd36 , Creg1  and Plk2  were expressed higher in B1-C5 cells (Fig. 4c, d and Additional file 1: Fig. S4c). Interestingly, the surface molecule CD36 is highly enriched in B1-C5 cluster (Fig. 4d), we confirmed this while also demonstrating that CD36 can be used to easily isolate cells of this aged B1 cluster via flow cytometry (Fig. 4e, f). We anticipate that this surface marker will greatly facilitate studies of the aging-related B cell dysfunction generally and in several specific disease contexts.
B1 cell aging is associated with clonal expansion
To characterize aging-related changes in the B-cell receptor repertoire of peritoneal B1 cells we performed scBCR-seq; that is, high-throughput single-cell B-cell receptor sequencing to obtain the paired full-length variable regions of heavy (VH) and light (VL) chains that determine the antigen specificity of antibodies . The single-cell resolution data allowed us to characterize BCR diversity on B1 cells, to identify particular B1 cell clones, and to quantify their expansion based on the numbers of cells for each clone.
We first characterized the BCR diversity of peritoneal B1 cells with age. D50 metric analysis quantifying the paired VH-VL amino acid sequence diversity revealed that the BCR repertoire of peritoneal B1 cells became progressively more restricted with advancing age (Fig. 5a). This phenomenon was in obvious contrast to peritoneal B2 cells, which stably maintained a diverse BCR repertoire at each of the examined ages (Fig. 5a). Consistent with our D50 analysis, the repertoire space occupied by the top 10 clones of B1 cells progressively increased with age, whereas that repertoire space of peritoneal B2 cells was stably maintained (Fig. 5b). Clonotype counting showed that almost all of the peritoneal B2 cell clones comprised individual cells, and this remained true as mice age (Fig. 5c). In contrast, the percentage of peritoneal B1 cell clones representing individual cells gradually decreased with age (Fig. 5c), and this decrease was accompanied by a concomitant increase in the proportion of B1 cell clones comprising many cells (> 100 cells) (Fig. 5c).
Given the BCR repertoire of peritoneal B1 cells becomes more restricted with age, we asked whether the identified B1 cell clones showed preferential usage of particular VH genes, VL genes, and/or VH-VL gene pairings. We observed broad usage of both VH and VL genes by young adult and elderly B1 cells (Additional file 1: Fig. S5a). However, despite the broad VH and VL gene usages, certain VH-VL gene pairings—notably Ighv12-3/Igkv4-91 and Ighv11-2/Igkv14-126—were obviously enriched in young adult (~ 20%) and elderly (~ 24%) B1 cells (Fig. 5d and Additional file 1: Fig. S5b). In contrast, no obviously biased usage of VH genes, VL genes, or VH-VL gene pairings was detected in neonatal B1 cells (Fig. 5d and Additional file 1: Fig. S5b).
Collectively, these data support that strong enrichment for certain adult B1 cell clones reflects clonal expansion with age [40, 65]. Indeed, clonotype tracking analysis revealed a number of B1 cell clones which obviously expanded with age (Fig. 5e). Of particular note, we found that seven of nine expanded B1 cell clones expressed either Ighv12-3/Igkv4-91 or Ighv11-2/Igkv14-126 pairing (Fig. 5e), both of which are known to recognize the phosphatidylcholine (PtC) present in plasma membranes of self-cells . Moreover, B1 cells that marked by the expression of cell proliferation regulators (Mki67, H2afz, Ran) were highly enriched with anti-PtC B1 cell clones in both young adult (~ 43%) and elderly (~ 36%) mice (Additional file 1: Fig. S5c, d). The proliferation capacities of anti-PtC B1 cell clones in young adult and elderly mice also reflected their clonal expansion phenomenon.
Therefore, our scBCR-seq results for neonates, young adults, and elderly mice provide direct evidence that the BCR repertoire of peritoneal B1 cells becomes progressively more restricted with age and demonstrates that certain self-reactive B1 cell clones accumulate with age.
Mouse peritoneal cavity contains many B1 cells [3, 4, 67, 68]. The fetal peritoneal cavities of both mice and humans were discovered as an enriched source of B1 cells three decades ago [22,23,24]. Single-cell study of peritoneal CD19+ B cells from neonates, young adults, and elderly mice here demonstrates that B1 cell represent a large majority (> 80%) of the total peritoneal B cells in unmanipulated wild-type C57BL/6 adult mice under physiological condition. The single-cell data also reveal B1 cell precursor subsets in the neonate peritoneal cavity and delineate distinct B1 progenitor developmental stages and pathways.
Neonate peritoneal cavity harbors several differentially developmental stages of B1 cell precursors, including a non-proliferating pro-B cell subset (C1), a proliferating pro-B cell subset (C2), a proliferating pre-B cell subset (C3), and a non-proliferating pre-B cell subset (C4) (Additional file 1: Fig. S3d). We have identified many potential cluster-specific genes, for example, Fc receptor-like 6 (Fcrl6) that was specifically expressed in non-proliferating pro-B cell cluster. Fcrl6 was recently demonstrated to be a marker for defining biologically distinct B cell progenitors . Consistent with the observation that Fcrl6+ pro-B cells exhibited protracted differentiation and proliferation , our single-cell transcriptome analysis showed no detectable cell cycle regulators was expressed in those Fcrl6+ pro-B cells (C1). Instead, Fcrl6+ pro-B cells had a high expression of the recombination activating gene, Rag1. These findings support that Fcrl6+ pro-B cells were undergoing Rag-mediated V(D)J gene segments rearrangement that has been demonstrated to couple with cell cycle state . Moreover, our pseudotime analysis indicates Fcrl6+ pro-B cells are in the earlier development stage comparing to proliferative Fcrl6− pro-B cells (C2). Early B1 progenitor subsets characterization will facilitate the studies in understanding the transcriptional networks that control the identity and function of B1 cells.
In the classical pathway of bone marrow B2 cell development, the developing B cell precursors first rearrange the Ig heavy-chain locus in pro-B cells, which results in the formation of pre-BCR and the expansion of large pre-B cells, and subsequently differentiate to small pre-B cells that rearrange the Ig light-chain (Igk or Igl) loci [51, 71, 72]. A recent study found that the fetal liver environment could promote a ‘premature’ Igk rearrangement at the pro-B cell stage leading to direct expression of mature BCRs at the surface, the B cell development in the fetal liver could bypass a requirement for pre-BCR dependent expansion and selection stage . Taken together, these studies indicate that the bone marrow environment favors pre-BCR dependent B cell development, whereas the fetal liver supports pre-BCR independent B cell differentiation. Our pseudotime analysis suggests neonate peritoneal B1 precursors differentiation may follow two trajectories: a pre-BCR dependent developmental trajectory (C1-C2 to C3-C4) and a pre-BCR independent developmental trajectory (C1 to C4). Given that the neonatal period describes a time window that encompasses the ontogenetic transition from pre-BCR independent to pre-BCR dependent B cell development, it is reasonable to speculate that peritoneal B1 cells are likely produced by both pathways in neonates. If this speculation is correct, it will be interesting to determine which types of B1 cells are generated via a pre-BCR independent pathway versus those generated via a pre-BCR dependent pathway.
In addition, we discovered a Lin−CD19−C-kitlowSca-1lowBcl11a-eGFP+ lymphoid progenitors in neonate peritoneal cavity, which had the potential to differentiate into B1 cells, but not B2 cells upon transplantation. A recent fate-mapping study revealed that a transient developmentally restricted hematopoietic stem cell (drHSC) is the cell origin of innate-like B and T lymphocytes . The drHSCs and Lin−CD19−C-kitlowSca-1lowBcl11a-eGFP+ lymphoid progenitors are both biased toward generating B1 cells and both of them are transiently present in the early life, indicating that Bcl11a marked peritoneal lymphoid progenitors might be differentiated from the drHSCs. Collectively, these findings support the insight that the B1 cell biased lymphoid progenitors might be already intrinsically committed to becoming B1 or B2 cells by the developmental stage before Ig chain rearrangement is initiated .
Our scRNA-seq analysis detected a small number of aged B1 cells in elderly mice. The identification of the surface molecule CD36 as a specific marker for aged B1 cells should be extremely helpful, enabling them to be isolated easily with flow cytometry for basic and potentially clinical investigations. While CD36+ B1 cells account for only a small fraction of peritoneal total B1 cells in advanced age mice, those cells can induce the senescence of nearby cells in a non-cell autonomous manner as demonstrated in a recent CD36 overexpression study . CD36 is a scavenger receptor  and has been linked to the pathogenesis of atherosclerosis and Alzheimer’s disease through its recognition of modified endogenous ligands, including oxidized low-density lipoprotein (oxLDL) , fibrillar amyloid-β , and soluble amyloid-β . Cd36-deficient mice are protected from the toxic effects of atherosclerosis , insulin resistance [80, 81], and fibrillar amyloid-β , so investigating the functions of CD36 of aged B1 cells, specifically in the context of those diseases, will likely yield etiopathologically and/or therapeutically relevant biological insights.
Our scBCR-seq analysis characterized the clonal expansion of B1 cells with advancing age in a highly resolved manner. Our data for neonatal mice shows that the BCR repertoire of peritoneal B1 cells is quite diverse, with a very large proportion of cells representing individual clones. In contrast, the BCR repertoire of peritoneal B1 cells in adult mice is restricted; the BCR repertoire of these cells is dominated by a few large B1 cell clones expressing particular BCRs (e.g., Ighv12-3/Igkv4-91 and Ighv11-2/Igkv14-126). Our data for both young adults and elderly mice enable us to detect an interesting trend: the extent of this BCR repertoire restriction increases with aging. The clonal outgrowth we observed is reminiscent of the behavior of cancer cells [83, 84], considering reports that transgenic mice expressing a B1 specific BCR spontaneously developed CLL in old age [85, 86], perhaps aged B1 cells may be prone to malignant transformation.
The clonal expansion of particular B1 cells may be relevant to other diseases: aged apoE-deficient mice with established atherosclerosis are known to harbor expansion of B1 cells with the CDR-H3 sequence AGDYDGYWYFDV . Interestingly, the two B1 cell clones that were most expanded in our aged mouse samples have the same CDR-H3 sequence (Fig. 5e). Another study recently showed that the PtC binding BCR with the CDR-H3 sequence MRYGNYWYFDV can also bind oxLDL , and again B1 clones with this CDR-H3 were expanded in our aged mouse sample. The clonal expansion of similar B1 cell clones in aged mice and mice with established atherosclerosis motivates speculation about potential pathogenic links between atherosclerosis and B1 clonal expansion phenomena during aging.
In summary, our combined scRNA-seq and scBCR-seq analyses of peritoneal B cells in differently aged mice enabled the identification of distinct B1 progenitor subsets, delineation of peritoneal B1 cell developmental trajectories, and characterization of distinct aged B1 subsets, as well as profiling of the age-related clonal evolution of the BCR repertoire in peritoneal B1 cells. Our study provides a high-resolution experimental basis and resource for studying the developmental heterogeneity and origins of B1 cells, and uncovers aging-related trends in the B cell distribution and BCR repertoire.
C57BL/6 wild-type mice were used for scRNA-seq and scBCR-seq. The Bcl11aeGFP reporter mice were described in a previous study . The Rag2−/− mice (008,449) were purchased from the Jackson Laboratory. Neonates (day 5 and/or day 6), young adults (week 10 and/or week 11), elderly mice (month 15 and/or month 19) were used for scRNA-seq and scBCR-seq experiments. All mice used were maintained under a specific pathogen-free condition at the Laboratory Animal Research Center of Tongji University.
Peritoneal B cells isolation and preparation for single-cell sequencing
Mice were sacrificed by cervical dislocation. Peritoneal B cells were isolated as previously described protocol . Briefly, a small fold of skin was grasped at the umbilicus using forceps, and a small incision was made with a scissors. Then the opposing skin edges were pulled in opposite directions by fingers to expose the peritoneum of peritoneal cavity. Ice-cold DPBS (AF29584116, Hyclone) with 2% (volume/volume) FBS (S711001S, Lonsera) was injected into the abdominal cavity. The peritoneum was massaged gently and a needle with a syringe was inserted into the peritoneum to collect the fluid. Next, an incision was made in the inner skin of the peritoneum. The remaining fluid was aspirated gently with a pipette. The fluid was collected in the tube on ice, and was centrifuged for collecting peritoneal cells. The peritoneal B cells were enriched with CD19 Microbeads (130–090-301, Miltenyi Biotech). The selected individual CD19+ cells (~ 80–90% purity) from differently aged mice were used for scRAN-seq and scBCR-seq.
scRNA-seq library preparation (10 × Genomics)
scRNA-seq libraries were generated with the 10 × Chromium platform 3’ v3 or 5’ library preparation kits (10 × Genomics) following the manufacturer’s instructions. Briefly, Single-cell Gel Beads-in-emulsion (GEMs) were formed on the Chromium Controller. The reverse transcription of barcoded RNA was performed. In this step, each cell was labeled with a specific barcode, and each transcript was labeled with a unique molecular identifier (UMI). Next, cDNAs were purified and amplified. Then the libraries were constructed by performing the following steps: fragmentation, ligation of sample index, adaptor ligation, SPRIselect cleanup, sample index PCR, and SPRIselect size selection. Purified libraries were sequenced with Illumina NovaSeq.
scBCR-seq library preparation (10 × Genomics)
The samples were processed for single B-cell receptor (BCR) V(D)J clonotype with Chromium Single Cell 5’ Library & Gel Bead Kit v2 (PN-1000006). In short, Gel Beads-in-Emulsion (GEMs) were formed in channels of a chip in the 10 × Chromium instrument. GEM reverse transcription reactions were performed and the GEMs were cleaned up. Purified cDNAs were subjected to amplification before cleaning up with SPRIselect beads. BCR target enrichments were performed on full-length cDNA. Amplified cDNA underwent two rounds of target enrichment using nested primer pairs specific for mouse BCR constant regions. Single-cell BCR V(D)J libraries were prepared as recommended by the manufacturer’s user guide. The cDNA libraries were sequenced by Illumina NovaSeq.
Read alignment and gene quantification
The fastq files were aligned to the mouse genome (GRCm38 version 3.0.0), which was downloaded from 10 × Chromium website (https://support.10xgenomics.com/single-cell-geneexpression/software/downloads/latest), using Cell Ranger (3.1.0) with the ‘count’ option and default parameters. After mapping, samples were grouped into individual AnnData objects by concatenating the raw_feature_bc_matrix_h5.h5 and adding the appropriate metadata information.
Scrublet  (version 0.2.1) was used to determine the doublets in the datasets. Each cell was assigned a doublet score (scrublet_score), which was a score between 0 and 1. These predicted doublet cells as well as the scrublet_score were clustered and visualized using the UMAP method.
For quality control, we filtered the cells with higher than 13% of mitochondrial contents and cells of fewer than 500 UMIs or fewer than 200 features. Besides, we removed the cells predicted as doublets. Features expressed in fewer than 3 cells were removed in the analysis. We also removed the mitochondrial genes and ribosomal protein coding genes in the analysis. Cells expressing macrophage (Cd68), T cell (Cd3d, Cd3e), or red blood cell (Hbb-bs, Hba-a1) markers were excluded and left B cells (~ 82.8%) were used for the further analysis.
Analysis and visualization of scRNA-seq data
SCANPY  was used for data analysis and visualization. The following steps were performed in an order: data normalization, log-transformation and highly variable genes selection using the ‘seurat’ flavor and principal component analysis. SCANPY as well as seaborn  was used for data visualization of heatmap, violin plot, etc.
Identification of highly variable genes and cell clustering
We selected the highly variable genes using the highly_variable_genes function in SCANPY with parameters: min_mean = 0.015, max_mean = 10, min_disp = 0.05. 3, 904 highly variable genes were used for the data visualization. Harmony  (version 0.0.5) was used to integrate data from different samples. ‘Sample’ information in the metadata was corrected in two rounds of harmony iterations. Uniform Manifold Approximation and Projection (UMAP)  visualization and Louvain clustering  were performed based on the armony corrected latent space. The Louvain function in SCANPY was used with default parameters, while the resolution of the clustering change from 0.5 to 2. SCCAF  was used to determine the discriminative ability of the clusters.
Differential gene expression analysis of scRNA-seq data
For differential gene expression analysis, the hurdle model in MAST  (version 1.12.0) was used to model the differential expression of the cells as well as the total number of counts. The genes with a false detection rate (FDR) lower than 0.01 were used for volcano plot. The genes with log fold change > 0.1 were considered as up-regulated genes, while others were the down-regulated genes.
For B1 cell development trajectory, we used Palantir  to derive the UMAP. In this process, Palantir first calculated the principal component space based on the highly variable genes and measured the diffusion map  using the principal component space. UMAP was then calculated on the resulted diffusion map space.
The DPT function  in SCANPY was used for pseudotime analysis, while the start time point was assigned as a cell in the precursor cluster. The DPT function infered the progression of cells through geodesic distance along the graph [98, 99]. The Diffusion Pseudotime was initially introduced by Haghverdiand . The implementation within Scanpy  used a further developed version, which was able to deal with disconnected graph .
Single-cell gene signature scoring
To understand the gene expression in certain pathways, we used the score_genes function in SCANPY. The score is the average expression of a set of genes subtracted with the average expression of a reference set of genes, while the reference set is randomly sampled. B1 and B2 cell precursors signatures were from the previous study —Table S1, which listed the genes differentially expressed in Lin−CD93+CD19+CD45R−/low B1 progenitors and Lin−CD93+CD19−CD45R+CD43+ B2 progenitors, the DEGs that upregulated in day 2 neonate B1 progenitors and B2 progenitors were used.
The analysis of scBCR-seq data
The fastq files from the BCR results were aligned to the mouse vdj reference (version 4.0.0) which was downloaded from the 10 × Chromium website using Cell Ranger (3.1.0) with the ‘vdj’ option and default parameters. The basic scBCR-seq statistics was then analysed using pyvdj v0.1.2. (https://github.com/veghp/pyVDJ). For BCR V(D)J data, only cells with at least one productive Igh chain and at least one productive Igk or Igl chain were considered for analysis. Immunarch  was used for further BCR analyses, including clonotype tracking, diversity estimation, repertoire overlap, gene usage analysis and chain pair analysis. We followed the same analysis workflow as presented in the tutorials at https://immunarch.com/.
To understand the presence of clonotypic expansions, we calculated the diversity 50 (D50) , which calculated percentage of dominant unique clones that account for 50% of the total clones. A higher D50 value means a higher diversity of the CDR3 sequence.
Heavy and light chain pairing
We considered contigs identified as either heavy (Ighv) or light (Igkv/Iglv) chain variable domains. Identified variable domains were considered complete if both of the chains were found in the same cell (with the same cell barcode). We estimated the numbers of Ighv chains and the Igkv/Iglv chains, while the top 20 chain types were used for Ighv or Igkv/Iglv analysis. The jointplot function in seaborn  was used to estimate the heavy and light chain pairs and to plot the heatmap for the distributions.
Flow cytometry and cell sorting
Cells were suspended in a solution of 2% (volume/volume) FBS in DPBS. Anti-mouse CD16/32 antibody (101,320, Biolegend) was used to block Fc receptors. Then cells were stained on ice for 30 min before washing. Intracellular staining was performed according to the manufacturer protocol of the FOXP3 Fix/Perm Buffer set (421,401, 421,402, Biolegend). Fluorochrome-labeled monoclonal antibodies (clones denoted in parentheses) against B220 (RA3-6B2), CD19 (6D5), CD3ε (145-2C11), CD4(RM4-5), CD8α (53–6.7), NK-1.1 (PK136), Nkp46 (29A1.4), CD11b (M1/70), CD11c (HL3), Gr1(RB6-8C5), Ter119 (TER-119), CD20 (SA271G2), IL-7Rα (A7R34), CD93 (AA4.1), MHCII (M5/115.15.2), CD138 (281–2), Sca-1 (D7), C-kit (2B8), CD43 (S11), CD5 (53–7.3), CD23 (B3B4), Ki67 (16A8), CD25 (PC61), Ly6D (49-H4) and CD36 (HM36) were purchased from BD Biosciences, Biolegend or eBioscience. Cells were analyzed on the LSRFortessa X-20 cell analyser (BD) or sorted on the FACSArial Fusion cell sorter (BD). Data was analyzed with FlowJo software (BD).
The cell populations were stained and sorted by flow cytometry, and were subsequently injected intravenously or intraperitoneally into sub-lethally irradiated (1 × 450 rads).
Rag2−/− mice. The drinking water was added antibiotics for 2 weeks. The peritoneal cells from recipients were collected for analysis after 2–4 weeks injection.
The statistical analysis was conducted with Microsoft Excel. P values were calculated using a two-tailed Student’s t-test.
Availability of data and materials
The scRNA-seq and scBCR-seq data described in this study have been deposited in ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) under accession number E-MTAB-10081. All other relevant data are available from the corresponding author on request.
Complementarity determining region 3 of the heavy chain
Arnold LW, Pennell CA, McCray SK, et al. Development of B-1 cells: segregation of phosphatidyl choline-specific B cells to the B-1 population occurs after immunoglobulin gene expression. J Exp Med. 1994;179(5):1585–95.
Yoshimoto M, Montecino-Rodriguez E, Ferkowicz MJ, et al. Embryonic day 9 yolk sac and intra-embryonic hemogenic endothelium independently generate a B-1 and marginal zone progenitor lacking B-2 potential. Proc Natl Acad Sci U S A. 2011;108(4):1468–73.
Stall AM, Farinas MC, Tarlinton DM, et al. Ly-1 B-cell clones similar to human chronic lymphocytic leukemias routinely develop in older normal mice and young autoimmune (New Zealand Black-related) animals. Proc Natl Acad Sci U S A. 1988;85(19):7312–6.
Haas KM, Poe JC, Steeber DA, et al. B-1a and B-1b cells exhibit distinct developmental requirements and have unique functional roles in innate and adaptive immunity to S. pneumoniae. Immunity. 2005;23(1):7–18.
Lee-Chang C, Bodogai M, Moritoh K, et al. Accumulation of 4–1BBL+ B cells in the elderly induces the generation of granzyme-B+ CD8+ T cells with potential antitumor activity. Blood. 2014;124(9):1450–9.
Kreslavsky T, Vilagos B, Tagoh H, et al. Essential role for the transcription factor Bhlhe41 in regulating the development, self-renewal and BCR repertoire of B-1a cells. Nat Immunol. 2017;18(4):442–55.
Rodig SJ, Kutok JL, Paterson JC, et al. The pre-B-cell receptor associated protein VpreB3 is a useful diagnostic marker for identifying c-MYC translocated lymphomas. Haematologica. 2010;95(12):2056–62.
Mallampati S, Sun B, Lu Y, et al. Integrated genetic approaches identify the molecular mechanisms of Sox4 in early B-cell development: intricate roles for RAG1/2 and CK1epsilon. Blood. 2014;123(26):4064–76.
Van Ness BG, Weigert M, Coleclough C, et al. Transcription of the unrearranged mouse C kappa locus: sequence of the initiation region and comparison of activity with a rearranged V kappa-C kappa gene. Cell. 1981;27(3 Pt 2):593–602.
Kyaw T, Tay C, Krishnamurthi S, et al. B1a B lymphocytes are atheroprotective by secreting natural IgM that increases IgM deposits and reduces necrotic cores in atherosclerotic lesions. Circ Res. 2011;109(8):830–40.
Chen M, Yang W, Li X, et al. Age- and brain region-dependent α-synuclein oligomerization is attributed to alterations in intrinsic enzymes regulating α-synuclein phosphorylation in aging monkey brains. Oncotarget. 2016;7(8):8466–80.
Honjo K, Won WJ, King RG, et al. Fc Receptor-Like 6 (FCRL6) Discloses Progenitor B Cell Heterogeneity That Correlates With Pre-BCR Dependent and Independent Pathways of Natural Antibody Selection. Front Immunol. 2020;11:89.
Kunjathoor VV, Febbraio M, Podrez EA, et al. Scavenger receptors class A-I/II and CD36 are the principal receptors responsible for the uptake of modified low density lipoprotein leading to lipid loading in macrophages. J Biol Chem. 2002;277(51):49982–8.
Wilkinson K, Boyd JD, Glicksman M, et al. A high content drug screen identifies ursolic acid as an inhibitor of amyloid beta protein interactions with its receptor CD36. J Biol Chem. 2011;286(40):34914–22.
Hayakawa K, Formica AM, Brill-Dashoff J, et al. Early generated B1 B cells with restricted BCRs become chronic lymphocytic leukemia with continued c-Myc and low Bmf expression. J Exp Med. 2016;213(13):3007–24.
Hayakawa K, Formica AM, Colombo MJ, et al. Loss of a chromosomal region with synteny to human 13q14 occurs in mouse chronic lymphocytic leukemia that originates from early-generated B-1 B cells. Leukemia. 2016;30(7):1510–9.
Upadhye A, Srikakulapu P, Gonen A, et al. Diversification and CXCR4-dependent establishment of the bone marrow b-1a cell pool governs atheroprotective igm production linked to human coronary atherosclerosis. Circ Res. 2019;125(10):e55–70.
Finak G, McDavid A, Yajima M, et al. MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol. 2015;16:56.
We would like to acknowledge the Frontier Science Center for Stem Cell Research flow cytometry core facility and Tongji Laboratory Animal Research Center.
This work was supported by the Natural Science Foundation of China (31870877), the National Key R&D Program of China (2021YFA1100903, 2021YFA1100900), Easter Scholar (TP2018018), Fundamental Research Funds for the Central Universities (22120190230). Z.M. is supported by the startup funding at Shanghai Fourth People's Hospital, Open Targets grant (OTAR2067) and the National Key R&D Program of China (2021YFF1200900, 2021YFF1200903).
Authors and Affiliations
Shanghai Key Laboratory of Maternal Fetal Medicine, Clinical and Translational Research Center of Shanghai First Maternity and Infant Hospital, Frontier Science Center for Stem Cell Research, School of Life Sciences and Technology, Tongji University, 1239 Siping Road, Shanghai, 200092, China
Translational Research Institute of Brain and Brain-Like Intelligence and Department of Anesthesiology, Shanghai Fourth People’s Hospital Affiliated to Tongji University School of Medicine, Shanghai, 200081, China
European Bioinformatics Institute, European Molecular Biology Laboratory, Wellcome Genome Campus, Hinxton, Cambridge, CB10 1SD, UK
Guangzhou Laboratory, Guangzhou International Bio Island, Guangzhou, 510005, Guangdong, China
School of Biomedical Sciences, Li Ka Shing Faculty of Medicine, Stem Cell and Regenerative Medicine Consortium, The University of Hong Kong, 5 Sassoon Road, Pokfulam, Hong Kong
Laboratory Animal Research Center, School of Medicine, Tongji University, Shanghai, 200092, China
YL designed and carried out experiments, analyzed the data and helped write the manuscript. JW, KL, SX, ML, XL, ZZ, XX, YZ, and JP performed experiments. ZM performed all bioinformatics analyses and helped write the manuscript. PL and SG provided intellectual inputs for the manuscript. YY designed research, analyzed the data and wrote the manuscript. YY supervised the research. All authors provided critical feedback on the manuscript. All authors read and approved the final manuscript.
Figure S1. a The violin plots showing the number of genes, the number of counts and the percentage of mitochondrial genes in all examined samples. b The violin plot showing the percentage of mitochondrial genes in each sample. c The scatter plot showing that the logarithm of the number of genes and the logarithm of the number of counts are linearly correlated. d UMAP plots showing the distribution of the six samples, the number of counts, the number of genes, and the percentage of mitochondrial genes. e The UMAP plots showing the distribution of the re-clustered samples after quality control. f The UMAP plots showing the marker gene expression in macrophages (Cd68, S100a8), T cells (Cd3d, Cd3e), erythrocytes (Hbb-bs, Hba-a1) in all examined samples. The color key indicates the expression level. g The UMAP plots showing the expression of the known B cell marker genes. The color key indicates the expression level. Figure S2. a The bar plot showing the percentages of day 5 and 6 neonate samples in each of the C1-C9 clusters. b The expression patterns of the MHCII molecules related genes and CD138 encoding gene in UMAP plot. The color key indicates the expression level. c Projection of the B2 cell progenitor signature (B2 cell progenitor signature was from a previous study, see more details in Methods). The color key indicates the average log normalized expression of the B1 signature genes over the randomly sampled reference genes. d The dot plot demonstrates cluster-specific gene expression in C1, C2, and C4 cells. e, f The expression pattern of the indicated genes in UMAP plot. The color key indicates the expression level. g FACS plots showing C1-C4 cells were present in neonate peritoneal cavity. Neonates (day 5–6, pool of 12 pups from 1 litter, n = 3) were used for analysis. h FACS plots showing the gating strategy of neonate peritoneal Lin−CD19−C-kitlowSca-1low cells. Neonates (day 6, pool of 12 pups from 2 litters, n = 3) were used. Lineage markers: CD3ε, CD4, CD8α, NK-1.1, Nkp46, CD11b, CD11c, Gr1, Ter119. Data are representative of three (g, h) independent experiments. Figure S3. a The C5-C9 cells are highlighted in the UMAP projection. b, c Distribution of the indicated gene expression in UMAP plots. The color key indicates the expression level. d Table shows the identification markers for B cell subsets (C1-C9) in neonate peritoneal cavity. Figure S4. a The proportion of B1-C5 cluster in B1 cells in each mouse sample. b Bar plot showing the pathway analysis of DEGs that up-regulated (log fold change > 0.1, FDR < 0.01, n = 1, 468 genes) in the B1-C5 cells comparing to B1-C3 cells. The top 10 pathways are presenting. The gene set analysis tool GeneAnalytics (LifeMap Sciences) was used for DEGs pathway analysis. c Distribution of the aging-related gene expression in UMAP plots. The color key indicates the expression level. Figure S5. a The table lists the usages of VH and VL genes of both B1 and B2 cells from differently aged mice. b Bar plot shows the percentages of B1 cell clones that recognize the phosphatidylcholine (both Ighv12-3/Igkv4-91 and Ighv11-2/Igkv14-126 pairings) in examined mature B1 cells from mice of different age. c The expression pattern of the proliferation-related genes in the re-projected mature B1 cell UMAP plot. The color key indicates the expression level. d Bar plot shows percentages and numbers of the B1 cell clones that recognize the phosphatidylcholine (PtC) in mature proliferating B1 subset (B1-C4) in mice of different age.
DEGs down-regulated in B1-C5 cells (B1-C5 vs. B1-C3).
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.