Dynamics of B cell repertoires and emergence of cross-reactive responses in patients with different severities of COVID-19

Abstract

Individuals with the 2019 coronavirus disease (COVID-19) show varying severity of the disease, ranging from asymptomatic to requiring intensive care. Although monoclonal antibodies specific to the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) have been identified, we still lack an understanding of the overall landscape of B cell receptor (BCR) repertoires in individuals with COVID-19. We use high-throughput sequencing of bulk and plasma B cells collected at multiple time points during infection to characterize signatures of the B cell response to SARS-CoV-2 in 19 individuals. Using principled statistical approaches, we associate differential features of BCRs with different disease severity. We identify 38 significantly expanded clonal lineages shared among individuals as candidates for responses specific to SARS-CoV-2. Using single-cell sequencing, we verify the reactivity of BCRs shared among individuals to SARS-CoV-2 epitopes. Moreover, we identify the natural emergence of a BCR with cross-reactivity to SARS-CoV-1 and SARS-CoV-2 in some individuals. Our results provide insights important for development of rational therapies and vaccines against COVID-19.

Keywords: B cell clonal expansion; B cell repertoires; BCR selection; BCR sharing; COVID-19; SARS-CoV-2; antibody; cross-reactivity.

Graphical abstract

Figure 1

Roadmap for analysis of BCR repertoires Top: we collected bulk blood IgG BCR samples from 3 healthy individuals and 2 individuals with mild, 12 with moderate, and 5 with severe symptoms of COVID-19 (different markers and colors). We also collected CD38⁺ plasma B cells from PBMC samples of 7 individuals in this cohort (6 moderate and 1 severe) and from 7 additional individuals (2 asymptomatic, 3 mild, and 2 moderate) and 3 healthy individuals (Data S1). Samples were collected at different time points during infection (shown in the center for bulk repertoires). We distinguished between productive and unproductive receptors that had frameshifts because of V(D)J recombination. Line segments of varying lengths represent full V(D)J rearrangements (colors). For each individual, we constructed clonal lineages for productive and unproductive BCRs and inferred the naive progenitors of lineages (STAR Methods). Bottom: (1) Using the set of unproductive inferred naive BCRs, we inferred a model to characterize the null probability for generation of receptors $P_{\mathrm{gen}}\left(\sigma \right)$(Marcou et al., 2018). We inferred a selection model (Sethna et al., 2020) to characterize the deviation from the null among inferred naive productive BCRs, with the probability of entry to the periphery $P_{\mathrm{post}}\left(\sigma \right)$ and selection factors $q_f\left(\sigma \right)$ dependent on receptor sequence features. (2) Based on temporal information of sampled BCRs, we identified clonal lineages that expanded significantly during infection. (3) We identified progenitors of clonal lineages shared among individuals and assessed the significance of these sharing statistics based on the probabilities to find each receptor in the periphery. The shared, expanding clonal lineages that contain plasma B cells are likely candidates for secreting responsive antibodies during infection. We verified the reactivity of receptors to SARS-CoV-2 antigenic epitopes using sorted single-cell data. We also identified previously characterized monoclonal antibodies (mAbs) specific to SARS-CoV-2 and SARS-CoV-1.

Figure 2

