Selection and Neutral Mutations Drive Pervasive Mutability Losses in Long-Lived Anti-HIV B-Cell Lineages.

High-affinity antibodies arise within weeks of infection from the evolution of B-cell receptors under selection to improve antigen recognition. This rapid adaptation is enabled by the distribution of highly mutable "hotspot" motifs in B-cell receptor genes. High mutability in antigen-binding regions (complementarity determining regions [CDRs]) creates variation in binding affinity, whereas low mutability in structurally important regions (framework regions [FRs]) may reduce the frequency of destabilizing mutations. During the response, loss of mutational hotspots and changes in their distribution across CDRs and FRs are predicted to compromise the adaptability of B-cell receptors, yet the contributions of different mechanisms to gains and losses of hotspots remain unclear. We reconstructed changes in anti-HIV B-cell receptor sequences and show that mutability losses were ∼56% more frequent than gains in both CDRs and FRs, with the higher relative mutability of CDRs maintained throughout the response. At least 21% of the total mutability loss was caused by synonymous mutations. However, nonsynonymous substitutions caused most (79%) of the mutability loss in CDRs. Because CDRs also show strong positive selection, this result suggests that selection for mutations that increase binding affinity contributed to loss of mutability in antigen-binding regions. Although recurrent adaptation to evolving viruses could indirectly select for high mutation rates, we found no evidence of indirect selection to increase or retain hotspots. Our results suggest mutability losses are intrinsic to both the neutral and adaptive evolution of B-cell populations and might constrain their adaptation to rapidly evolving pathogens such as HIV and influenza.

Introduction

High-affinity antibodies arise during the adaptive immune response from the very process that gave vertebrates an adaptive immune system in the first place: adaptation by natural selection. In response to infection or vaccination, mutagenic enzymes and error-prone polymerases cause somatic hypermutation of B-cell receptors and thus create variation in their ability to bind antigen (Teng and Papavasiliou 2007). B cells with high-affinity receptors are more likely to receive survival and replication signals from helper T cells, and thus selection for improved antigen binding drives the development of high-affinity B-cell receptors that are later secreted as antibodies (Griffiths et al. 1984; MaClennan and Gray 1986; Liu et al. 1989; Victora and Nussenzweig 2012). Understanding how the immune system evolved to facilitate the rapid adaptation of B-cell receptors during infection may provide general insights into the evolution of adaptability.

Two features of B-cell receptor genes suggest their long-term evolution has been shaped by selection for adaptability during infection. First, the germline genes that recombine to produce B-cell receptors are enriched for nucleotide motifs that are targeted with high frequency by the mutagenic enzymes involved in somatic hypermutation (Rogozin and Kolchanov 1992; Rogozin and Diaz 2004; Wei et al. 2015). High mutation rates provide the genetic variation required for B-cell adaptation, and low B-cell mutation rates have been linked to immunodeficiency disorders (Levy et al. 1998; Bonhomme et al. 2000) and to the decline in immune function with age (Frasca and Blomberg 2009). Second, mutational “hotspots” occur where mutations are most likely to be beneficial. Hotspots are concentrated in loops of the B-cell receptor protein that are directly involved in antigen recognition (complementarity determining regions, CDRs) (Oprea and Kepler 1999). In contrast, structurally important regions of the B-cell receptor (framework regions, FRs), which are usually less directly involved in antigen binding, are enriched with motifs that have low mutability (Oprea and Kepler 1999). In addition, mutations in the CDRs are more likely to be nonsynonymous and involve amino acids with distinct biochemical properties (Wagner et al. 1995; Kepler 1997; Hershberg and Shlomchik 2006; Saini and Hershberg 2015). The differential mutability of CDRs and FRs appears to focus mutations to regions where they are likely to produce variation in antigen affinity without destabilizing the protein. The frequency and distribution of mutational hotspots in B-cell receptor genes therefore seem to contribute to their adaptability during immune responses.

As B cells mutate during the immune response, however, changes in the frequency and distribution of mutational hotspots might affect the subsequent adaptability of B-cell receptors. This change in adaptability may be especially important in B-cell lineages that coevolve with pathogens like HIV and influenza. Experimental removal of hotspots decreases both the mutation rate of the altered site relative to others and the overall mutation rate of the sequence (Wei et al. 2015; Hwang et al. 2017). Loss of highly mutable motifs has been hypothesized to occur during the immune response due to motifs’ propensity to mutate (Hoehn et al. 2017; Stamper and Wilson 2017), and decreased mutation rates due to such “hotspot decay” might explain declines in the evolutionary rates of B-cell lineages over several years of HIV infection (Wu et al. 2015; Sheng et al. 2016). In addition, changes in the distribution of highly mutable motifs across FRs and CDRs might increase the frequency of deleterious mutations in the former and decrease the frequency of beneficial mutations in the latter.

