How many different clonotypes do immune repertoires contain?

TL;DR

This paper reviews estimates of immune receptor diversity, highlights limitations of current methods in capturing small clone populations, and proposes combining statistical and mechanistic models to better estimate total clonotype numbers.

Contribution

It identifies a key limitation in current diversity estimation methods and suggests a novel approach integrating statistical and mechanistic models.

Findings

Current methods underestimate total clonotypes due to ignoring small clones.

The tail of small clones significantly impacts diversity estimates.

Combining models offers a promising strategy for accurate estimation.

Abstract

Immune repertoires rely on diversity of T-cell and B-cell receptors to protect us against foreign threats. The ability to recognize a wide variety of pathogens is linked to the number of different clonotypes expressed by an individual. Out of the estimated $\sim 10^{12}$ different B and T cells in humans, how many of them express distinct receptors? We review current and past estimates for these numbers. We point out a fundamental limitation of current methods, which ignore the tail of small clones in the distribution of clone sizes. We show that this tail strongly affects the total number of clones, but it is impractical to access experimentally. We propose that combining statistical models with mechanistic models of lymphocyte clonal dynamics offers possible new strategies for estimating the number of clones.

Figures

Full text

I Introduction

The diversity of immune repertoires plays an important role in the host’s ability to recognize and control a wide range of pathogens. While actual recognition of an antigen depends on having a relatively specific T-cell or B-cell receptor (TCR or BCR), multiple experimental examples show that reduced receptor diversity may limit the efficacy of adaptive immune repertoires Yager2008 ; Muraro2014 . This effect becomes especially pronounced in individuals infected with cytomegalovirus Khan2002 ; Souquette2017 ; Smithery2018 , or with immunosenescence Qi2014 , when naive clonotypes are significantly reduced, and the organism is left to rely on reduced immune diversity.

The number of different clonotypes, called species richness in ecology, is thus an important quantity to estimate, both biologically and clinically. Here we review some of the experimental and theoretical approaches that have been used to estimate the number of distinct clonotypes in TCR and BCR repertoires, in naive, memory, or unfractionned repertoires. Based on existing repertoire data and computational models, we demonstrate that no statistical method can overcome the limitations of small sampling. We argue that this problem, which is inherent to all existing methods, could be overcome by combining repertoire data with stochastic models of lymphocyte population dynamics, taking into account the caveats of convergent recombination and experimental noise.

II Past estimates

The variable part of each TCR $\beta$ (and BCR heavy) chain is composed by putting together variable (V), diverse (D) and joining (J) regions. TCR $\alpha$ (and BCR light) chains only have V and J genes and no D genes. Additionally each chain experiences additions and deletions of nucleotides at the gene junctions, which increases the diversity. This junctional rearrangements have been identified as the main contributor to sequence diversity Davis1988 ; Murugan2012 . The $\alpha$ and $\beta$ chain generations in TCR are separated in time, and have been shown to be independent Dupic2019a ; Grigaityte2017 . The total number of different $\alpha\beta$ pairs that the generation machinery can produce is much greater than the total number of receptors in the whole human population Mora2016 ; Elhanati2018a ; Dupic2019a . Thus, each person harbors only a small fraction of the potential diversity of receptors. To estimate the number of distinct receptors, we need to count them.

An early quantitative direct estimate of the size of the TCR repertoire dates back to Arstila et al Arstila1999 , following earlier considerations Davis1988 ; Davis1998 . Their approach was to focus on a subclass of receptors (either $\alpha$ or $\beta$) with a specific V-J class and length, and sequence them using low-throughput methods. The number of different sequences in that subclass is then extrapolated back to the get the full diversity by dividing by the known frequency of V-J and length usage. The authors were very careful to verify their estimate in different V-J classes and donors, and to account for rare clonotypes that might have escaped sequencing. A total of $\sim 10^{6}$ different TCR $\beta$ chains were thus estimated in a sample of $10^{8}$ T-cells. Each $\beta$ chain paired with 25 different $\alpha$ chains, resulting in $\sim 25\cdot 10^{6}$ distinct TCR$\alpha\beta$. Much smaller TCR $\alpha\beta$ diversity was reported in the memory subset, $\sim 2\cdot 10^{5}$, consistent with the idea that memory cells form a selected and thus restricted subset.

