Heterogeneous clonal diversity across GCs

We performed LCM on 10 individual GCs from a human LN via using five serial tissue sections (Fig S10) and we applied Rep-seq to obtain sequencing data of the IGHV of the BCR, obtaining ∼100,000 sequences per sample (Fig 1A and E). Two replicate PCR amplifications were performed for each GC and analyzed to verify the reliability and reproducibility of our PCR amplification and NGS approach, and to quantify the variance intrinsic to our measurement protocol, as multiplex PCR approaches may be subject to amplification bias, especially when PCR substrate is limiting. To that end, we studied the effect of the DNA concentration on PCR reproducibility by comparing the replicate analysis performed on only one isolated tissue section per GC with the replicate analysis from five (serial) tissue sections per GC. We observed much lower Sørensen–Dice similarity indices when we used only one tissue section, indicating variegated bias in each individual amplification and sequence run when PCR substrate is limiting (Section 1 in Supplemental Data 1). However, when using ∼fivefold more PCR input, replicate analyses showed excellent concordance (Fig 1B), suggesting that PCR bias has a minor impact on our results. Moreover, mutation analysis revealed that all samples were very similar in terms of mutation rate and the nature of these mutations (Section 2.1 in Supplemental Data 1). Following previous conventions (3, 15, 21 Preprint), sequences were grouped together into clones by shared V, J gene segments and CDR3 length, and more than 84% junction nucleotide sequence identity, as optimized from the distance to the nearest distribution model (22) (Section 3 in Supplemental Data 1). We analyzed the similarity between the clonal repertoire of each sample and observed a high similarity between samples from the same GC (Fig 1B and C), confirming the reliability of our sequencing protocol, showing a low overlap between IGHV sequences obtained from different GCs (<0.1), suggesting that GCs are distinct evolutionary environments and thus relatively independent from each other.

We then classified sequences into three categories according to their frequency, corresponding to the dominant clone (i.e., most abundant clone), an expanded clone (frequency >1%), or a non-expanded clone (frequency <1%, typically involving less than 1,000 sequences belonging to that clone in our data). As highlighted in Fig 1D, the clonal abundance within different GCs is heterogeneous.

To quantify the heterogeneity between different GCs, we computed the clonal diversity in terms of dominance (proportion of the most abundant clone), evenness (homogeneity of the clonal abundance), richness (number of clones), and Shannon entropy, an alternative measure of diversity, that is, less sensitive to singletons (clonal group with a unique sequence). In agreement with previous mouse models (3, 8) that identified heterogeneous clonal dominance across GCs, we show that the clonal diversity within each GC takes a wide range of values (Fig 1C and F). As an example in our setting, the clonal dominance ranged from 5% (GC9) to 30% (GC10). To study the role of sample size in the GC heterogeneity quantification, we paired each sample with an artificial sample of equivalent size, obtained from randomly selecting sequences from all the GCs combined. As highlighted in Fig 1F, differences related to sample size are not statistically significant compared with the diversity variation across samples. These findings highlight the importance of studying individual GCs and prove that the heterogeneity is not a result of sub-sampling. In Fig S11A–D, we show that the variability in the diversity metrics between samples is consistent across different clonal identification methods, where the same conclusion can be obtained when quantifying diversity in terms of CDR3 abundance only.

Heterogeneous V gene repertoire usage across GCs

Analyzing the abundance of V genes in each sample reveals that GCs differ significantly in terms of their V gene usage (Fig 2A–C and Section 4 in Supplemental Data 1). Interestingly, three genes in particular stand out with regards to their frequency across all samples: IGHV1–2, IGHV2–5, and IGHV1–18 (Fig 2E). Although these genes were shown to be the most frequently found in the peripheral blood of adult patients suffering from end-stage renal disease (23), they were still found in a higher abundance in our data than in other LNs and the bone marrow of healthy donors analyzed by a BIOMED-2–based multiplex PCR amplification and sequencing approach similar to ours (24) (Fig 2E), suggesting that these V genes may have been positively selected during the GC reaction in the reactive lymph node used in our study. This is further supported by the fact that these genes are found at a higher frequency than the most abundant V genes in non-functional BCR sequences (Fig 2F), which are not subjected to selection. The same figure with functional and non-functional V gene abundances for the same V gene label is provided Fig S13.