Although hotspots have been shown to spontaneously decay through random mutations (Hoehn et al. 2017), a full picture of the factors influencing the evolution of B-cell receptor mutability is missing. First, highly mutable motifs might be regained through mutation. Second, selection for affinity and protein stability might favor nonsynonymous mutations that incidentally increase or decrease mutability. Finally, selection might act on somatic hypermutation rates themselves. Although selection in theory can favor low mutation rates due to the reduced frequency of deleterious mutations (Kimura 1967; Lynch et al. 2016), rapidly changing environments may indirectly select for a higher mutation rate through its association with beneficial mutations (Leigh 1970; Taddei et al. 1997; Giraud et al. 2001; Pal et al. 2007). Thus, although short-term differences in fitness among B cells arise from differences in the affinity of their receptors, B cells with more mutable CDRs might have a higher probability of producing high-affinity descendants able to keep up with evolving antigens in the long term. Indirect selection for mutability might therefore retain or increase the frequency of highly mutable motifs in CDRs, whereas mutability losses in FRs might be selected due to the reduced frequency of destabilizing mutations.

Understanding changes in mutability may reveal constraints on B-cell adaptation, but the contributions of different mechanisms to changes in the frequency and distribution of mutational hotspots during B-cell responses are largely unknown. We investigated the evolution of B-cell receptor mutability by fitting phylogenetic models to sequences from long-lived anti-HIV B-cell lineages. In characterizing mutability, we considered two sequence-based features that appear to have been strongly selected in the evolution of the adaptive immune system: overall mutability (the density of mutagenic nucleotide motifs) and also changes in the mutability of CDRs relative to FRs. First, we examined B-cell mutability in the unmutated common ancestors of anti-HIV antibodies. Next, we investigated the effects of random mutations and positive selection for amino acid substitutions that increase affinity for antigen. Finally, we tested for selection to increase, retain or decrease the frequency of highly mutable motifs.

Results

Ancestral B Cells Have Higher Mutability in CDRs than FRs

To characterize changes in mutability during B-cell evolution, we inferred the evolutionary histories of previously reported B-cell lineages from three HIV-1 patients. The CH103 and VRC26 lineages comprise heavy and light chain B-cell receptor sequences obtained from high-throughput sequencing over 144 and 206 weeks of infection in two patients, respectively (Liao et al. 2013; Doria-Rose et al. 2014). We also analyzed heavy chain sequences of three lineages from the VRC01 data set, which was sampled from a third patient over a 15-year period (Wu et al. 2015). The lineages we analyzed were originally investigated for having evolved the ability to neutralize diverse HIV strains.

To infer changes in mutability over time, we used Bayesian phylogenetic analyses (see Materials and Methods) to obtain a sample of time-resolved trees from the posterior distribution of each lineage’s genealogy for the heavy and light chains separately, and to estimate the nucleotide sequences of all internal nodes. Mutabilities of the observed and the inferred internal sequences were estimated using the S5F model. This model assigns relative mutation rates to all five-nucleotide DNA motifs and is based on a large independent data set of antigen-experienced B cells (Yaari et al. 2013). The mutability of each sequence was defined as the geometric mean of the S5F scores across all sites in the B-cell receptor sequence. We estimated the number, magnitude, and distribution of mutability changes on all branches by computing the difference in average log-S5F-score for all pairs of parent–descendant nodes. Figure 1 illustrates mutability evolution in the heavy chain of lineage CH103.

. 1.

Evolution of S5F-mutability in the heavy-chain CH103 B-cell lineage. (a) Long-term declines in mutability across the maximum clade-credibility tree. Nodes are colored according to the geometric mean of the S5F mutability scores across their sequences. Branches measured in weeks are annotated to indicate gains (+) or losses (–) of mutability. (b) Mutability over time for a combined sample of 100 trees from the posterior distribution. Blue points correspond to terminal nodes (observed sequences), and black points correspond to inferred internal nodes. The red line represents an average of regression lines calculated for each tree in a sample of 1,000 trees. (c) For each tree in the posterior distribution of trees obtained for each lineage, we computed the magnitude of mutability changes on all branches, producing one distribution per tree. The distributions for a sample of 100 trees are shown in the main plot as overlaid densities. We also computed, for each tree, the fraction of changes in that tree that were losses, producing one value per tree. A distribution of such values for the full sample of 1,000 trees is shown in the inset plot, with the 95% highest posterior density interval shown in gray. Red lines indicate the means of the distributions.

To investigate the potential for mutability changes during the response, we first characterized the mutability of each lineage’s ancestor. In the ancestors of all heavy and light chains, sites in CDRs had higher average mutability than sites in FRs. On an average, across the seven heavy and light chain ancestors, mutability was ∼35% higher in CDRs than in FRs (range: 26–54%; supplementary fig. S1, Supplementary Material online). Previous analyses of germline V genes (which, together with D and J genes, recombine to produce mature B-cell receptors) showed that mutability is lower in FRs and higher in CDRs than expected based on their amino acid sequences, suggesting selection to increase the frequency of highly mutable motifs in CDRs and decrease their frequency in FRs (Oprea and Kepler 1999; Saini and Hershberg 2015). Consistent with those analyses, we found that the FRs of lineage ancestors had lower S5F mutability than expected based on their amino acid sequences. On an average, across all heavy and light chains, the mean FR mutability of ancestral B-cell receptors was lower than 99% of sequences obtained by randomizing their codons (according to usage frequencies in humans; Nakamura et al. 2000) while keeping the amino acid sequences constant (range: 96–100%; supplementary fig. S1, Supplementary Material online). However, while previous studies found that the CDRs of V genes are typically more mutable than expected based on their amino acid sequences (Oprea and Kepler 1999; Saini and Hershberg 2015), we found the S5F mutability of CDRs in ancestral B-cell receptors to be, on an average, greater than only 60% of randomized sequences with the same amino acid sequence (range: 22–92%, supplementary fig. S1, Supplementary Material online). This difference in the results for CDRs may be due to the different mutability metrics used in those studies. Both studies defined mutability as the propensity for nonsynonymous mutations. In contrast, we used the S5F model to quantify the propensity for all mutations. These results suggest that loss of S5F mutability in FRs is only possible if those regions undergo amino acid changes.