With the onset of high-throughput sequencing of immune receptor repertoires Weinstein2009 ; Robins2009 ; Boyd2009a ; Six2013 ; Robins2013a ; Georgiou2014 ; Heather2017 ; Minervina2019a ; Bradley2019 came the realization of the importance of the sampling problem. In sequenced repertoires, many clonotypes are seen just once, suggesting that there are possibly many more that have similar or slightly smaller sizes but were not sequenced, simply by chance. This issue does not only affect “small” clones. A clone of $10^{5}$ cells among a total of $3\cdot 10^{11}$ T cells will often not be seen even once in a typical sample of $10^{6}$ cells. To deal with this issue, a commonly adopted approach has been to use statistical estimators (see Laydon2015 for a overview).

Using the Poisson abundance statistical method Fisher1943 , Robins et al Robins2009 obtained estimates of $\sim 10^{6}$ TCR$\beta$ nucleotide clonotypes for CD8 and CD4 naive cells, and $\sim 5\cdot 10^{5}$-$10^{6}$ for CD8 and CD4 memory TCR$\beta$. Qi et al. Qi2014 used another method called the Chao2 estimator Chao2002 , which uses multiple replicates of the sequencing experiment, to estimate TCR species richness. They obtained much larger estimates than Robins et al.: $\sim 2\cdot 10^{8}$ TCR$\beta$ nucleotide clonotypes in CD4 and CD8 naive repertoires, $\sim 1.5\cdot 10^{6}$ clonotypes in CD4 memory repertoires, and about $5$-$10$ times fewer in CD8 memory. All those numbers decreased with age. Using the Poisson abundance method, diversities of $\sim 1$-$2\cdot 10^{9}$ for naive and $\sim 5\cdot 10^{7}$-$10^{8}$ for memory BCR heavy chains were reported DeWitt2016 . Recent estimates of BCR heavy-chain species richness using an advanced statistical estimator Kaplinsky2016 yielded smaller diversities, ranging from $\sim 10^{7}$ Soto2019 to $10^{7}$-$10^{9}$ Briney2019 , presumably because they focused on amino acid rather than nucleotide clonotypes and ignored hypermutations in the V and J segments.

III The sampling problem

The approaches described above share the common problem that it is impossible to extrapolate what happens for small clones from small samples, which capture the largest clones Laydon2015 ; Kaplinsky2016 . Getting information on the small clonse is in fact impractical: it would require sequencing essentially all lymphocytes in an organism. Humans harbor of the order of $3\cdot 10^{11}$ T cells (and roughly the same order of B cells). Of these, only a few percents are contained in blood, of which a small fraction ($\sim 10^{6}$) is sampled in typical experiments (Fig. 1A). Even in mice, which contain fewer lymphocytes ($\sim 10^{8}$ T cells) and can be sacrificed to isolate all the body’s lymphocytes, cell loss during the experiment hampers this approach.

In Fig. 1B-D, we illustrate with simulations what happens when one analyses samples of $10^{6}$ cells from three synthetic repertoires. These repertoires are described by different clone size distributions, corresponding to a widely different number of clones: a pure power law (Fig. 1B), a mixture of a power law and neutral model Hubbell2001 (Fig. 1C), and a power law with a low-frequency cut-off (Fig. 1D). Their species richness are widely different, ranging from $N\sim 7\cdot 10^{5}$ to $1.6\cdot 10^{10}$. Yet, the sampled repertoires show similar clone size distributions, and comparable observed diversity ($10^{5}$-$10^{6}$), because they behave similarly for large clone sizes, but drastically differ in the tail of small clone sizes.

Any statistical method that extrapolates from observations assumes, knowingly or implicitly, an underlying model for how the clone size distribution behaves for the smallest clones. The Poisson abundance and Chao estimators Fisher1943 ; Chao1984 ; Chao2002 discussed earlier both assume a well peaked distribution of clone sizes, which is not the case in our examples. As a result, the Chao1 Chao1984 estimator can underestimate species richness by up to a 1,000-fold factor (Fig. 1B). A more advanced estimator such as Recon Kaplinsky2016 gives even worst underestimates.

Real repertoires are likely affected by this problem. The sampled (large clone) part of their clone size distribution has been shown to follow a power law both for TCR and BCR Mora2016 . Naive subsets display shorter tails of large clones Oakes2017 , suggesting that the power-law behaviour in unfractioned repertoires is dominated by memory clones. Our three synthetic examples are consistent with power laws for large clones, but differ greatly for small clones, yielding very different species richness. Extrapolating the distribution of clone sizes is the key idea behind DivE Laydon2015a and Recon Kaplinsky2016 , which were proposed to estimate diversity in TCR subsets. However, these approaches assume that the behaviour at large clones is informative for small clones, which may not always be true.

IV Proposed solution: stochastic modeling