Next, we analyzed the V gene usage in the top 15 dominant clones from each individual GC. We found that the frequencies of the IGHV1–2, IGHV2–5, and IGHV1–18 genes varied significantly between the GCs (Fig 2D). As an example, in GC6, the IGHV1–18 or IGHV2–5 genes were not found in the top 15 clones, whereas they were present in seven and five of the top 15 clones from GC9, respectively. Similar results were obtained when using different thresholds for dominance (Section 4 in Supplemental Data 1). Moreover, we analyzed the clones most abundantly encountered in our analysis and their features with two different sequence-grouping algorithms. We used the TRIP tool (25) to group the sequences with identical CDR3 amino acid sequences into clones. In parallel, we used the hierarchical agglomerative clustering algorithm to cluster the junctional sequences (26 Preprint), which was previously proven to be effective at identifying clones (27). This approach requires identical V and J gene segments, identical CDR3 lengths, and more than 84% junction nucleotide sequence overlap to assign two sequences to the same clone (27) (see the Material and Methods Section “Grouping clone sequences” for more details). Both algorithms result in a similar identification of the most abundant clones. We observed that the clones most abundantly encountered in our analysis exhibit common features in terms of CDR3 length and V, D, and J gene usage (Section 5 in Supplemental Data 1). These data reveal different V genes expansions in individual GCs and convergent evolution of the expanded and most abundant clones.

Heterogeneous clone functionality across GCs and SHM-induced crippling mutations

We classified the sequences into four categories based on the identification of a frameshift or a stop codon. The categories are the following: (i) functional sequences, (ii) out of frame non-functional and in frame non-functional due to an early stop codon induced by (iii) V(D)J recombination or (iv) SHM. Out of frame sequences were defined as having an out of frame junction, with a frameshift detected by IMGT-V-Quest (28). In frame non-functional sequences were assumed to arise either from V(D)J recombination when having at least one stop codon in the IMGT N-region or from SHM otherwise. More specifically, the BCR junction is formed when the germline V, D, and J genes are associated during V(D)J recombination. During this process, additional nucleotides are inserted between the V, D and D, J genes, which are referred to as the N-regions (29). Consequently, we assumed that stop codons in N-regions are mostly derived from V(D)J recombination errors.

We observe that most of the sequences is functional (≈90%), and 78% of the non-functional sequences result from a frameshift. The remaining sequences are roughly equally divided between SHM-induced and V(D)J recombination–related stop codons. The distribution appears consistent between GCs, with the proportion of non-functional sequences ranging between 10% and 25% across samples (Fig 3B). The frequency of non-functional sequences was found to be relatively similar across samples from the same GC, except for GC5, where non-functional sequences were found in a low abundance in one of the two NGS samples. The high degree of similarity of the non-functional sequences in GCs 4–6 and GCs 2–8 (Fig 3E) was related to shared clones in these GCs, which also carried the non-functional allele, supporting the notion that these are due to V(D)J recombination errors. Studies in mouse models focusing on the role of crippling mutations during SHM identified a limited number of stop codons in the GC with a range of 2–5%, where the crippling mutations were strikingly enriched in B cells, expressing low surface BCR levels, being seven times more abundant (17.2%) (14) compared with our current analysis of a human lymph node.

We computed the SHM mutational spectrum by regrouping the position-wise mutation frequencies across all sequences. Our data highlighted an expected increase in the mutation frequency in the CDR regions, but no significant differences between functional and non-functional sequences were found in this context (Fig 3C). The same comparison cannot be performed at the clonal level because each clone has different V or J genes and thus the mutation positions cannot be compared across different clones. Furthermore, we also studied the SHM spectrum with regards to the type of mutations induced by AID in each GC for the functional and the non-functional clones. Mutation analysis revealed that all samples were very similar in terms of mutation rate and the nature of these mutations (Fig S2A). Interestingly, no difference between functional and non-functional clones regarding the SHM spectrum was observed suggesting that the SHM process seems to be stochastic. However, the selection pressure of the functional dominant clones was assessed by studying the replacement to silent mutation (R/S) ratio in CDR for the functional expanded clones versus the singletons (Fig S2B). The R/S ratio indicates that replacement mutations are being selected over silent mutations and it positively correlates to the selection pressure undergoing in a GC environment. We observed increased R/S ratio for the functional dominant clones compared with the “singleton” sequences, verifying the difference of selection pressure between the expanded and the non-expanded clones (Fig S2C).

Pairing between functional and non-functional alleles