Mutability Is More Often Lost than Gained

Next we investigated long-term changes in mutability as B-cell receptors evolved from their ancestors. Consistent with previous analyses (Sheng et al. 2016; Hoehn et al. 2017), average mutability decreased with time in four of the seven heavy and light chains (fig. 1 and supplementary fig. S2, Supplementary Material online). To investigate the factors contributing to those net long-term trends, we analyzed the frequency of mutability gains and losses across branches. Mutability losses should arise from hotspot decay, positive selection of the amino acid changes that incidentally decrease mutability, selection for lower mutability due to the reduction in the frequency of deleterious mutations, or a combination of those factors. Mutability gains should reflect spontaneous hotspot gains through mutation, positive selection of amino acid changes that incidentally increase mutability, indirect selection for higher mutation rates by association with beneficial mutations, or a combination of those factors. To summarize the net contributions of mutability-decreasing and mutability-increasing mechanisms, we computed the fraction of branches with mutability losses out of all branches with mutability changes. By computing the frequency of mutability losses for each tree in the posterior sample, we estimated the posterior distribution (fig. 1).

Across the entire BCR sequence, mutability losses occurred more frequently than gains in four of the seven heavy and light chains (figs. 1 and 2). On an average, across the seven data sets, ∼61% of changes in mutability were losses (range: 46–72%). The four chains where mutability losses were significantly more frequent than gains include three chains where mutability decreased in the long term (heavy and light chains of CH103 and light chain of VRC26, supplementary fig. S2, Supplementary Material online) and one chain where mutability increased in the long term despite the high frequency of mutability losses (VRC01-13). This last result illustrates how mutability can increase in the long term despite the propensity for mutability to be lost on any given branch, potentially reflecting the effects of genetic drift, selection for mutations that incidentally increase mutability, or selection for mutability itself.

. 2.

Frequency of losses relative to the total number of changes in mean log-S5F mutability during the evolution of anti-HIV B-cell lineages. Results are shown for the entire analyzed region of the B-cell receptor, and separately for framework regions (FRs) and complementarity determining regions (CDRs). Each point denotes the fraction of changes in mean log-S5F mutability that were losses, averaged across a sample of 1000 trees from the posterior distribution. Vertical red lines indicate the 95% highest posterior density interval.

The frequency of mutability losses was not significantly different from the frequency of gains in the heavy chains of VRC26 (95% highest posterior interval 39–73%) and VRC01-19 (44–79%) and was slightly lower than the frequency of gains in lineage VRC01-01 (43–50%).

Mutability changes in FRs and CDRs were consistent with the changes observed across the entire B-cell receptor sequence. Long-term declines in mutability occurred in the FRs of four of the seven heavy and light chains (supplementary fig. S3, Supplementary Material online) and in the CDRs of five of the seven heavy and light chains (supplementary fig. S4, Supplementary Material online), and the difference between CDR and FR mutabilities changed little (supplementary fig. S5, Supplementary Material online). Consistent with the net long-term trends, both CDRs and FRs had similar frequencies of mutability losses (FR average 56%, range: 39–70%; CDR average 54%, range: 42–63%; fig. 2). Despite considerable amino acid divergence from the unmutated ancestors in FRs (16–67% for different heavy and light chains), sequences from the last sampling time in each data set had lower FR mutability than 94% of randomizations with the same amino acid sequence (range: 61–99.9%, supplementary fig. S6, Supplementary Material online). The lower mutability of ancestral FRs relative to their amino acid sequences was therefore retained throughout the evolution of the B-cell lineages. CDRs, however, became less mutable relative to their amino acid sequences than the ancestors. On an average, across the seven heavy and light chains, CDRs of sequences from the last sampling times were more mutable than 35% of their corresponding randomizations (range: 7–64%, supplementary fig. S6, Supplementary Material online), down from 60% in ancestral CDRs.

Hotspot Decay and Selection for Amino Acid Substitutions Contribute to Mutability Losses

Highly mutable motifs have been hypothesized to decay due to their propensity to mutate (Hoehn et al. 2017; Stamper and Wilson 2017), but motifs should also be influenced by positive selection. The first case is straightforward. For instance, the hypothetical sequence CAGCTT contains the highly mutable cytosine at the center of the AGCTT motif (Rogozin and Diaz 2004; Yaari et al. 2013; Wei et al. 2015), and a C →T mutation in the underlined position would decrease the mutability of the site ∼5-fold (Yaari et al. 2013). Positive selection on amino acid substitutions should influence motifs in two ways. Mutational hotspots can be disrupted if selection favors amino acid sequences whose codons happen to be, on an average, less mutable than codons in the ancestral sequence. Selection can also affect mutability in neighboring codons. For example, if selection led to the replacement of CAG (glutamine) for CGG (arginine) in the sequence CAGCTT, the mutability of the underlined C nucleotide (not involved in the substitution) would decrease ∼13-fold (Yaari et al. 2013).