To access small clones that cannot be directly probed experimentally, we need to explicitly model the biological processes that shape these distributions, without having to take a leap of faith. Unlike extrapolation, such models might predict behaviours for small clones that are quantitatively different than the trend suggested by large clones. Of course, model assumptions should be tested experimentally, their parameters estimated from measurements, and confidence intervals put on their predictions. We now briefly review two simple models that have been proposed to describe the dynamics and clone size distributions of naive and memory repertoires.

Cells in naive repertoires have not experienced strong proliferation due to antigen recognition. Nevertheless not all clones are of the same size, in part because clones leave the thymus with different initial sizes, and in part because they undergo stochastic division and death. The simplest model of naive repertoires is Hubbell’s neutral model of ecology Hubbell2001 , which assumes constant division and death rates ($\nu<\mu$) for each cell, with new clones introduced with rate $\theta$ and constant initial size $k$ (Fig. 2A, left). More complex variants of that model may include intrinsic fitness differences between clones or cells, e.g. through competition for self-antigens Lythe2016 or cytokines Desponds2016 .

Under that simple model, the steady state distribution of clones can be computed analytically Desponds2016 ; Greef2019 ; Altan-Bonnet2019 , and falls off exponentially for clones larger than $k$ (Fig. 2A, right), meaning that large clones are rare. The total number of clones $N$ can also be calculated analytically as a function of the model parameters $\mu$, $\nu$, and $k$, as well as the total number of cells $T$ (Fig 2C). Unless cell division almost exactly balances death ($\nu\sim\mu$), or the introduction clone size $k$ is large, the typical naive clone size is fairly small. This means that the total number of clones is very large, and comparable to the total number of cells. To get a more precise estimate would require to measure the division rate of naive T cells $\nu$, and initial clone size $k$.

A limitation of this approach is the assumption that the clone size distribution quickly reaches a steady state. Naive T cells are very long-lived $\mu^{-1}\sim 3\ {\rm years}$ DeBoer2013 , and the size of the naive pool changes with age, so that steady state may never be reached. Transient models of naive repertoires remain to be explored in more detail. A recent experimental study suggests that some T-cell naive clones are much larger than predicted by the neutral theory Greef2019 . However the origin of these outliers is not yet well understood, and may have to do with the inadequacy of our current definitions of naive and memory cells through surface markers.

Modeling memory repertoires requires taking into account the expansion and then contraction dynamics after an infection. These dynamics are driven by new pathogens that infect the host, are recognized and then cleared DeBoer1995 , which leads to a constantly changing antigenic landscape Desponds2016 . This random encounter with antigens can be simply modeled by bursts of division events for all cells of the same clone, with rate $r$, causing each cell to effectively multiply $m$ times into memory following antigen clearance (Fig. 2B, left). Again, this model can be solved exactly in the continuous limit at steady state. The clone size distribution follows a power law for large clones (Fig. 2B, right). The predicted number of clones $N$ depends critically on the power law exponent $\alpha$, and can be calculated as a function of the model parameters (Fig. 2D). $N$ drops to zero for power law exponent $\alpha$ close to 1. Interestingly, measured exponents $\alpha$ from unfractioned T cell $\beta$ chain repertoires (whose large-clone tail is believed to be dominated by memory clones) range from $1$ to $\sim 1.5$ Mora2016 . This high sensitivity to parameters makes estimates from data very difficult.

Extensions of this model include the emergence of antigenic “niches” Desponds2017 , where clonal expansion is limited by antigen availability, leading to diminishing returns upon multiple stimulation events Mayer2019 . Such mechanisms would limit the size of the largest clones and would cut off the power law behaviour, which is not observed in data.

V Connecting models to data

Several caveats and corrections must be taken into account when linking stochastic models of population dynamics such as discussed above to repertoire data. Of importances are the issues of convergent recombination and experimental noise.

Population dynamics models focus on clones, defined as the set of cells originating from a common recombination event. However, two recombination events can lead to exactly the same sequence. The two corresponding clones would be indistinguishable, and form a single clonotype in the repertoire. This effect can be corrected for by using models of recombination Marcou2018 ; Sethna2019 . These models, which are inferred from data, can predict the distribution of generation probabilities of full receptors, or of single chains, both at the level of nucletoide or amino acid sequences, as shown in Fig. 3A for TCR. From this distribution, the probability of convergent recombination can be computed to predict the number of distinct clonotypes as a function of the number of “clones”, defined as independent recombination events (Eq. 7 of Elhanati2018a ). In Fig. 3B we plot that prediction for the $\alpha$, $\beta$, and $\alpha\beta$ TCR in humans, for both nucleotide and amino acid clonotypes. These computations show that, for the full TCR$\alpha\beta$ clonotypes, convergent recombination is so rare that it hardly affects species richness. For the $\alpha$ and $\beta$ chains alone, however, the effect is substantial. Since most repertoire data are of single chains, this correction should be applied when linking data to the type of models discussed above, as was done in Ref. Greef2019 .