The expansion of non-functional IGH alleles in the GC may be derived from selected B cells that harbor both a functional and a non-functional allele due to V(D)J recombination errors (30, 31) (Fig 3A). Based on this notion, we paired each non-functional IGH allele with a functional IGH allele by assuming that equal/similar abundance indicates that they come from the same B cell. We show in the obtained mapping that the paired functional and non-functional alleles were observed in all the GCs, consistently in both replicates. Furthermore, the frequency of the functional alleles that were paired with non-functional ones varies between individual GCs (Fig 3D).

Identifying non-functional clones derived from SHM

We performed functionality analysis also at clonal level across GCs and separated the clones into three categories. We defined a clone as functional if it contained more than 95% functional sequences (F), non-functional (NF) if it contained less than 5%, and F and NF otherwise, that is, clones consisting of both functional and non-functional sequences. The F and NF clones corresponded to the 4% of all the clones studied, without significant correlation with their abundance (i.e., being F and NF seems to be independent of being selected). All of them were in frame and likely caused by SHM because the stop codons were outside of the junction region. We found at least 5 F and NF clones with 1,000 different unique sequences and investigated these in more detail. After inference of their phylogenetic trees, we can observe a separation between functional and non-functional branches (Figs 3F and S13). Thus, we can check for mutations that were selected through affinity maturation, as performed previously in mouse models (15). For example, we observe several selected mutations of the functional sequences in the FWR2 region of the studied F and NF clone (Fig 3G). Still, we did not identify specific mutations consistently selected across clones (Fig S13). On the other hand, our analysis revealed that specific crippling mutations occurred independently in different F and NF clones, when they shared the same V genes (Table S1). As an example, the mutation Y59* occurred consistently in five independent clones with the IGHV3 gene.

Regarding the timing of these crippling mutations, the distance between the non-functional sequences and the root in the inferred tree can inform us about the number of mutations the BCR sequences underwent before the crippling mutation occurred, which correlates with the number of cell divisions (roughly one mutation per two cell divisions (9)). It was estimated from our generated trees that an average of 3.6 mutations occurred before the crippling mutation (Fig S13).

Expansions of specific V genes are frequently involved in stereotypic rearrangements found in B-cell malignancies (IGHV1–8, IGHV1–2, IGHV3–23, IGHV4–34 (32, 33)), autoimmune diseases, and infectious diseases (IGHV4–34, IGHV5–51, IGHV1–69, IGHV1–46 (34)). We observed that most of the F and NF clones with crippling mutations (22/35) use some of those genes. More specifically, the genes IGHV1–2, IGHV3–23, IGHV4–34, and IGHV2–5 were used by the F and NF clones. The SHM can act as a double-edged sword for the organism because it is necessary for an effective immune response, but at the same time it introduces mutations that can induce the recognition of self-antigens and consequently can lead to autoimmune disease. The GCs not only rely on the selection of the antibodies with the highest affinity but also on autoreactivity checkpoints that are needed for the counter-selection of B cells that can bind to self-antigens (1). As studies in mouse models demonstrated that self-reactive GC B cells are counter-selected or inactivated by SHM (35), our finding supports this hypothesis of counter-selection mechanism of self-reactive B cells in the human GC. Still, the role of negative selection in the GC is not yet completely clear, despite the deepened knowledge obtained from recent studies using mouse models (14, 36).

B-cell reactivation in the GC

Several studies in both rodents and humans by us and others have shown various frequencies of shared clones between different GCs (3, 19, 20), which might be indicative of B-cell reactivation in ongoing GC reactions (8). We investigated the number of shared clones in a human LN. To increase the robustness of our analysis, we considered a clone to belong to a GC only if it was consistently identified in both NGS replicates. We found 10.8% (396/3,650) of functional clones to be present in at least two GCs (Fig 4B). This is in line with a recent study in the mouse, where ∼10% of the clones were found to be shared between individual GCs induced by CGG immunization (20). Moreover, some of the shared clones exhibit common features with other shared clones across GCs, regarding V, D, and J gene usage and CDR3 length (Section 5 in Supplemental Data 1). Shared clones between different GCs indicate a frequent B-cell recirculation from one GC to another (Fig 4A). Nevertheless, clones that were expanded in more than two GCs were rare and most of the shared clones were expanded only in one GC, with only 5% of clones found with dominance >0.1% in more than one GC (Fig 4B). The fact that shared clones exhibit different degrees of dominance in different GCs fit with the notion that GCs are distinct competitive environments, leading to clonal expansions varying by several orders of magnitude across GCs. In fact, a recirculating B cell enters different GCs at a disadvantage compared with other already established dominant clones, where its initially low abundance makes it more difficult to compete for T cell help (9).