To test if the observed losses of S5F mutability were consistent with the decay of highly mutable motifs expected under the S5F model, we simulated B-cell receptor evolution under a model that allows for variation in mutation rates across motifs based on their S5F mutability scores (Yaari et al. 2013) (see Materials and Methods). We compared changes in mutability simulated under the S5F-based model to changes under models that do not allow for motif-driven mutation rate variation. Instead, these models assume that the mutation rate is identical across all sites (“uniform” model) or depends on the position within a codon (“codon-position” model), with the relative rate for each position estimated from the data (see Materials and Methods). For each branch on the maximum clade credibility (MCC) tree of each lineage, we started from the inferred nucleotide sequence of the branch’s parent node and simulated 100 replicates of a descendant sequence. We constrained the simulations to produce the same number of synonymous and nonsynonymous changes as inferred from the data but allowed them to occur at different positions. Thus, the overall amount of selection is expected to be the same in the data and in the simulations, and the models differ in the locations of synonymous and nonsynonymous mutations. We then partitioned observed and simulated changes in mean log-S5F mutability into changes caused by synonymous mutations and changes caused by nonsynonymous mutations.

Observed synonymous changes in mutability were consistent with the decay of mutational hotspots simulated under the S5F model (fig. 3). In contrast, synonymous changes simulated with constant mutation rates across motifs (uniform and codon-position models) increased mutability in all lineages. Those results suggest that hotspot decay explains most of the mutability loss caused by neutral mutations, and that the S5F model, in particular, accurately describes neutral sequence evolution in B-cell lineages.

. 3.

Changes in mean log-S5F mutability due to synonymous and nonsynonymous changes in anti-HIV B-cell lineages. (a) Changes in mutability due to synonymous changes, averaged across all branches. Gray points indicate values inferred from the data, and colored points indicate values obtained by simulation under different models. Red indicates an S5F-based model where different nucleotide motifs mutate with different rates, dark blue indicates a model with no mutation rate variation across sites, and light blue indicates a model with different mutation rates for each position of a codon. Simulations were performed independently for each branch on the MCC tree of different anti-HIV B-cell lineages, starting from the inferred sequence of the parent node. Each simulated sequence was constrained to have the same number of nonsynonymous and synonymous changes as observed in the branch. Vertical bars indicate the 95% range obtained from 100 simulations per model. (b) Changes in mean log-S5F mutability due to synonymous and nonsynonymous changes summed across all branches of the B-cell genealogies. (c) Changes in mean-log S5F mutability in the CDRs of 13 ancestor–descendant B-cell receptor sequence pairs where binding affinity to HIV-1 increased (Liao et al. 2013).

However, the average nonsynonymous loss simulated under the S5F model was greater than the observed loss in 13 of 14 regions/lineages (except in the CDRs of light-chain VRC26, supplementary fig. S7, Supplementary Material online). Although those differences are within the 95% range of values simulated under the S5F model, the consistent overestimation of nonsynonymous mutability losses by the S5F model suggests the distribution of nonsynonymous substitutions is affected by mechanisms other than across-motif variation in mutation rates, such as positive or purifying selection on specific amino acid mutations. Purifying selection may counter mutability losses if it eliminates nonsynonymous changes that happen to decrease mutability. As previously reported by Sheng et al. (2017), we found that several nonsynonymous changes occurred with high probability under the S5F mutability model but were scarce in the MCC trees inferred from the data (supplementary fig. S8, Supplementary Material online), which might indicate they were under purifying selection during the evolution of the B-cell lineages. Mutability losses caused by the ten most common transitions simulated under the S5F model were on an average 75% greater than the losses caused by the ten most common transitions observed in the data (which in fact caused mutability gains on an average in two data sets). Purifying selection against some of the most likely amino acid transitions under the S5F model would therefore contribute to smaller nonsynonymous mutability losses than expected.

To evaluate the strength of positive selection in the B-cell lineages, we used BASELINe (Yaari et al. 2012) to quantify deviations in the frequency of nonsynonymous substitutions from its expected value in the absence of selection (and under the mutational biases captured by the S5F model). In line with previous analyses of the same data sets (Sheng et al. 2016) and of B-cell lineages from healthy donors (Yaari et al. 2015), we detected an enrichment of nonsynonymous changes in the CDRs of four of the seven heavy and light chains (lineages CH103 and VRC26) and a lower frequency of nonsynonymous mutations than expected in the FRs of all heavy and light chains (fig. 4). This result contrasts with whole-repertoire analyses showing predominantly purifying selection across both types of regions (McCoy et al. 2015), and suggests positive selection predominates in the CDRs. However, interpretation of such dN/dS ratios as selection strengths is complicated by the fact that most nonsynonymous changes observed in the B-cell sequences are polymorphisms, not fixations.

. 4.