Sequence features of immune receptors in the bulk repertoire across cohorts (A) The relative counts for IGHV gene usage for inferred naive progenitors of clonal lineages in healthy individuals and those with mild, moderate, and severe COVID-19 symptoms. The bars indicate the use frequency averaged over individuals in each cohort, and dots indicate the variation in V gene frequencies across individuals (biological replicates) within each cohort. (B and C) Statistics of length of HCDR3 amino acid sequence for different individuals (biological replicates) in each cohort. The violin plots in (B) show the mean HCDR3 length of each individual (dots) in a given cohort (color), with the violin plot cut parameter set to 0.1. The mean HCDR3 lengths of the sorted single cells and verified mAbs (axis) for RBD-reactive (pink squares) and NTD-reactive (purple plus symbols) receptors are shown on the right. Full lines in (C) show distributions averaged over individuals (biological replicates) in each cohort (color), and shading indicates regions containing one standard deviation of variation across individuals within a cohort. One-way ANOVA statistical tests were performed, comparing the mean HCDR3 of all COVID-19 cohorts and the healthy repertoires from the Great Repertoire Project (GRP) dataset (Briney et al., 2019) with the healthy control from this study: healthy-mild: F_1,3 = 12.0, p = 0.04; healthy-moderate: F_1,13 = 15.7, p = 0.0016; healthy-severe: F_1,6 = 37.5, p = 0.00087; healthy-GRP: F_1,11 = 0.9, p = 0.359. Significance cutoffs: n.s. p > 0.01, ^∗p ≤ 0.01, ^∗∗p < 0.001. (D) The relative counts for IGHJ gene usage for inferred naive progenitors of clonal lineages in cohorts of healthy individuals and COVID-19 cohorts of individuals with mild, moderate, and severe symptoms. The bars indicate the use frequency averaged over individuals in each cohort, and dots indicate the variation in J gene frequencies across individuals (biological replicates) within each cohort. See Data S1 for details regarding biological replicates.

Figure 3

Differential statistics of immune repertoires across cohorts (A) The distribution of the log probability to observe a sequence $\sigma$ in the periphery ${\log }_{10}{P}_{\mathrm{post}}\left(\sigma \right)$ is shown as a normalized probability density function (PDF) for inferred naive progenitors of clonal lineages in cohorts of healthy individuals and the mild, moderate, and severe cohorts of individuals with COVID-19. Full lines show distributions averaged over individuals (biological replicates; Data S1) in each cohort, and shading indicates regions containing one standard deviation of variation among individuals within a cohort. (B) Clustering of cohorts based on their pairwise Jensen-Shannon divergences (D_JS) as a measure of differential selection on cohorts (STAR Methods). (C) The bar graph shows how incorporating different features into a SONIA selection model contributes to the fractional D_JS between models trained on different cohorts. The error bars show the standard deviation of these estimates, using five independent sets of 100,000 generated BCRs for each selection model (STAR Methods). (D–F) Logo plots show the expected differences in the log-selection factors for amino acid usage, $⟨\text{Δ}\log Q_{\mathrm{cohort}}\left(a\right)⟩=⟨\log Q_{\mathrm{cohort}}\left(a\right)-\log Q_{\mathrm{healthy}}\left(a\right)⟩$, for the (D) mild, (E) moderate, and (F) severe COVID-19 cohorts. The expectation values $⟨\cdot ⟩$ are evaluated on the mixture distribution $\frac{1}{2}\left({\mathit{P}}_{\mathrm{post}}^{\mathrm{cohort}}+{\mathit{P}}_{\mathrm{post}}^{\mathrm{healthy}}\right)$. Positively charged amino acids (lysine, K; arginine, R; and histidine, H) are shown in blue, and negatively charged amino acids (aspartate, D, and glutamate, E) are shown in red. All other amino acids are shown in gray. Positions along the HCDR3 are shown up to 10 residues starting from the 3′ (positive values) and 5′ ends (negative values). (G) The bar graph shows the average mean difference between the log-selection factors for IGHV gene usage for the mild (green), moderate (yellow), and severe (red) COVID-19 cohorts, with the mean differences computed using the mixture distribution $\frac{1}{2}\left({\mathit{P}}_{\mathrm{post}}^{\mathrm{cohort}}+{\mathit{P}}_{\mathrm{post}}^{\mathrm{healthy}}\right)$, and the average is taken over the 30 independently trained SONIA models for each cohort. Error bars show standard deviation of these estimates across the inferred SONIA models (STAR Methods).

Figure 4