Interestingly, the distribution of shared clones among GCs fits a Poisson distribution well (Fig 4B), suggesting that the reactivation of a B-cell in a different GC is a memoryless stochastic process where the past evolutionary history of the B cell does not play any role. According to the Poisson distribution and denoting as λ the occurrence rate (in our case the rate of B-cell reactivation), the number of occurrences N at the exposure time t is a random variable of mean μ = λt and distribution: P[N=n]=µne−µn!.(1)Fitting Eq. (1) to the observed distribution gives μ = 0.25. From this, we can compute the probability that a given clone seeds at least one other GC during the whole GC reaction P(N ≥1) = P(N ≥0) − P(N =0) = 1 − e^−μ = 1 − e^−0.25 ≈ 22%.

We can also estimate the probability of seeding per clone and day. When using 20 d as a typical lifetime of a GC, generated in response to immunization with a simple hapten-protein conjugate such as NP-CGG (3), we can estimate the GC B-cell reactivation rate to be λ ≈ 0.2520 = 0.0125 seeding per day per clone. However, it is worth keeping in mind that the calculated reactivation rate is only an upper limit, as the GC lifetime most likely is much longer in chronically activated LNs. For instance, lifetimes of up to 100 d have been observed for some viral infections (37), and GCs formed in response to synthetic antigens can persist for up to 1.5 yr (38). Such extended GC lifetimes would result in much smaller reactivation rates. As an example, if the data would have been collected at day 200 of the GC reaction, the estimated reactivation rate would decrease to 0.1% (Fig 4A). We are interested not only in how many clones (sequences sharing the same V, J gene segments and CDR3 length, and more than 84% junction nucleotide sequence identity) are shared across GCs but also in how many identical CDR3s are observed between different GCs. We observed that only 2% (38/1,885) of the CDR3s are shared among GCs, which is considerably less than for clones. This indicates that the CDR3 typically continues to mutate after B-cell recirculation in the new GC environment (Fig 4F). To gain insight into the parallel evolution of a given clone in different GCs, we inferred phylogenetic trees from clones found in multiple GCs. In Fig 4C, we show the phylogenetic tree inferred from clone IGHV1-2_63_IGHJ6, which was found expanded with dominance >1% in two different GCs. The hierarchical reconstruction of the tree clearly separated the two GCs, where sequences in the same GC were more similar to each other than sequences in different GCs (Fig 4D). In this example, sequences in GC8 are generally closer to the root (unmutated V gene, J gene, and consensus CDR3) than GC1 in terms of sequence similarity (Fig 4E), suggesting that GC8 yielded a clone that initiated GC1. This apparent temporal separation (Fig 4D and E) was also observed in the Brainbow mouse model after 20 d of subcutaneous immunization with alum-adsorbed CGG (3, 39), suggesting that mutations may be selected through affinity maturation at different rates and through different evolutionary processes in different GCs. Such a clonal tree pattern was found to different extents in other shared clones and is thus not anecdotal (other shared clones are provided in Fig S14).

Convergent epitope reactivity across GCs

We studied the structural similarity of dominant clones by considering their paratope structure via paratope modeling. To study the paratope similarity of dominant clones, we assigned to each clone in each GC a representative sequence, defined as the most abundant sequence within that clone. Then, we combined three metrics, the CDRs similarity, Paratype (40), and Ab-Ligity (41), to identify clones that are likely to bind similar epitopes from their heavy chain sequences (Section 6 in Supplemental Data 1). We refer to those predicted clones as clone pairs. If the similarity of the clone pairs is above the optimized threshold for two of the three aforementioned metrics, the clone pair is predicted to have evolved towards binding the same epitope. We show a few illustrative examples of predicted clone pairs in Fig 5A. Although all pairs have a relatively low CDR3 sequence similarity (ranging from 0.63 to 0.74), the pairs exhibit identical (or nearly identical) CDR1 and CDR2, and similar paratope (>0.76) and high structure similarity. The first two clone pairs in the Fig 5A (first and second column starting from the left) could be the same clone misclassified as a different clone because of too many newly accumulated mutations in the CDR3. Alternatively, they could have originated from highly similar naive B cells. The fact that they exhibit differences in the nucleotide sequence of the N region of the V and J genes suggests that they are different clones exhibiting common epitope reactivity. The last two clone pairs (third and fourth column starting from the left) have different V and J genes, respectively, with a 4(7) nucleotides difference between the two V(J) regions and considerable CDR3 dissimilarity, but their paratope and CDR3 loop structure were predicted to be highly similar (Fig 5A).