Selection in framework regions (FRs) and complementarity determining regions (CDRs) of B-cell receptors from anti-HIV B-cell lineages. Selection strength is measured as the log-odds ratio between the observed ratio of nonsynonymous to synonymous substitution frequencies, , and the ratio expected under the S5F mutability model in the absence of selection, (Yaari et al. 2012, 2013). Selection strength values >0 indicate positive selection, and values <0 indicate purifying selection.

The enrichment of nonsynonymous mutations in CDRs suggests that positive selection is concentrated in those regions. To estimate the potential contribution of positive selection to the mutability loss in the CDRs, we summed nonsynonymous changes in the mean log-S5F mutability of CDRs across all branches of the MCC trees. Nonsynonymous changes caused a net loss of CDR mutability in the same four heavy and light chains where an enrichment of nonsynonymous changes was detected (fig. 3). Similarly, synonymous substitutions caused a net loss of CDR mutability in the same lineages. On an average across those four lineages, nonsynonymous substitutions accounted for ∼79% of the inferred mutability loss in CDRs (range: 54–91%). On an average, selection for amino acid substitutions might therefore contribute to as much as 79% of the loss of mutability in the CDRs (in the extreme case where all nonsynonymous changes in the CDRs were under positive selection).

To investigate changes in mutability due to nonsynonymous changes for which there is evidence of positive selection, we analyzed 13 pairs of ancestor–descendant sequences experimentally characterized by Liao et al. (2013). In all of those pairs, binding affinity increased to at least one of the two tested HIV-1 envelope proteins. In ten of the pairs, the increase in affinity was associated with one to ten nonsynonymous changes in the CDRs. At least some of those nonsynonymous changes were thus likely under positive selection for their effect on affinity. Our analysis showed that, in eight of the ten pairs where CDRs underwent nonsynonymous changes, those nonsynonymous changes caused CDR mutability to decrease (fig. 3). Positive selection of amino acid mutations that incidentally disrupted highly mutable motifs therefore likely contributed to the observed mutability loss in the CDRs.

No Evidence of Selection on Mutability Itself

Under persistent or recurrent selection, alleles that increase the mutation rate may increase in frequency by hitchhiking with beneficial mutations, which themselves arise more frequently due to the increased mutation rate (Leigh 1970; Taddei et al. 1997; Giraud et al. 2001; Pal et al. 2007). Such indirect selection for increased mutation rates might therefore lead to the conservation of highly mutable motifs in CDRs, where mutations that improve affinity are selected for during B-cell evolution. In contrast, selection to reduce the frequency of deleterious mutations (Kimura 1967; Lynch et al. 2016) might directly favor mutability losses in FRs.

To test if mutability is subject to direct or indirect selection over the long term, we compared the frequency of synonymous changes in mutability between terminal and internal branches of the B-cell trees. Terminal branches are expected to be enriched for deleterious mutations (Williamson and Orive 2002; Edwards et al. 2006; Kosakovsky Pond et al. 2006; Lemey et al. 2007). If mutability losses are deleterious in the long term, branches where mutability losses occurred should be less likely to contribute descendants to the B-cell population, and would therefore be more likely terminal, rather than internal branches. In contrast, mutability losses should be more frequent on internal branches if mutability losses are beneficial. Because changes in mutability due to nonsynonymous substitutions may be driven by selection for B-cell receptor affinity and stability, we first restricted the analysis to mutability changes from synonymous substitutions.

We found no consistent evidence of selection to increase, retain, or decrease mutability (fig. 5). In one instance (the CDRs of light chain VRC26), mutability losses were more frequent on terminal branches than on internal branches, possibly indicating selection against mutability losses. In another instance (FRs of heavy chains VRC26), mutability losses were more frequent on internal branches than on terminal branches, suggesting selection for mutability losses. In the remaining cases, the frequencies of synonymous mutability losses were similar in terminal and internal branches for both FRs and CDRs. No general trends were seen when we compared terminal branches, internal branches that belong to the main “backbone” of the tree, and other internal branches (Lemey et al. 2007) (supplementary fig. S9, Supplementary Material online), or when we included mutability changes due to nonsynonymous substitutions (results not shown).

. 5.

Frequency of losses relative to the total number of changes in mean log-S5F mutability caused by synonymous substitutions during the evolution of anti-HIV B-cell lineages. Blue indicates changes on terminal branches, and black indicates changes on internal branches. Results are shown separately for framework regions (FRs) and complementarity determining regions (CDRs). Each point denotes the frequency of changes in mutability that were losses, averaged across a sample of 1000 trees from the posterior distribution. Vertical lines indicate the 95% highest posterior density interval.

Results Are Consistent for Three of Four Mutability Metrics

In addition to the S5F model (Yaari et al. 2013), we repeated the analyses using three other mutability metrics. The different mutability metrics were estimated from databases of somatic mutations not subject to the effect of selection, such as synonymous mutations, mutations in noncoding flanking regions of V genes, and mutations in unproductive B-cell receptor genes that are not expressed by B cells but still undergo somatic hypermutation. Two of the alternative metrics are based on a discrete classification of DNA motifs into either hotspots or regular motifs: the number of WRCH and DGYW hotspots (Rogozin and Diaz 2004) or the number of “overlapping” (WGCW) hotspots (Wei et al. 2015), where and . We also quantified mutability using the 7-mer model, which, similar to the S5F model, assigns relative mutation rates to different motifs on a continuous scale, but does so for seven-instead of five-nucleotide motifs (Elhanati et al. 2015).