Dynamics of BCR repertoires during infection (A) The binding level (measured by OD₄₅₀ in ELISA) of the IgM (left) and IgG (right) repertoires to SARS-CoV-2 (RBD) epitopes increases over time in most individuals. (B) The log ratio of BCR (mRNA) abundance at late time versus early time is shown for all clonal lineages that are present in at least two time points (STAR Methods). Each panel shows dynamics of lineages for a given individual, as indicated in the label. The analysis is shown for individuals for whom the binding level (OD₄₅₀) of the IgG repertoire increases over time (shown in A). The count density indicates the number of lineages at each point. Lineages that show a significant expansion over time are indicated in red (STAR Methods). (C) IGHV gene use of lineages is shown for non-expanded (left) and expanded (center) lineages in all individuals (colors). The right panel shows, for each individual (colors), the fraction of expanded lineages with a given IGHV gene as the number of expanded lineages divided by the total number of lineages with that given IGHV gene. The size of the circles indicates the total number of lineages in each category. (D) Boxplot of log₁₀ relative read abundance in the plasma B cell repertoire (STAR Methods) for expanding (red) and non-expanding (cyan) lineages that contain reads from plasma B cells in different individuals. Receptors from the plasma B cell repertoire are significantly more abundant in expanding lineages in four individuals based on ANOVA test statistics: individual 3: F_1,42 = 5.4, p = 0.02; individual 5: F_1,31 = 0.5, p = 0.5; individual 7: F_1,49 = 0.01, p = 0.91; individual 9: F_1,42 = 4.1, p = 0.04; individual 10: F_1,42 = 2.9, p = 0.1; individual 13: F_1,64 = 7.7, p = 0.007.

Figure 5

Sharing of BCRs among individuals (A) The histogram shows the number of clonal lineages that share a common progenitor in a given number of individuals, indicated on the horizontal axis. (B) The density plot shows the distribution of ${\log }_{10}{P}_{\mathrm{post}}$ for progenitors of clonal lineages shared in a given number of individuals, indicated on the horizontal axis. The histogram bin size is 0.5. The scaling of sequence counts sets the maximum of the density in each column to one. Sharing of rare lineages with ${\log }_{10}{P}_{\mathrm{post}}$ below the dashed line is statistically significant (STAR Methods). Green diamonds indicate clonal lineages below the dashed line with significant expansion in at least one of the individuals. Orange triangles indicate clonal lineages below the dashed line that contain reads from the plasma B cell repertoire in at least one of the individuals. (C and E) Histograms showing the numbers of clonal lineages that share a common progenitor in a given number of individuals that have expanded significantly during infection in at least one of the individuals (C) or contained reads from the plasma B cell repertoire in at least one of the individuals (E). (D and F) Scatterplots with transparent overlapping markers show ${\log }_{10}{P}_{\mathrm{post}}$ for progenitors of clonal lineages shared in a given number of individuals that have expanded (D) or contain reads from the plasma B cell repertoire (F) in at least one individual. The dashed line is similar to (B).

Figure 6

Statistics of BCRs reactive to RBD and NTD epitopes (A) Relative counts for IGHV gene usage for known mAbs (Data S3) reactive to RBD (pink) and NTD (green) epitopes of SARS-CoV-2 and for receptors obtained from single-cell sequencing of the pooled sample from all individuals (STAR Methods), sorted for RBD (yellow) and NTD (blue) epitopes. (B) The histogram shows the number of NTD-sorted receptors from single cell sequencing (Data S2) and RBD- and NTD-specific verified mAbs (Data S3) found in the bulk+plasma B cell repertoires of a given number of individuals (STAR Methods), indicated on the horizontal axis. (C) The distribution of the log probability to observe a sequence $\sigma$ in the periphery ${\log }_{10}{P}_{\mathrm{post}}\left(\sigma \right)$ is shown as a normalized PDF for inferred naive progenitors of known RBD- and NTD-specific mAbs and for RBD- and NTD-sorted receptors from single-cell sequencing. $P_{\mathrm{post}}\left(\sigma \right)$ values were evaluated based on the repertoire model created from individuals with moderate symptoms. The corresponding ${\log }_{10}{P}_{\mathrm{post}}$ distribution for bulk repertoires of the moderate cohort (similar to Figure 3A) is shown in black as a reference. (D) Similar to (C) but restricted to receptors that are found in the bulk+plasma repertoire of at least one individual in the cohort (Data S2 and S3). Colors are consistent between panels, and the number of samples used in each panel is indicated in the legend.