We classified the predicted clone pairs with common epitope reactivity within one GC and across GCs into the following categories:

• • Same clone present in different GCs: this is the most straight forward way to study GC convergent evolution (blue curve in Fig 5C).

• • Same V and J gene, but CDR3 similarity <0.84: This category needs to be carefully analyzed because some clones could be misclassified as different because they have accumulated too many mutations in the CDR3 and their CDR3 similarity is thus below 84%, the threshold we used for clone identification in our study (green curve on Fig 5C and D).

• • Different V or J gene: This category is based on literature showing that different gene rearrangements can yield convergent epitope reactivity (40, 42). For this category, however, there is a chance that some pairs are misclassified because of a misalignment in the V or J genes, which is especially likely to happen when their CDR3s are nearly identical (yellow curve on Fig 5C and D).

We also investigated the number of predicted clone pairs within and across GCs based on each of these categories. We considered multiple dominance ranking thresholds, ranging from the top 10 to the top 300 clones in each GC and showed the predicted number of identified clone pairs in each of the three categories with an increased number of dominant clones (Fig 5C and D). We observe an increased number of predicted clone pairs across and within GCs for all three categories, with the number of identified clone pairs approaching 1,000 when the top 300 dominant clones per GC are considered (red line in Fig 5C and D). Most of the predicted clone pairs come from the second category (same V and J genes and CDR3 similarity <0.84, depicted by the green curve in Fig 5C and D).

Regarding the paratope commonalities between dominant clones across and within GCs, we characterized paratope similarity by studying the paratope distance between clone pairs, defined as (paratype distance + CDRs distance)/2. We show that the distribution of paratope distance is very similar across and within GCs, with only a slightly lower average paratope distance within GCs than across GCs (Section 7 in Supplemental Data 1). This finding suggests that different evolutionary forces can lead to convergent immune responses.

IGH next-generation sequencing/Rep-seq

In brief, LCM of the individual GCs from human LN was combined with the extraction of gDNA and an NGS approach to study both the functional and the non-functional clones.

Tissue preparation

Patient frozen LN tissue from the pathology department of the Amsterdam UMC hospital was used in this study. (LN studied was a cervical LN resected out of a 46-yr-old woman suffering from chronic sialadenitis. Patient material was obtained according to the ethical standards of our institutional medical ethical committee, as well as in agreement with the 1975 Declaration of Helsinki as revised in 1983). The LN was fresh-frozen in liquid nitrogen shortly after surgical resection. Five serial cryosections (10 μm per cryosection, CryoStar NX-70; Thermo Fisher Scientific) of the LN tissue on PEN slides (1 mm; Zeiss) were used to obtain higher DNA concentrations. Immunostaining was performed to visualize the GCs. The PEN slides were incubated with drops of hematoxylin (KLINPATH; VWR Life Science) for 1 min at RT, after a quick and gentle wash with demi water and tap water, they dried overnight at RT. In parallel, H&E staining of the frozen tissue was performed on normal (TOMO) slides at Amsterdam UMC diagnostics.

LCM of individual GCs

After tissue preparation, we used the objective 5x of the Leica PALM-LMD6 to draw and laser capture individual GCs from five serial cryosections of the human LN studied. The laser capture was performed five times for each GC in the same 0.2 ml Eppendorf tube.

gDNA extraction from individual GCs and DNA concentration

2 μl 1x proteinase K (recombinant PCR grade, 25 mg; Roche) dissolved in 2.5 ml TE (10 mg/ml) and 20 μl 1x proteinase K lysis buffer 50 mM Tris–HCl (pH 8, 5-Enzo), 100 mM NaCl (Sigma-Aldrich), 1 mM EDTA (Sigma-Aldrich), 0.5% Tween-20 (Bio-Rad), 0.5% NP-40 (Sigma-Aldrich) were added in the tubes used for the LCM and the samples are incubated at 56°C overnight. The proteinase K was inactivated by incubation for 5 min at 95°C. The Qubit dsDNA HS Assay kit (Invitrogen) was used to quantify the extracted gDNA, according to manufacturer’s protocol.

