Bone marrow and peripheral blood was obtained from 19 healthy adult donors (aged 24–86 years) with no known disease affecting the immune system and undergoing total hip replacement surgery at Guy’s Hospital, London, UK. The samples were collected with informed consent under the REC number 11/LO/1266.

The B cells were isolated and sorted as previously published (28). Briefly, BM material was removed from the head of the femur and filtered into RPMI-1640 (Sigma-Aldrich). Bone marrow mononuclear cells (BMMCs) and peripheral blood mononuclear cells (PBMCs) were isolated using Ficoll-Paque PLUS (GE Healthcare Life Sciences) according to the manufacturer’s instructions. For the BMMCs, CD19⁺ B cells were then enriched to >98% using CD19 microbead magnetic separation (Miltenyi).

Bone marrow mononuclear cells were stained using PE anti-human Ig light chain lambda (MHL-38, BioLegend), APC anti-human Ig light chain kappa (MHK-49, BioLegend), PE/Cy7 anti-human CD38 (HIT2, BioLegend), PerCP/Cy5.5 anti-human IgD (IA6-2, BioLegend), Pacific Blue anti-human IgM (MHM-88, BioLegend), APC/Cy7 anti-human CD10 (HI10a, BioLegend), and FITC CD27 (M-T271, Miltenyi Biotec). PBMCs were stained using CD19 APC (HIB19, BD BioScience), IgD PerCP/Cy5.5 (IA6-2, BioLegend), CD27 FITC (M-T271, Miltenyi Biotec), and CD10 APC/Cy7 (HI10a, BioLegend).

B cells were sorted into Sort Lysis Reverse Transcription (SLyRT) (21) buffer using the FACS Aria (BD BioSciences). B cells were sorted into four cell types: large pre-B (IgK⁻IgL⁻CD38⁺IgM⁺), immature (IgK⁺ or IgL⁺CD27⁻IgM⁺IgD⁻CD10⁺), transitional (IgD⁺CD27⁻CD10⁺), and naïve (IgD⁺CD27⁻CD10⁻) as shown in Figure 1. Due to the lytic (RNA stabilizing) nature of the sort buffer and the rarity of some of the cell populations, we were unable to check post-sorting purity. We set the collection gates well away from the FMO control gates as a precautionary measure (Figures 1B,C).

High-throughput sequencing was carried out as previously described (21, 29). Briefly, reverse transcription was performed directly on the sample immediately after sorting and then a semi-nested PCR was performed, adding multiplex identifiers (MIDs) to distinguish patients (29). High-throughput sequencing was carried out using the Roche 454 GS FLX system (LGC Genomics), and data cleanup was performed as before (29). In addition, for analysis of the CDR3 peptide sequence character, the data were cleaned to remove sequences where the CDR3 was likely inaccurate as a result of sequencing error, i.e., CDR3 regions outside the normal distribution of CDR3 lengths (1–35 amino acids for heavy chain and 1–20 amino acids for light chain) and/or sequences identified by IMGT as unproductive.

V(D)J gene assignment was carried out using IMGT/HighV-QUEST (30, 31). The physicochemical properties of the CDR3 amino acid sequences were calculated using the R package Peptides (32, 33), and clustering analysis of the Ig gene sequences was carried out using Levenstein distance on the CDR3 regions using R scripts available on our website (34).

As all of the repertoires were antigen-naïve, then true clonal expansions would be negligible. Therefore, in order to remove biases caused by PCR amplification, only unique gene rearrangements were used for this analysis. Where the clustering identified more than one related sequence, a modal sequence was used to represent the gene rearrangement. The data were stored in CSV files, and data analysis was performed using Microsoft Excel, GraphPad Prism, and R.

The frequency of each gene (both at the individual gene and at the gene family level) observed in the data was calculated for each cell subset from each donor. The frequency (in percentage) of each VDJ family combination (heavy chain) or VJ family combination (light chain) was also calculated for each cell subset from each donor. The mean values of gene combination frequencies were calculated for each cell subset, and 3D bubble plots were created using the R package plot3D (35). Statistical analysis (Mann–Whitney, Wilcoxon test, and ANOVA, with post-test analysis where appropriate) was performed using R or GraphPad Prism.

The physicochemical properties of CDR3 regions at heavy and light chains were compared between different cell types. These properties consisted of length, hydrophobicity indicated by GRAVY index (36), Boman index (37), molecular weight (Mr), isoelectric point (pI) (38), aliphatic index (39), frequency of amino acid classes in the CDR3 region, and Kidera factors (40). The R package lem4 (41) was used for fitting and analyzing the mixed model of our data, describing the fixed-effect (cell types) and the random-effect (patients) in a linear predictor expression. The likelihood ratio test was calculated with the statistical method ANOVA to estimate the statistical significance between populations, i.e., a pair of cell subsets.

Principal component analysis (PCA) and clustering, using Minkowski distance, were applied to the Kidera factors and gene usage frequencies from the CDR3 data as follows. First, the mean values of the Kidera factors and gene usage frequencies were computed for each donor. Second, the mean values and frequencies of all donors were grouped and then analysed by PCA and clustering.