Consistent with the results for S5F mutability, hotspots and overlapping hotspots decreased in number over time (supplementary figs. S10 and S11, Supplementary Material online) and were lost more frequently than gained (supplementary fig. S13, Supplementary Material online). In contrast, average 7-mer mutability increased over time in four of five lineages (supplementary fig. S12, Supplementary Material online), and losses of 7-mer mutability were approximately as frequent as gains across the entire B-cell receptor sequence (supplementary fig. S13, Supplementary Material online).

Discussion

We used phylogenetic methods to reconstruct changes in B-cell receptor sequences and found consistent loss of mutational hotspots over the course of the adaptive immune response. Our analyses shed light on previous studies of long-term trends in mutability (Sheng et al. 2016; Hoehn et al. 2017) by quantifying the contributions of different mechanisms to these losses. Selection for amino acid substitutions appears to have contributed to the mutability loss in the CDRs—precisely the regions where high mutation rates contribute the most to adaptation. However, mutability changes caused by synonymous substitutions throughout the sequence demonstrate mutability loss from the spontaneous decay of highly mutable motifs through neutral mutations. Disruption of highly mutable motifs through selection and hotspot decay therefore seems to be intrinsic to the evolution of B cells during affinity maturation. These factors counteract the high mutation rate selected in the evolution of immunoglobulin genes, and they suggest adaptability may inevitably become compromised in evolving B-cell lineages.

Observed mutability losses due to nonsynonymous changes were on an average smaller than expected under the S5F model. Purifying selection against amino acid changes that occur with high probability under the S5F model may have contributed to this result. For instance, we found tyrosine-to-cysteine to be one of the most common amino acid mutations under the S5F model, and that mutation typically resulted in a loss of mutability in simulations (supplementary fig. S8, Supplementary Material online). Yet Y-C mutations occurred with much lower frequency in the data, suggesting that purifying selection against Y-C mutations reduced the nonsynonymous mutability loss relative to the S5F expectation.

Although mutability losses occurred both in FRs and CDRs, the higher mutability of CDRs relative to FRs observed in ancestral receptors persisted over at least several years of B-cell evolution, suggesting that some degree of adaptability may still be maintained. This observation may reflect survival bias: lineages in which mutability differences between CDRs and FRs decrease over time may be less likely to persist in the long term.

Theoretical simulations (Taddei et al. 1997) and experimental evolution of bacteria (Giraud et al. 2001; Pal et al. 2007) suggest recurrent adaptation in asexual populations can select for increased mutation rates over the long term, but we found no consistent evidence of selection to increase or retain mutability during B-cell evolution. Any long-term benefits of high mutation rates (in terms of increased chance of producing beneficial mutations) appear insufficient to overcome the rapid decay of highly mutable motifs through mutation and positive selection on amino acid substitutions. We also did not detect selection for reduced mutability in FRs, despite the potential for mutability losses to decrease the rate of destabilizing mutations. Because mutability is already low in the FRs of V genes (Oprea and Kepler 1999; Saini and Hershberg 2015) and ancestral B-cell receptors (this study), it is possible that the fitness effect of mutability losses in FRs is so small that such losses are nearly neutral given the effective size of the B-cell population. Under this “drift-barrier” hypothesis (Lynch et al. 2016), selection for decreased mutability would therefore not be able to fix mutability losses.

Losses of mutational hotspots in CDRs may reduce the adaptability of experienced B cells since lower mutation rates reduce the supply of genetic variation available for selection. Declines in mutation rates caused by hotspot losses may have contributed to reported declines in the evolutionary rates of broadly neutralizing B-cell lineages during chronic HIV infection (Sheng et al. 2016). However, these losses may not be important for lineages that bind to conserved sites. As B-cell receptors evolve high affinity for conserved sites, beneficial substitutions become rare, and the corresponding transition from positive to purifying selection should cause substitution rates to fall (Sheng et al. 2016). In line with the previous analysis by Sheng et al. (2016), we found the oldest lineages to be under the strongest purifying selection in the CDRs. Those lineages also had minimal cumulative changes in mutability over the sampling period, suggesting most of the mutability loss occurred during early adaptation. However, as a result of hotspot loss, these lineages might adapt poorly if formerly conserved antigenic sites suddenly acquire mutations.