Multiplex-nested PCR protocol

The gDNA extracted from individual GCs was used for the multiplex-nested PCR, consisting of two rounds. Round 1: amplification of the leader genes region (LD). Round 2: amplification of the FR1 region with primers containing Illumina adapter overhang nucleotide sequences. The input used for round 2 was the round 1 PCR product. We used technical replicates to check the method reproducibility (duplicates of PCRs for the same DNA/GC).

For round 1: PCR amplification of gDNA from individual GCs (2 ng input) was performed with 1 μl of JH reverse primer (10 μM, provided by MERCK) and 1.8 μl of LD forward primer set pools (10 μM per primer, provided by MERCK) using 25 μl of Q5 Hot Start Master Mix (2X) (New England Biolabs) for a total volume reaction of 50 μl. For round 2: PCR amplification of the round 1 PCR product was performed (1 μl) with 1 μl of JH reverse primer (10 μM, provided by MERCK) and 1 μl of FR1 forward primer set pools (10 μM per primer, provided by MERCK) using 25 μl of Q5 Hot Start Master Mix (2X) (New England Biolabs) for a total volume reaction of 50 μl.

The following PCR program was used for both round 1 and round 2: 5 min at 95°C; three cycles of 5 s at 98°C and 2 min at 72°C; three cycles of 5 s at 98°C, 10 s at 65°C, and 2 min at 72°C; and 25 cycles of 20 s at 98°C, 30 s at 60°C, and 2 min at 72°C; with a final extension cycle of 7 min at 72°C on a Biometra T-ADVANCED PCR machine. Primers and adaptors used in round 1 and round 2 are provided in Fig S10C.

Round 2 PCR product purification

Electrophoresis of the round 2 multiplex-nested PCR product was performed on a 1% agarose (Molecular grade; Bioline) gel, followed by gel excision of the 400 bp PCR product with Macherey-Nagel PCR clean up and gel extraction kit, according to manufacturer’s protocol.

NGS library preparation and sequencing

The concentration of the gel excised round 2 multiplex-nested PCR product was quantified by using the Qubit dsDNA HS Assay kit (Invitrogen) according to manufacturer’s protocol. 100 ng of this PCR product were used for the preparation of the 16S metagenomic sequencing library for the Illumina MiSeq system, according to manufacturer’s protocol. More specifically, for the library construction, the Illumina indexes were attached to the amplicon with 8 PCR cycles using Nextera Index Primer. The final purified (AMPure XP beads) product was then quantified using qRT-PCR according to the qRT-PCR Quantification Protocol.

Guide (KAPA Library Quantification Kits for Illumina Sequencing platforms) and qualified using the TapeStation D1000 ScreenTape (Agilent Technologies). The paired-end (2 × 300 bp) sequencing was then performed by Macrogen using the MiSeq platform (Illumina).

Sample reads preprocessing