Principal component analysis was performed using the prcomp function in R. The Minkowski distances (with power of 4) were calculated using dist(method = “minkowski”) function in R based on all CDR3 properties. The distances were then plotted with dendrograms (trees) using the dendrapply function in R.

Randomise datasets were generated by randomly shuffling the sequences across four cell subpopulations. PCA analysis was then performed to be compared with the original dataset in order to show that our observations of differences between cell subpopulations were not random events.

Peripheral blood mononuclear cells were stained with FITC anti-human CD14 and APC anti-human CD3 (clone M5E2 and UCHT1, respectively), and a population of enriched B cells (CD3⁻CD14⁻) was collected into 50% FCS (Biosera) and 50% RPMI-1640 (Gibco). The CD3⁻CD14⁻ enriched B cells were labeled with a rhodium intercalator (Rh103, DVS Sciences) followed by intracellular and extracellular staining with a panel of 30 different metal-tagged antibodies (DVS Sciences, BD BioSciences, and BioLegend). Cells were fixed, iridium stained (Ir193, DVS Sciences), and normalization beads (DVS Sciences) were added before analysis on the mass cytometry system (DVS Sciences). Between 1 and 5 × 10⁵ stained cells were analysed per sample.

Data were normalized and files were concatenated and cleaned up to remove debris (by gating on cell length and DNA⁺ cells), to exclude normalization beads (Ce140⁻ cells), to positively select intact cells (Ir191⁺Ir193⁺), to positively select live cells (Rh103⁻Ir193⁺), and to identify CD19⁺ and/or CD20⁺ B cells. CD38^hiCD24^hi B cells were identified and exported as a new group prior to performing SPADE (Spanning-tree Progression Analysis of Density-normalized Event) analysis (42). SPADE analysis groups cells into “nodes” based on the expression of all 30 markers to produce a two-dimensional tree. Using a color coded expression scale, the nodes in the tree were manually grouped into larger “bubbles” to collect together nodes, and therefore cells, which had similar expression, i.e., all those with high IgM expression were grouped together in one bubble.

Pre-B (large pre-B) and immature B cells, from BM samples, and matched transitional and naïve B cells, from PB samples, were sorted (Figures 1B,C) prior to high-throughput sequencing using an IgM-specific constant region primer. Both the heavy and light chain (kappa and lambda) Ig genes were amplified with a total of 96,593 heavy and 49,101 light chain sequences generated after initial data cleanup. These B cell populations are all thought to be exogenous antigen-naïve and therefore will not have been activated to undergo somatic hypermutation and expansion. We do not see evidence of somatic hypermutation in the gene sequences (data not shown), and therefore, we have assumed that any sequences with the same CDR3 region arise from PCR duplication. Therefore, only one example sequence of any unique gene rearrangement was used in this analysis, resulting in 39,577 heavy chain and 42,542 light chain sequences grouped by donor and cell type. Sequencing error does not substantially affect the assignment of germline Ig genes to the sequences; however, for the CDR3 peptide analysis we further removed sequences where the peptide sequence may be inaccurate due to sequencing error. This resulted in 29,074 heavy chain and 29,128 light chain sequences (Supplementary Tables). Sequences can be accessed on the National Center for Biotechnology Information’s Sequence Read Archive in raw format (BioProject: PRJNA39946; Sequence Read Archive accession: SRP081849) or in processed format with metadata at www.bcell.org.uk.

Heavy chain V, D, and J family usage did not show any significant differences in repertoire between pre-B and immature cells from the BM. There were, however, significant differences between these BM cells and the peripheral transitional and naïve cells (Figure 2). IGHV3 family genes are the most predominant genes in the human peripheral repertoire. It was interesting that in the BM this was particularly the case, with IGHV3 cells actually being removed from the repertoire during B cell maturation: there is a highly significant >13% decrease in the use of IGHV3 family genes in naïve cells with small increases in all other families to compensate (Figure 2A). Naïve cells also showed a significantly decreased use of IGHJ6 and, together with transitional cells, a >6% reduction in use of IGHD2 family genes.

Since we had expected that peripheral transitional cells would fall between immature BM cells and peripheral naïve cells in the development pathway, and that any changes in repertoire we saw would reflect this, we were surprised to see that this was not always the case. There was a significant 5% increased frequency of IGHD3 family usage in transitional cells compared to all other cell types. Furthermore, there was a significant >9% increase in IGHJ6 usage, compensated for by decreases in IGHJ3, 4, and 5 usages, in transitional cells compared to all other cell types. This is reflected in the different size of bubble V3D3J6 in the bubble plots (Figure 2B). The different repertoire of transitional and naïve cells compared to the BM cells (p < 0.05, Wilcoxon) and compared to each other (p < 0.001, Wilcoxon) is illustrated by a PCA analysis of gene family usage (Figure 2C).