An additionnal issue complicates the comparison of models to data: experimental noise in the observed frequencies of clonotypes. In practice the number of reads (or unique molecular identifiers when they are used) $n$ observed in data for a given clonotype is not simply the result of random sampling, and is not distributed according to a Poisson law, as was assumed in Fig. 1 and in all previous work on diversity estimation. Instead, noise is over-dispersed, due to additional noise caused by DNA amplification and library preparation prior to sequencing. This noise model can be fitted using replicates of the repertoire sequencing experiment. This inference is impossible to separate from the inference of the clone size distribution $\rho(f)$. The two must thus be learned simultaneously from replicates by maximizing the likelihood of observed abundances, which depend on both the clone size distribution and the noise properties Pogorelyy2018c . Applying this approach to the $\beta$ chain of unpartioned T cells with $\rho(f)\propto f^{-1-\alpha}$, yields species richness of $N\sim 10^{8}-10^{9}$, with power-law exponent $\alpha\approx 1$-$1.2$ Pogorelyy2018c . In this estimate, the power law is taken as a given, and not linked to a model of clonal dynamics. A full mechanistic model treatment combined with the statistical model remains an interesting direction to explore, which could help shed light on the differences between memory and naive repertoires.

VI Sampling and repertoire sharing

Several recent papers have focused on shared immune receptors from high throughput BCR repertoire data Greiff2017a ; Briney2019 ; Soto2019 . These high profile analyses report absolute percentages of shared clonotypes. However, it was shown in the context of TCR (but the same holds for BCR) that these fractions are not absolute properties of the repertoires, but rather depend on sampling depth and the number of individuals that share the clonotypes Mora2016 ; Elhanati2018a . Sharing estimates based on samples of the repertoire are bound to grossly underestimate the true sharing fraction. Therefore, reporting sharing percentages without appropriate information about sample sizes is meaningless.

To assess the true overlap between the repertoires of two or more individuals, we would need to sequence all their lymphocytes, which is impractical. However, statistical model of sequence probabilities can be used to extrapolate sharing estimates to the full repertoire size $N$, provided that number is known Elhanati2018a . For instance, clonotypes whose probability $p$ is larger than $1/N$ are expected to be present in $1-e^{-pN}>63\%$ of individuals, and can be considered “public”. Recombination models such as the one of Fig. 3A can be used to estimate the fraction of clonotypes that are public. For example, for $N=10^{10}$, the model predicts that about $15\%$ of TCR$\beta$ amino acid clonotypes expressed by human individuals are public Elhanati2018a .

VII Conclusions

While we have focused our review on the number of distinct clonotypes, what really matters for biological function is the number of different specificities. Due to cross-reactivity, each receptor can recognize many antigens and each antigen can be recognized by many receptors with different strengths. To account for this degeneracy, we would need to define a functional coverage of the antigenic space Zarnitsyna2013 . However, we currently do not have a comprehensive sequence-to-function maps for TCR and BCR that would allow us to estimate such a quantity.

Simply counting clonotypes also ignores their relative abundances. Clonotypes expressed by very few cells may not be as relevant for immune protection as very frequent clonotypes. Other diversity measures such as Hill numbers account for differences in frequencies Yaari2015 ; Mora2016 . Some of these measures are in fact more robust than species richness, because they put more focus on large clones and are less susceptible to sampling noise. Depending on the question, these measures may be better suited than species richness.

Our discussion has focused mostly on T cells, and has ignored the complications of hypermutations in BCR, which cause lineages to split into many clonotypes. Whether diversity is defined at the level of lineages or clonotypes will lead to different answers DeWitt2016 ; Soto2019 ; Briney2019 . Developping specialized population dynamics models of B cell development and affinity maturation that include hypermutations is an interesting research direction.

We emphasized that estimating species richness cannot be disentangled from estimating the full distribution of clone sizes. As we gain insight into various aspects of lymphocyte dynamics, from thymic output to infection and memory formation, better mathematical descriptions can be leveraged to propose refined forms for the clone size distribution, and to fit their parameters to observations. Only with such a combination of modeling and data will we be able to get a better picture of repertoire diversity and immune coverage.

**Acknowledgements. ** This work was partly supported by ERC CoG 724208.