Two samples were acquired for each of the 10 GCs. For each sample, the reads were merged with BBMerge (48) (using the recommended command for optimal accuracy bbmerge-auto.sh). We compared the reads obtained from BBMerge with the ones from PEAR (49) (using a quality index of 30) and no significant differences were found. Then, the primers’ and adaptors’ sequences were removed by using a customized cutadapt script to remove seven primers from the merged reads (50). The quality of the preprocessed reads was tested by the FASTQC tool available in the Galaxy platform (https://usegalaxy.org). The sequences were transformed to FASTA format. Finally, the alignments to V and J germline sequences, and the identification of CDRs/FWRs and frameshifts were obtained by submitting the reads to the IMGT-High-V-Quest web portal (28) (http://www.imgt.org/HighV-QUEST/search.action).

Grouping sequences into clones

Following previous conventions (3, 15, 24, 51, 52, 53), sequences were grouped together into clones if they shared the same V, J gene segments and CDR3 length, and more than 84% junction nucleotide sequence identity. The threshold was optimized using the distance to nearest distribution model (22) and a negation table approach (54), as detailed in the Section 3 in Supplemental Data 1. The junction clustering was performed with a hierarchical agglomerative clustering algorithm (26 Preprint), previously proven to be effective at identifying clones (27). As our data are especially susceptible to incorrect V-gene assignment, due to part of FRW1 missing in our sequencing data (the first 15–20 residues), clones with different V genes from the same subgroup were merged together when the CDR3 of their most common sequence were identical.

As a way of checking for the consistency of our analysis, we alternately grouped sequences into CDR3 groups irrespectively of their V and J genes, which can be directly performed with the TRIP tool (25), a software framework implemented in R shiny.

Similarity between samples

To quantify the similarity between samples X and Y, we first normalized the number of sequences belonging to each clone by the total number of sequences in that sample such that |X| = |Y| = 1. Then, the similarity was defined as the Sørensen–Dice coefficient, known to be more robust than the Jaccard index when some information is missing from the dataset (55) Similarity=2|X∩Y||X|+|Y|=|X∩Y|=1−∑i∈clones|xi−yi|2,(2)where x_i and y_i are the abundance of clone i in sample X and Y, respectively. A similarity of one indicates perfectly identical samples and one of zero that they do not share any clone.

Quantifying clonal diversity

To quantify the diversity of GCs, we used different statistical representations of GC clonal composition: richness, dominance, entropy, and evenness (53). The diversity indices are denoted with the unified framework of Hill (^q D) (56), where values of q that are smaller than unity disproportionately favor rare species (thus more sensitive to singletons), whereas values greater than unity disproportionately favor the most common species.

• • Dominance (1/^∞ D): the clonal dominance is defined as the number of cells belonging to the most abundant clone divided by the total number of observed cells in the GC.

• • Richness (⁰ D): the clonal richness is defined as the number of clones in the GC. To provide an estimate of the total number of clones in a GC based on an incomplete sample, we use the bias-corrected Chao1 formula (57, 58):

NChao=Nobs+f1 (f1−1)2(f2+1),(3)where N_obs is the number of observed clones in the GC, f₁ is the number of clones detected exactly once, and f₂ the number of clones detected exactly twice. In our setup, the second term of Eq. (3) typically equals between 50% and 100% of the number of observed clones.

• • Shannon entropy (log[1 D]): the Shannon entropy has been a popular diversity index in the ecological literature. It is a measure of diversity that accounts “fairly” for the abundance of both rare and frequent clones (59) and is thus less affected by the presence of singletons. Defining the sample size n, and p_i = Xin as the normalized occurrence of clone i in the sample, the Shannon entropy is defined as:

H=−∑i∈clonespiln(pi).(4)

Similar to the Chao richness, we also use an estimator to infer Shannon entropy with incomplete sample information (60). Defining the sample coverage C = 1−f1n, we adjust the relative species abundances with pj = C,HChao=−∑i=1Nobspj ⁡logpj1−1−pjn.(5)

Evenness=exp(HChao)NChao.(6)

The evenness is bounded between 0 and 1, where an evenness of one corresponds to a perfectly homogeneous sample, that is, all clones have the same occurrence. For these calculation, we used the Python library diversity (53) (https://github.com/Aurelien-Pelissier/cdiversity).

Classifying non-functional and functional clones

Each sequence was labeled either as non-functional when having a frameshift or a stop codon or functional otherwise. For each clone, a productivity index, defined as the number of productive sequences divided by the total number of sequences in that clone (thus between 0 and 1), was computed.

We further classified non-functional sequences into three categories: (i) out of frame, (ii) early stop codon induced by VDJ recombination, and (iii) early stop codon induced by SHM. Out of frame sequences were defined as having an out of frame junction, with a frameshift detected by IMGT. In frame non-functional sequences were assumed to arise either from V(D)J recombination when having at least one stop codon in the IMGT N region or from SHM otherwise.

Phylogenetic tree representation

We used phylogenetic trees to visualize the evolution of the antibody receptor sequences. In such representation, each founder cell defines the unmutated germline of a new tree, and newly acquired mutations are ideally represented as downstream nodes. We defined the root by taking the unmutated V, J germline and filling the remaining junction region with the consensus sequence of all available sequences within the clone. To compute the trees, the grouped sequences from each clone were first aligned with ClustalW (61). This alignment was necessary to infer a mutation matrix, as tree inference algorithms typically require sequences to be aligned and of the same length, whereas experimentally determined sequences have different lengths and include insertions and deletions. To infer trees from a large amount of sequences, we used a hierarchical clustering approach, where the most similar sequences were grouped together and progressively aggregated with other groups until the root node (maximum distance) was reached (neighbour joining method in ClustalW). Such a method was preferred over typical Bayesian maximum likelihood estimation (39) or Markov chain Monte Carlo sampling (62) because the latter were intractable for the amount of sequences per clone in our data (>5,000 sequences) and would thus require significant sub-sampling leading to an important loss of information in our phylogeny analysis.

Sequence similarity

The Levenshtein distance (63), defined as the minimum number of edits required to transition from one sequence to the other, is a common metric to quantify sequence similarity. To reduce the bias caused by length differences, we used a normalized Levenshtein distance (64) that incorporates the length of both sequences and satisfies the triangle inequality. Given two strings s₁, s₂ and Lev (s₁, s₂) the Levenshtein distance between s₁ and s₂, the normalized Levenshtein distance Lev_norm is defined as:Levnorm(s1,s2)=2·Lev(s1,s2)|s1|+|s2|+Lev(s1,s2).(7)

Paratope and antibody modeling

The paratope residues were predicted by submitting the antibody sequences to Parapred (65) (https://github.com/eliberis/parapred). To convert the output of Parapred, binding probabilities, into a binary label, we selected a threshold of 0.67, shown to be optimal by the authors of the original article (65). Antibody sequences were submitted to repertoire builder (66) (https://sysimm.org/rep_builder/) to construct the antibody structure from a homology model. As part of the FRW1 region was missing from our sequencing data (the first 15–20 residues), missing nucleotides were filled with the unmutated sequence of their respective IGHV. The obtained structures were visualized with PyMol (67).

Functional convergence

Although antibodies within the same clones are likely to target a common epitope (40), it has also been observed that clones from different V and J gene backgrounds express common epitope reactivity (68). To identify different clones that are potentially binding to the same target epitope, we combine three approaches.

• • CDR similarity: as most of the paratope residues are localized on the CDR loops of the antibody (≈90% according to (69)), antibodies that exhibit high CDR similarity are likely to express common reactivity (although expressing different V or J genes). Given two sequences s₁ and s₂ with CDRs denoted as H1, H2, and H3, we define the CDR similarity as the average Levenshtein distance between their respective CDRs.

CDR similarity(s1,s2)=13∑i=[1,2,3]Levnorm(Hi1,Hi2).(8)

The threshold for common reactivity prediction was set to 84% for this metric, as optimized in Section 6 in Supplemental Data 1.

• • Paratype (40) simplifies the complex phenomenon of antibody–antigen interaction into sets of shared paratope residues. It was shown to help in identifying antibodies from different clones binding to the same target epitope. Two antibodies are said to be paratype if they have a sequence identity across the predicted paratope regions greater than a given threshold. Paratope identity is defined as the number of identical paratope residues (residues that are predicted to be in the paratope in both cases) divided by the smallest number of paratope residues of either sequence being compared. The threshold for common reactivity prediction was set to 76% for this metric, as optimized in Section 6 in Supplemental Data 1.

• • Ab-Ligity (41) is a structure-based similarity measure tailored to the antibody–antigen interfaces. Using predicted paratopes on modeled antibody structures, it allows for the identification of sequence-dissimilar antibodies that target highly similar epitopes. In short, the Ab-Ligity framework enumerates all sets of three residues in the paratope structures, and tokenizes them in terms of both the distance between the residues and their intrinsic chemical properties (aliphatic, hydroxyl, sulphur, aromatic, acidic, amine, basic). The obtained set is then compared among other antibodies to obtain a similarity score, and predicted as common epitope if above 0.26, as optimized in Section 6 in Supplemental Data 1.

To find pairs of antibodies that are likely to bind the same epitope, we combine these three metrics. An antibody pair is predicted to bind to a common epitope if it has a similarity above the given threshold for two of the three metrics. In addition, because a different CDR3 length generally implies a different binding mode, we add the condition that the two antibodies should have the same CDR3 length to be predicted as binding the same epitope. Although some cases of antibodies with different CDR3 length expressing common epitope reactivity were observed (41), such a situation is anecdotal. In fact, removing this constrain considerably reduces the accuracy of the three metrics described above, and existing tools are currently unable to accurately identify such cases, mainly because of the difficulty of modeling the CDR3 loop accurately (70).

Comparison to public V gene database

We compared the V genes frequencies with a public database of 6 human donors in lymph node (24). LN samples were retrieved form tissues labeled as submandibular lymph nodes (MLN) and inguinal lymph nodes.