In contrast to the heavy chain repertoire, the light chain gene family repertoire does not distinguish between cell types. There are no significant changes in kappa family usage (Figure 3A). Some differences were seen in lambda families (Figure 3B). The IGLV2 family usage is significantly increased by 10–15%, at the expense of all other families, and IGLJ1 family usage is significantly increased by 2–5%, at the expense of IGLJ3. As a result of this, an ANOVA analysis of the combinatorial lambda family repertoire showed a significant difference between the immature and the transitional and naïve stages of development (p < 0.001) (Figures 3C,D). However, clustering by PCA showed that any differences in light chain VJ gene usage were not strong enough to be able to distinguish between the different cell types (Figure 3E). Nor were there any obvious differences between the different cell types in lambda CDR3 amino acid sequence, since PCA of the Kidera factors to assess the physicochemical character of the CDR3 did not distinguish between the groups (Figure 3F).

As the above analysis of gene family repertoire indicated that there were repertoire changes between cell types, we analyzed all the genes individually to check if we had missed any significant gene selection due to the averaging effect of looking at the family level (Figure 4). Not all the IGHV3 family genes are decreased in naïve cells compared to BM cells. While there are significant decreases in IGHV3-15, IGHV3-30, and IGHV3-33 in particular, IGHV3-9 is actually increased (Figure 4A). Other notable increases occur in the two main IGHV1 family genes: IGHV1-18 and IGHV1-69, and in the IGHV6 gene. The IGHD2 family decreases are contributed by IGHD2-15 and IGHD2-2, and while the compensatory increase in other IGHD genes seemed unremarkable across the board, IGHD1-7 and IGHD4-17 did show significant differences (Figure 4B). In spite of the significant change in IGHD3 family use in transitional cells, this did not show up at the individual gene level, implying that the increase occurs throughout the IGHD gene family. Despite the lack of significant changes in IGK family repertoire, there was a small (~3.8%) but significant increase in IGKV3-11 gene use in naïve cells compared to immature cells. This appeared to be at the expense of small (<3%) decreases in IGKV3-20 and IGKV4-1 genes. The increase in IGLV2 family during development seemed to be mainly due to significant increases of 12.8 and 7% in IGL2-14 and IGL2-23, respectively (Figure 4C).

There is a certain amount of interindividual variation that occurs in these analyses, but the trends for selection of these genes in the repertoire are consistent, as illustrated in Figure 4, where the individual donors are shown separately for genes that are removed from the repertoire (Figure 4D) or that are increased in the repertoire (Figure 4E) during early development.

Although much of the CDR3 region is comprised of contributions from the individual IGHV, IGHD, and IGHJ genes, reflecting some of the repertoire selection effects that are captured in the analysis above, the actual amino acid sequences encoded by CDR3 varies tremendously even within the same VDJ combinations. In addition to the direct effects of endonuclease action on the genes, and N region insertion by terminal deoxynucleotidyl transferase, the reading frame of the IGHD region can also vary. Since the CDR3 region encodes a crucial part of the antibody binding site, and key functional aspects of its structure are dependent on the primary sequence (43), we also analyzed the biophysical characteristics of the CDR3 amino acid sequence. Initially we used Kidera factors, which are a set of 10 orthogonal factors that encapsulate information from ~140 different measurable biophysical characteristics of peptides. The data from PCA analysis of the CDR3 Kidera factors are in accordance with that for the VDJ gene analysis, showing that the characteristics of pre-B and immature cells are found in overlapping clusters (Figure 5A). Naïve cells and transitional cells, however, form separate yet non-overlapping clusters. The data from heavy chain CDR3 Kidera analysis separate the groups of cells better than the gene usage data, with 30% of the data contributing to PC1. To elucidate which characteristics were mainly responsible for the differences, we analyzed some of the most common ones individually. The numbers of charged, basic, and aromatic amino acids in each sequence, and the sequence Boman index, were significantly increased in naïve cells compared to pre-B cells (Figure 5B). Conversely, the number of small amino acids per sequence, the hydrophobicity (GRAVY index), aliphatic index, and overall length of sequence were all disfavored characteristics that were removed from the repertoire during development (Figure 5C). Interestingly, the selection on the size of CDR3 region did not seem as strong in the older donors as it did in the younger ones (Figure 5D).

The heavy chain gene and CDR3 PCA analysis (Figures 2C and 5A) indicated that transitional cells, in addition to being distinctive from pre-B cells and immature cells, also had a different repertoire to naïve cells. We used cluster analysis (based on Minkowski distances) to investigate the relationships further, which confirmed, by both VDJ usage (Figure 6A) and Kidera factors (Figure 6C), that transitional cells have a different repertoire to the other cell types. Naïve cells formed a sub-branch of the cluster containing pre-B and immature cells suggesting that the naïve repertoire is more similar to the BM cells than to the transitional cells. Clear examples of individual genes where the usage in transitional cells differs from the rest of the cells can be seen in Figure 6B, and biophysical characteristics showing the significantly different character of the heavy chain CDR3 in transitional cells are shown in Figure 6D. Since this subset of cells has been reported to contain Bregs, as well as being the precursor to naïve B cells, we investigated the heterogeneity of the population by mass cytometric analysis of surface markers. Although the population is small, it does appear to contain a number of different potential subpopulations, as illustrated by the IgM SPADE plot in Figure 6E.