A direct link between motif-based mutability scores and rates of B-cell evolution is still lacking. Hotspot loss is associated with reduced substitution rates (Wu et al. 2015; Sheng et al. 2016), which are influenced not only by the underlying mutation rate but also by changes in selective pressures and generation times (Sheng et al. 2016). Using robust counting (O’Brien et al. 2009) and a random local clock model (Drummond and Suchard 2010), we found no evidence for consistent declines in the synonymous or nonsynonymous substitution rates of these lineages (supplementary figs. S14–S17, Supplementary Material online). Simulation-based power analysis (Supplementary Methods) revealed that synonymous and nonsynonymous substitution rates had to decline at least 50% to be detected with this method (supplementary figs. S18–S20, Supplementary Material online). Observed long-term declines in S5F mutability were <50% (supplementary fig. S2, Supplementary Material online) and suggest their effects on substitution rates would not have been detected, assuming a one-to-one relationship between the geometric mean of S5F scores and the actual mutation rate. Additionally, inference of B-cell genealogies may be affected by nonequilibrium nucleotide frequencies, distorting inferences of absolute rates (Hoehn et al. 2017). Standard phylogenetic methods such as the ones used in this study further assume that sites evolve independently, an assumption violated by the context-dependency of mutations in B-cell receptor sequences. More elaborate substitution models that account for the particularities of B-cell evolution (McCoy et al. 2015; Cui et al. 2016; Hoehn et al. 2017) might overcome some of these limitations and help quantify the relationships between mutability metrics and absolute mutation and substitution rates.

Our results suggest a trade-off between the short- and long-term adaptability within B-cell lineages. Repeated infections by antigenically related pathogens, such as influenza viruses, often recall memory B-cell populations (Wrammert et al. 2011; Li et al. 2012; Jiang et al. 2013; Andrews et al. 2015). Preferential recruitment of experienced, less mutable B cells at the expense of naive B cells may compromise the long-term adaptability of the immune repertoire to these pathogens. Protection against such pathogens may rely on a robust naive response that can complement preexisting B-cell lineages once they fail to adapt to new antigens. Weaker naive responses, for instance in the elderly (Frasca and Blomberg 2009), may thus be problematic. Finally, some strategies for eliciting universal responses against HIV and influenza attempt to recapitulate the evolution of long-lived, highly mutated B-cell lineages that were found to produce broadly neutralizing antibodies in infected patients (Haynes et al. 2012). Our results suggest those approaches may be hindered by a decrease in the adaptability of highly mutated antibodies.

Materials and Methods

Sequence Data

We analyzed BCR sequences from three published studies of B-cell lineages sampled longitudinally in individual HIV-1 patients not subject to antiretroviral therapy. The CH103 lineage comprises broadly neutralizing antibodies (bnAbs) isolated from individual B cells and clonally related sequences bioinformatically isolated from high-throughput sequencing over 144 weeks of HIV-1 infection in a single donor (Liao et al. 2013). Likewise, the VRC26 lineage comprises both bnAbs isolated from individually sorted B cells and clonally related heavy-chain sequences obtained from high-throughput sequencing over 206 weeks of HIV-1 infection in a different patient (Doria-Rose et al. 2014). In both cases, BCR sequences obtained from high-throughput sequencing were classified as clonally related to the isolated antibodies based on V and J gene usage. A third lineage, VRC01, was sampled longitudinally from a third HIV patient for 15 years starting from ∼5 years after the date of infection. VRC01 is a large lineage that possibly consists of multiple independent B-cell lineages (Wu et al. 2015; Sheng et al. 2016). We therefore used PARTIS, a recently developed hidden-Markov-model-based method (Ralph and Matsen 2016) to partition the heavy-chain VRC01 data set into sets of sequences likely to have descended from the same naive B cell, and to identify their most likely germline V and J genes. We analyzed the three largest lineages identified in this way, hereafter VRC01-13, VRC01-01, and VRC01-19. We aligned each set of sequences in MACSE v1.01 b (Ranwez et al. 2011) along with their concatenated V and J genes (for the heavy chains) or along with their V genes only (light chains) (table 1). Including the J genes in the light chain data sets produced bad alignments across the J region.

Table 1.

BCR Alignments Analyzed.

Lineage (chain)	V Gene	J Gene	N Seqs.	N Sites	Sampling Times (w.p.i.)	Seqs. per Point
CH103 (heavy)	IGHV4-59*01	IGHJ4*01	460	321	14, 53, 92, 66, 136, 140, and 144	2–232
CH103 (light)	IGLV3-1*01	–	175	285	14, 53, 92, 136, and 144	3–83
VRC26 (heavy)	IGHV3-30*18	IGHJ3*01	681	489	38, 48, 59, 119, 176, and 206	4–273
VRC26 (light)	IGLV-51*02	–	464	294	38, 48, 59, 119, 176, and 206	6–212
VRC01-13 (heavy)	IGHV1-2*04	IGHJ1*01	157	351	240, 544, 580, 784, 808, 832, 856, 884, 924, and 948	2–43
VRC01-01 (heavy)	IGHV-ORF15-1*04	IGHJ4*03	124	372		0–35
VRC01-19 (heavy)	IGHV1-2*04	IGHJ1*01	110	438		1–28

Note.—The number of sequences includes the V + J or V germline sequences. w.p.i., weeks postinfection.

Phylogenetic Inference

Using BEAST v.1.8.2 (Drummond et al. 2012), we fit Bayesian phylogenetic models to the BCR sequence data in order to estimate time-resolved genealogies and internal node sequences for heavy and light chains separately. We used a GTR nucleotide substitution model, assumed a random local clock model to account for potential variation in substitution rates (Drummond and Suchard 2010), and enabled robust counting (O’Brien et al. 2009) to estimate the numbers of synonymous and nonsynonymous changes on each branch (supplementary table S1, Supplementary Material online). To reduce the number of parameters, we used empirical base frequencies and assumed a shared nucleotide transition matrix across codon positions for the robust counting inference. The inferred dynamics of mutability losses and gains were qualitatively robust to the choice of demographic model (constant population size, logistic, or exponential growth) used to calculate coalescent prior probabilities. We set the V + J or V germline sequence of each alignment as the outgroup and assigned them a sampling time of zero to represent the assumption that, at the start of the infection, the ancestral BCR sequence was close to its corresponding germline genes (except for insertions and deletions at the junctions). For each data set, four independent MCMC chains were set to run for 500 million steps and sampled every 1,000 steps (for parameter values) and every 10,000 steps (for trees). We downsampled the set of trees recovered by the MCMC chains to obtain 1,000 trees per chain, sampled at regular intervals between the end of the the burn-in and the end of the entire chain.

Because of the long computation times for the larger data sets (lineages VRC26 and CH103), their MCMC chains were interrupted before 500 million steps had been reached. With the exception of heavy-chain VRC26, interrupted chains had ESSs close to or greater than 200 for most parameters and for the likelihood, prior and posterior probabilities. Although MCMC chains failed to converge for the heavy chain data of lineage VRC26, estimates of the parameters of interest (mean mutability changes and fraction of mutability losses across the tree) were numerically close across replicate MCMC chains, and we therefore present the results in the main text.

To investigate their ability to detect consistent changes in mutability, we repeated the phylogenetic analyses on alignments simulated using different population genetic models (Supplementary Methods). For alignments simulated under a model where all sites are equally likely to mutate, mutability losses were estimated to be approximately as frequent as gains. As expected, mutability losses were estimated to be more frequent than gains in data sets simulated under an S5F-based model, where motifs with high S5F score have a high probability of mutating.

Quantifying Mutability

We quantified the mutability of different sites based on the S5F mutability model by Yaari et al. (2013), which assigns relative mutation probabilities to each of the 1,024 five-nucleotide DNA motifs. Briefly, the scores were derived by counting the number of synonymous mutations associated with each motif in a curated database of mutations from high-throughput sequencing studies, followed by normalization of mutational counts based on the frequencies with which the motifs occurred in the data and statistical estimation of the scores for absent or underrepresented motifs. The use of synonymous mutations presumably eliminates the influence of selection and thus captures the intrinsic propensity of the different motifs to mutate.

Mutability of Randomized Sequences

To estimate the expected mutability of B-cell receptors given their amino acid sequence under random codon usage, we randomized each sequence by sampling, for each amino acid, a codon from the set of codons encoding that same amino acid. The probability of sampling each possible codon was equal to its relative usage frequency in humans (Nakamura et al. 2000). We computed the geometric mean of S5F mutability scores for 1,000 randomizations of each sequence.

Simulations of Sequence Evolution on MCC Trees

For each branch on an MCC tree, we compared the nucleotide sequences of the parent and descendant nodes to identify the number of codon sites with nonsynonymous and synonymous differences. Starting from the parent sequence, we simulated an alternative descendant sequence constrained to have the same number of nonsynonymous and synonymous differences from the ancestor. To introduce each mutation, we sampled a nucleotide site according to different models for variation in mutation rates across sites (see below). We then mutated the nucleotide at the sampled site, and either kept the mutation if the number of mutations of the same type (synonymous or nonsynonymous) was less than the corresponding number inferred from the data, and otherwise rejected it.

To simulate sequence evolution under an S5F-based model, we sampled mutated sites with sampling probabilities proportional to the S5F-scores of the five-nucleotide motifs centered at that each site, as in Yaari et al. (2013). In addition to the S5F-based model, we performed simulations under two models that assume mutation rates do not depend on the local nucleotide sequence around each site. First, we simulated sequence evolution under a uniform model where all sites are equally likely to mutate. Second, we used a codon-position-based model, where codon positions 1, 2, and 3 have different mutation rates. Estimates of the relative mutation rates of each codon position were obtained from the robust counting inference performed in BEAST and used to parameterize the simulations.

To model nucleotide transitions, we made use of the fact that, in addition to estimating the relative mutation rate of different motifs, the S5F model includes the probability of transitions between nucleotides, given the motif where a mutation occurs. All simulations used the S5F-based nucleotide transition probabilities from (Yaari et al. 2013), regardless of whether they assumed variable or constant mutation rates across motifs.

Quantifying dN/dS

Because of variation in relative mutation rates across sites, standard dN/dS tests are unreliable when applied to B-cell receptor sequence data. For example, an excess of nonsynonymous mutations relative to synonymous mutations may result from high mutability in sites where mutations tend to be nonsynonymous, leading to incorrect inference of positive selection. To estimate the strength of selection in B-cell receptor sequences, we followed (Sheng et al. 2016) and used BASELINe (Yaari et al. 2012) to compare the observed ratio of synonymous and nonsynonymous mutations to the ratio expected under motif-specific variation in mutation rates. Briefly, selection strength is quantified as the log-odds ratio between and , where π and are the observed and expected frequencies of nonsynonymous mutations, respectively. The numbers of synonymous and nonsynonymous mutations for each sequence are obtained by comparing the sequence to a reference sequence, which, in our case, was the ancestral sequence of the corresponding lineage.

Click here for additional data file.