Deep immunoglobulin repertoire sequencing depicts a comprehensive atlas of spike-specific antibody lineages shared among COVID-19 convalescents

ABSTRACT

Neutralizing antibodies are a key component in protective humoral immunity against SARS-CoV-2. Currently, available technologies cannot track epitope-specific antibodies in global antibody repertoires. Thus, the comprehensive repertoire of spike-specific neutralizing antibodies elicited by SARS-CoV-2 infection is not fully understood. We therefore combined high-throughput immunoglobulin heavy chain (IgH) repertoire sequencing, and structural and bioinformatics analysis to establish an antibodyomics pipeline, which enables tracking spike-specific antibody lineages that target certain neutralizing epitopes. We mapped the neutralizing epitopes on the spike and determined the epitope-preferential antibody lineages. This analysis also revealed numerous overlaps between immunodominant neutralizing antibody-binding sites and mutation hotspots on spikes as observed so far in SARS-CoV-2 variants. By clustering 2677 spike-specific antibodies with 360 million IgH sequences that we sequenced, a total of 329 shared spike-specific antibody clonotypes were identified from 33 COVID-19 convalescents and 24 SARS-CoV-2-naïve individuals. Epitope mapping showed that the shared antibody responses target not only neutralizing epitopes on RBD and NTD but also non-neutralizing epitopes on S2. The immunodominance of neutralizing antibody response is determined by the occurrence of specific precursors in human naïve B-cell repertoires. We identified that only 28 out of the 329 shared spike-specific antibody clonotypes persisted for at least 12 months. Among them, long-lived IGHV3-53 antibodies are likely to evolve cross-reactivity to Omicron variants through accumulating somatic hypermutations. Altogether, we created a comprehensive atlas of spike-targeting antibody lineages in COVID-19 convalescents and antibody precursors in human naïve B cell repertoires, providing a valuable reference for future vaccine design and evaluation.

KEYWORDS:

• SARS-CoV-2
• spike
• antibodyomics
• antibody repertoire
• epitope mapping
• lineage tracking

Introduction

The ongoing coronavirus disease 2019 (COVID-19) pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has already resulted in over 771 million reported cases and more than 6.9 million deaths worldwide as of 17 October 2023 (https://covid19.who.int). The development of more effective vaccines may represent one potential avenue for alleviating the COVID-19 pandemic. The surface spike (S) glycoprotein of SARS-CoV-2 engages the cellular receptor angiotensin-converting enzyme 2 (ACE2) via the receptor-binding domain (RBD), which is believed to be a critical target to block viral entry [1]. Spike, especially RBD, is also the major target for inducing neutralizing antibodies and has been commonly used as an immunogen in various vaccines [2]. Neutralizing antibodies targeting the N-terminal domain (NTD) and S2 have also been reported [3, 4]. Neutralizing antibody levels are highly predictive of immune protection from SARS-CoV-2 infection or vaccination [5, 6]. In response to gradually established population immunity, SARS-CoV-2 continues to evolve, generating numerous Variants of Concern (VOCs) or Variants of Interest (VOIs) [7, 8]. Emerging evidence supports that neutralizing antibody response shared among populations drives the genesis of SARS-CoV-2 variants [7, 9, 10], indicating convergent selection pressure on viral evolution. Thus, a comprehensive understanding of the magnitude, breadth, and longevity of spike-specific antibody responses in naturally infected populations could be of significance for guiding the development of vaccines.

Given the critical role of the spike-specific antibody response in protective immunity against SARS-CoV-2 infection, this area is subject to intense research providing valuable knowledge for vaccine improvement [11, 12]. However, these studies are based primarily on serological and cell biological approaches. Specific antibody responses were analyzed at the global (spike) level or low resolution at the subdomain (RBD, NTD, or S2) level, making it difficult to identify effective antibody components that contribute to protective immunity. Previously, attempts were made to use advanced serology to map SARS-CoV-2 spike epitopes using linear peptides [13–15], these approaches are only applicable to investigate linear epitopes and likely to miss conformational epitopes. Therefore, the comprehensive repertoire of spike-specific neutralizing antibodies elicited by SARS-CoV-2 infection needs to be well defined. To overcome the above obstacle, we and other groups have started to profile antibody repertoire using high-throughput immunoglobulin heavy chain repertoire sequencing (IR-seq), providing a unique perspective for the systemic B cell response after SARS-CoV-2 infection [16–18]. In this study, we further improved the downstream analysis methodology of IR-seq by combining it with structural and bioinformatics analysis. We established an antibodyomics method that enables quantifying the abundance and prevalence of antibody lineages targeting certain neutralizing epitopes on the SARS-CoV-2 spike protein. After running the whole pipeline, we obtained a comprehensive repertoire of spike-specific neutralizing antibody responses in COVID-19 convalescent patients. This full repertoire gives in-depth insights into neutralizing antibody-mediated protective immunity against SARS-CoV-2.

Materials and methods

Study cohort and blood sample collection

The cohort in this study included 33 COVID-19 convalescents and 24 SARS-CoV-2-naïve individuals. Some donors in this cohort have been previously described [17, 19, 20]. An additional 23 donors were newly included in this study. COVID-19 participants were recruited between 31 March 2020 and 20 July 2020 from patients attending Guangzhou Eighth People’s Hospital with a confirmed diagnosis of COVID-19. Healthy controls were recruited before the COVID-19 pandemic and had a negative serology result. This study was approved by the Review Committee of Guangzhou Eighth People’s Hospital of Guangzhou Medical University. Blood samples from COVID-19 inpatients were sampled at study entry and then 1–3 time points before discharge (approximately within 4 weeks). Discharged patients were invited to provide a follow-up sample 3–12 months after study enrollment. In total, 78 blood samples were collected and used for further IR-seq. These blood samples from COVID-19 convalescents were classified into 4 groups: Visit 1 (22 samples, within 2 weeks post symptom onset), Visit 2 (26 samples, 2–4 weeks post symptom onset), Visit 3 (10 samples, 3–6 months post symptom onset), and Visit 4 (20 samples, 12 months post symptom onset).

Immunoglobulin heavy chain library preparation and sequencing

Peripheral blood mononuclear cells (PBMCs) were isolated from blood samples using Opti-Prep lymphocyte separation solution (Axis Shield Poc As, Norway) following the manufacturer’s instructions. Total RNA was extracted using TRIzol^TM LS reagent according to the manufacturer’s protocol (Life Technologies). A iR-RepSeq-plus Heavy-Chain Cassette based on dam-PCR (iRepertoire, Inc., Huntsville, AL, USA) was used to amplify the IgH repertoire sequences as previously described [20, 21]. Briefly, multiplex primers covering the human IGHV genes (forward primers) and constant region primers (reverse primers) were designed. The forward primers Fi (forward-in) and reverse primers Ri (reverse-in) also included Illumina paired-end sequencing communal primers, respectively (Illumina, USA). Unique barcodes were introduced in the first round by the constant region primers. After gel purification using a QIAquick® gel extraction kit (Cat No. 28,704; Qiagen), the product was pooled and sequenced on an Illumina NovaSeq 6000 with paired-end 250 bp read mode (Novogene, China).

Bioinformatics analysis

Processing of raw sequencing data was first performed. Reads were filtered for base quality, and those with scores over 20 at 3’ ends were retained using Trimmomatic (v 0.36) [22]. After filtering, the paired-end reads were separated based on the unique barcodes at the 5’ end of the reads. The separated reads were merged using FLASH (v 1.2.11) [23] if their overlapping regions were more than 30 bp. Finally, the whole-length antibody sequences with 300–470 bp were annotated within the V(D)J germline genes using MIXCR (v 3.0.3) [24], and the reference V(D)J sequences were downloaded from the IMGT database (http://www.imgt.org/). Somatic hypermutation rates (SHM) were defined as the frequency of nucleotide substitutions in the V gene region (FR1-FR3) compared to putative germline genes. Productive sequences with the same IGHV gene, the same IGHJ gene, and identical HCDR3 (amino acid) were defined as one IgH clone (Figure S5). Clone size was normalized to the number of reads for a clone per 1,000,000 reads in each IgH repertoire to correct for bias in sequencing depth across samples. To identify spike-specific lineages, a list of known SARS-CoV-2 mAbs was curated manually from 171 studies (Table S1). All these mAbs required well-defined IGHV/IGHJ genes and complete HCDR3 sequences. After compiling, a total of 2677 spike-specific mAbs were obtained for downstream analysis. Using the heavy chains of collected spike-specific mAbs as references to compare with IgH clones from COVID-19 convalescents and healthy donors, we defined the IgH clones that have the same IGHV and IGHJ germline genes and their HCDR3 amino acid similarity ≥80% (hamming distance) with the known spike-specific mAbs as clonally related spike-specific lineages. Here, the hamming distance was calculated using R package stringdist (v 0.9.10). The abundance of a spike-specific lineage was calculated by summing the size of clones that clonally related to the given spike-specific mAb in each IgH repertoire. The clonally related B-cell lineages were then clustered into the same clonotype if they have the same IGHV and IGHJ germline genes and their HCDR3 amino acid similarity ≥80%, which was considered to recognize similar epitopes [25–27]. If a clonotype contains clones from different individuals, it was defined as a shared clonotype. The antigen specificity of each clonotype was annotated according to the antigen specificity of the reference mAb that belongs to the same clonotype as itself. For example, if the reference mAb of a clonotype is an RBD-targeting neutralizing mAb, this clonotype is therefore defined as RBD-targeting neutralizing clonotype. The prevalence of a spike-specific clonotype was calculated by dividing the number of individuals detected in the COVID-19 or healthy group by the total number of COVID-19 or healthy individuals. The identification of Omicron BA.1 cross-reactive clonotypes followed the same pipeline, and another 799 mAbs (Table S7) isolated from Omicron BA.1 breakthrough-infected donors were included for this analysis. All the visualization was conducted in R platform using packages including ggplot2 (v 3.3.3), GraphPad Prism (v 8.0.1), ggseqlogo (v 0.1), circlize (v 0.415), and ComplexHeatmap (v 2.14.0).

Unsupervised classification of neutralizing epitopes and spike sequence entropy analysis

Human SARS-CoV-2 spike-specific mAbs with available structures were downloaded from the PDB (https://www.rcsb.org/). A total of 251 SARS-CoV-2 antigen–antibody complexes were collected. Epitope residues and buried surface area (BSA) for each epitope residue of the 251 spike-specific mAbs were determined using the PDBePISA server (https://www.ebi.ac.uk/msd-srv/prot_int/). BSA for each epitope residue is considered a feature of a certain antibody and used to construct a feature matrix M_AxB for downstream analysis, where A is the number of antibodies and B is the number of features (amino acid length of RBD, NTD, or S2). Therefore, three BSA matrixes for 221 RBD mAbs, 25 NTD mAbs, and 5 S2 mAbs were obtained, which were subsequently used as inputs for epitope classification with the R package UMAP (v 0.2.9.0). Uniform Manifold Approximation and Projection (UMAP) is an algorithm for dimensional reduction. Next, antibodies were further embedded into two-dimensional space for visualization with UMAP and clustered into different epitope groups. Clustering was performed using the K-means algorithm. UMAP and K-means clustering were conducted independently. All UMAP plots were generated by the R package ggplot2 (v 3.3.3). Per-site variation of SARS-CoV-2 spike was calculated by sequence entropy using the COVID-19 Viral Genome Analysis Pipeline (https://cov.lanl.gov/components/sequence/COV/int_sites_tbls.comp?t = 2) with default parameters. A total of 15,489,226 SARS-CoV-2 sequences that deposited in the Global Initiative for Sharing Avian Influenza Database (GISAID; https://gisaid.org) from 1 January 2020 to 17 October 2023 were included in this analysis.

Quantification and statistical analysis

Statistical analysis was performed in R or GraphPad Prism. Statistical tests for each analysis are indicated in the respective figure legends. P values were considered to be significant when < 0.05. Error bars represent the standard deviation (SD).

Results

Overview of the antibodyomics pipeline

Firstly, 2677 published human-derived spike-specific mAbs were collected (Table S1), 251 neutralizing mAbs of them with determined structures (Table S2) were used for epitope classification using an unsupervised clustering algorithm. Then the neutralizing epitopes on the spike were mapped with the substitutions in SARS-CoV-2 variants. Secondly, PBMC from 33 COVID-19 convalescents and 24 healthy controls were used for IgH repertoire sequencing. More than 360 million IgH sequences were obtained. Thirdly, the sequenced IgH sequences were clustered with the published spike-specific mAbs to identify spike-specific antibody lineages. Lastly, the identified spike-specific antibody lineages were used for downstream analysis, including identification of shared antibody lineages, epitope mapping, and persistent or cross-reactive antibody lineage tracking (see methods, Figure 1).

Figure 1. Schematic diagram of the antibodyomics pipeline. This pipeline is based on the combination of high-throughput immunoglobulin heavy chain (IgH) repertoire sequencing, structural and bioinformatics analysis. Firstly, the published spike-specific mAbs were collected, which of them have determined structures were used for epitope classification. Secondly, PBMC from COVID-19 convalescents and healthy controls were used for IgH repertoire sequencing. Then the sequenced IgH sequences were clustered with the published spike-specific mAbs to identify spike-specific antibody lineages. Lastly, the identified spike-specific antibody lineages were used for downstream analysis, including identification of shared antibody lineages, epitope mapping, persistent or cross-reactive antibody lineage tracking.

High-resolution structure-based unsupervised classification of spike epitopes

To obtain a full spectrum of neutralizing epitopes known on the SARS-CoV-2 spike, 2677 spike-targeting mAbs from 171 publications were compiled and analyzed, including 1792, 167, and 161 targeting the RBD, NTD, and S2, respectively (Figure S1A, Table S1). The IGHV usage of RBD, NTD, and S2 mAbs showed distinct profiles (Figure S1B). Among them, 251 experiment-confirmed human neutralizing mAbs with characterized structures were used for epitope mapping (Table S2). According to the unsupervised clustering, 251 spike-targeting mAbs can be divided into 12 major classes: RBD I-VII, NTD I-III, and S2 I-II (Figure 2A–C, S1–3).

Figure 2. Epitope mapping of SARS-CoV-2 RBD antibodies. (A) UMAP and unsupervised clustering of 221 RBD-targeting mAbs. Seven epitope groups were identified (left). Pie charts show the epitope distribution of mAbs encoded by different IGHVs (right). (B) IGHV germline gene usage was projected onto UMAP (left). Pie charts show the IGHV usage distribution of different classes of mAbs (right). (C) Heatmap shows epitope residues of ACE2 and 7 classes of RBD epitopes. The heatmap is divided into 4 blocks by row. The first block shows the epitope residues of ACE2 and the percentage of each class of mAbs that have interacted with each epitope residue on the RBD. The second block shows the buried surface area (BSA) of ACE2 and the average BSA of each class of mAbs at each epitope residue on the RBD. The third block summarizes the per-site sequence entropy of the RBD. The bottom block shows the substitution sites in the RBD region. (D) Footprints of ACE2 and 7 classes of epitopes on RBD. RBD surface presentation is shown as different orientations. Epitope residues are colored based on the percentage of mAbs that have interacted with these sites.

We defined 3 anatomic surfaces on RBD, namely, the outer surface (orientation 1), the inner surface (orientation 2), and ACE2 binding surface (orientation 3 and 4) to help understand their molecular properties (Figure 2D). The RBD I and RBD II epitopes are targeted by antibodies that are predominantly encoded by IGHV3-53/3-66 and IGHV1-58, respectively, which largely overlap with the ACE2-binding site and are accessible only when the RBD is in the “up” conformation (Figure 2B and D). Compared to the RBD I epitope, the footprint of the RBD II epitope is less extensive appearing to cover only part of the ACE2-binding site towards the tip of the RBD. While RBD III and IV epitopes overlap significantly, RBD III buries more inner surface (orientation 2) than RBD IV. Although RBD III and IV epitopes overlap with the ACE2-binding site (Figure 2D), they are accessible in both “up” and “down” conformations for antibody binding. Most antibodies targeting the RBD III epitope are mainly derived from IGHV1-2, IGHV1-69, IGHV3-30, and IGHV3-53, accounting for 64.7% of total RBD III-directed mAbs (Figure 2B). RBD V is a semi-cryptic epitope when the RBD adopts the “down” conformation, which is located opposite the ACE2-binding site (Figure 2D). The footprint of the RBD VI epitope is solvent-exposed on the RBD, accessible in both “up” and “down” conformations. RBD VII is a cryptic epitope that faces toward the center of the trimer, whose exposure requires the RBD in the “up” conformation (Figure 2D). Most of the RBD VII mAbs shared epitope residues 369–385 (Figure 2C, S1C), which show broad cross-reactivity across sarbecoviruses [28, 29]. The neutralizing epitopes on the NTD are relatively limited and concentrated. The 25 structurally characterized NTD-targeting mAbs can be divided into 3 classes (Figure S1C, S2A). NTD I and NTD II epitopes serve as a neutralizing supersite [3]. Quantitatively, all 14 NTD I and II mAbs targeting the site of vulnerability are encoded by 5 germline genes (Figure S2B). However, the NTD III epitope defines distinct regions around the supersite (Figure S2C and 2D). In contrast to the RBD and NTD, only two linear neutralizing epitopes on S2 have been found thus far. One is the fusion loop (epitope residues 812–825) and the other is a conserved stem helix (epitope residues 1143–1158) (Figure S3A and 3B). Collectively, a total of 12 classes of neutralizing epitopes were identified on the SARS-CoV-2 spike, and preferential IGHVs were used for targeting certain neutralizing epitopes.

SARS-CoV-2 spike mutation hotspots highly overlap with immunodominant neutralizing epitopes

We next mapped the substitutions in SARS-CoV-2 variants to the identified neutralizing epitopes on spike. According to 15,489,226 SARS-CoV-2 genome sequences deposited in the GISIAD database from 1 January 2020 to 17 October 2023, we analyzed the per-site variation in spikes by calculating sequence entropy (Figure 2C, S2C, S3A). On the RBD, a total of 34 substitutions emerged in the prevailing SARS-CoV-2 VOI or VOCs. These were defined as 5 substitution hotspots (Figure 2C). Notably, most of these substitutive positions largely interact with neutralizing antibodies (Figure S4A). Among them, several substitutions including R346, L452, L455, F456, E484, and Q493 are located within multiple epitopes (Figure 2D), allowing them to escape a broad range of neutralizing antibodies [7, 30]. On the NTD, 5 substitution hotspots were identified (Figure S2C). Most of the substitutive residues are located at the NTD I and NTD II defined supersite (Figure S2D, S4B). Unlike the RBD and NTD, all substitutions on S2 are outside of the two neutralizing epitopes except for substitution at residue P1143 that observed on recently emerged BA.2.86, which located at the boundary of stem helix epitope (Figure S3A). Overall, the above results reveal that the mutation hotspots, particularly RBD substitutions, highly overlap with the immunodominant neutralizing mAb-directed residues, highlighting strong neutralizing antibody-mediated immune pressure.

Spike-specific igH lineages are commonly present in both SARS-CoV-2-exposed and -unexposed antibody repertoires

To obtain a full repertoire of spike-specific antibody lineages induced by SARS-CoV-2, longitudinal PBMC samples collected within two weeks to one-year post symptom onset (Visit 1 to Visit 4: Visit 1 (within 2 weeks), Visit 2 (2–4 weeks), Visit 3 (3–6 months), and Visit 4 (12 months)) from a cohort of 33 COVID-19 convalescents and 24 healthy donors (Table S3) were subjected for IgH sequencing. The 2677 spike-targeting mAbs (Table S1) were used as references to align with the sequenced IgH repertoires to identify spike-specific lineages (Figure S5). Here, the IgH sequences encoded by the same IGHV and IGHJ genes and their HCDR3 amino acid similarity ≥80% with the known spike-specific mAbs are considered as clonally related lineages [19, 31, 32], which generally recognize similar epitopes [25, 26, 27] (Figure S5). With this parameter, 591 spike-targeting mAbs were found to have clonally related lineages in COVID-19 patients’ IgH repertoires (Figure 3A and -B). Meanwhile, 191 spike-targeting mAbs were found to have clonally related lineages in SARS-CoV-2-unexposed IgH repertoires (Figure 3A). Notably, 157 mAbs were found to have clonally related lineages in both SARS-CoV-2-exposed and -unexposed IgH repertoires, indicating that precursors of spike-specific antibodies are commonly present in the repertoires of SARS-CoV-2-naïve individuals (Figure 3A and B). Indeed, we found that most of the detected spike-specific IgH lineages in SARS-CoV-2-unexposed repertoires are IgM (59.2%), whereas the major isotype in SARS-CoV-2-exposed individuals are IgG (43.6%). Class-switched IgG and IgA accounted for 62.9% in SARS-CoV-2-exposed individuals but only 23.5% in unexposed donors (Figure 3C). This result indicates that spike-reactive naïve B cells expressing IgM switched to IgG⁺ or IgA⁺ antibody-expressing cells after encountering spike antigens. Analysis of the somatic hypermutation (SHM) of heavy chain V gene region of the 2677 spike-targeting mAbs revealed that the spike-targeting mAbs generally carried limited SHMs, with a mean SHM of 3.88 ± 0.11% (Figure 3D). Consistent with these isolated mAbs, their clonally related lineages identified in the SARS-CoV-2-exposed repertoires also showed limited SHM (Figure 3D). The above observations suggest that the majority of spike-specific antibody lineages elicited by SARS-CoV-2 infection are derived from naïve B cells that are constitutively present in a broad human population.

Figure 3. Identification of spike-specific IgH lineages. (A) Heatmap shows the 625 spike-specific IgH lineages. The heatmap is divided into 4 parts by row according to the different spike domains. The top annotation represents grouping information (Visit 1 to Visit 4). Left annotation represents IGHV usage of the spike-specific IgH lineages. Right annotation represents the prevalence of each lineage. The most prevalent lineages are marked with their related mAbs. (B) Venn diagram shows the number of overlapped spike-specific IgH lineages between SARS-CoV-2-exposed and -unexposed individuals. (C) Comparison of antibody isotype distribution of the spike-specific IgH lineages in SARS-CoV-2-exposed and -unexposed individuals. A two-tailed chi-square test was performed (*** p < 0.001). (D) Comparison of somatic hypermutation (SHM) of the spike-specific IgH lineages in different antibody isotypes in SARS-CoV-2-exposed and -unexposed individuals. SHM of the known spike-specific mAbs is also shown. The Wilcoxon rank-sum test was performed (** p < 0.01, *** p < 0.001).

Shared antibody responses target not only neutralizing epitopes on the RBD and NTD but also non-neutralizing epitopes on the S2

Next, we clustered the spike-specific mAbs with their clonally related IgH sequences to define shared IgH clonotypes (Figure S5). We found that the 591 spike-targeting mAbs-related lineages that present in COVID-19 patient’s repertoires (Figure 3A and B) can be clustered into 329 clonotypes (Figure 4A, Table S4). Of the 329 clonotypes, 202 targeted the RBD, 38 targeted the NTD, 33 targeted the S2, and 56 targeted the undefined domain of spike (S-undefined) (Figure 4B, Table S4). We summarized all the clonotypes in Table S4 and described the representative clonotypes in detail below. Notably, the top 8 highly shared clonotypes (cl.95, cl.276, cl.142, cl.143, cl.144, cl.220, cl.43, and cl.44) were detectable in more than 35% of COVID-19 donors (Figure 4A, Table S4), representing the most common IgH lineages against SARS-CoV-2. Spike-specific IgH clonotypes with identical HCDR3s were frequently observed in different COVID-19 donors and even in some SARS-CoV-2-unexposed individuals (Figure 4C). Clonotype cl.95 encoded by IGHV3-7 is derived from the same lineage as S2-targeting non-neutralizing mAbs such as COV2-2333 and COV2-2002, which recognize the heptad repeat (HR1) region of S2 [33] (Figure 4A and Figure 4C). Clonotype cl.276 encoded by IGHV3-15 is derived from the same lineage as the non-neutralizing NTD-targeting mAbs CV21 [34], DH1118, and DH1119 [4] (Figure 4A and C). The most common clonotype cl.142, and two other clonotypes cl.143 and cl.144, encoded by IGHV3-30, are S2-targeting non-neutralizing antibody lineages [34] (Figure 4A and C; Figure S6). IGHV3-30 also encodes a prevalent RBD-targeting clonotype cl.220, which is related to the neutralizing mAbs CC12.17, COV2-2007, and G32B1 (Figure 4A; Figure S5). Of note, the most common RBD-targeting antibody clonotypes cl.43 and cl.44 are both encoded by IGHV3-53, which correspond to the two well-documented public clonotypes described previously in our and other studies and recognize the RBD I epitope (Figure 2A) [19, 27]. In addition, another 24 shared clonotypes were observed in more than 20% of COVID-19 donors in our cohort, including 13 RBD-targeting neutralizing clonotypes, 2 NTD-targeting neutralizing clonotypes, 6 S2-targeting non-neutralizing clonotypes, and 3 S-undefined clonotypes (Figure 4A and B, Table S4). Among them, clonotype cl.86 relates to RBD I mAbs BD-236, BG4-25, CV30, and E4, representing the third most common clonotype after cl.43 and cl.44 encoded by IGHV3-53 (Figure 4A, S6). Clonotype cl.209 encoded by IGHV1-58 relates to mAbs 13G9, COV2-2196, S2E12, ect that recognize the RBD II epitope (Figure 2A, S6). Clonotype cl.208 encoded by IGHV1-24 relates to the NTD-supersite (NTD I) targeting mAb N11 and COVA2-25 (Figure 4A, S6). Collectively, this above observation indicates that SARS-CoV-2 infection also induced massive non-neutralizing antibody responses.

Figure 4. Identification of shared spike-specific clonotypes. (A) Population prevalence and IGHV usage of the 329 shared clonotypes. The size of the dot indicates the number of clonotypes. 11 representative clonotypes are marked. (B) Distribution of shared clonotypes that target different spike domains. The bars are partitioned and colored based on the prevalence of each clonotype. (C) Circos plots and logo plots show the clonal relationship of the top 8 highly shared clonotypes. Interconnecting links within circos plots indicate the relationship between IgH sequences that share IGHV and IGHJ genes and have identical HCDR3 aa sequences. Green links represent clonally related IgH sequences shared by two SARS-CoV-2-exposed individuals, and orange links represent clonally related IgH sequences shared by SARS-CoV-2-exposed and -unexposed individuals. The width of the link is determined by the number of unique IgH sequences shared by two individuals. (D) Logo plots show HCDR3 sequences of each clonotypes. Representative published mAbs clonally related to the clonotypes are listed below the logo plots.

The immunodominance of the neutralizing antibody response against SARS-CoV-2 is determined by the occurrence of specific precursors in human naïve B-cell repertoires

To further define the prevalence of antibody lineages that target different classes of neutralizing epitopes on spikes, we characterized the expression profile of 251 mAbs with determined structures and clustered them with the identified spike-specific antibody clonotypes (Figure 4A). Overall, 77 out of the 251 epitope-defined human mAbs showed population-level convergence (Figure 5A, Table S5). The 77 mAbs corresponded to 44 shared clonotypes according to our categorization. The population prevalence of those spike-specific antibody clonotypes that target different epitopes showed distinct profiles (Figure 5A). IGHV3-53-encoded cl.43 and cl.44 targeting the RBD I epitope and IGHV1-58-encoded cl.209 targeting the RBD II epitope represent the most prevalent neutralizing antibody lineages against SARS-CoV-2 (Figure 5A). IGHV1-69-encoded cl.210 targeting RBD III represents the third most common neutralizing antibody lineage, corresponding to our previously identified L452-contacting antibodies [35]. However, antibody responses targeting other RBD epitopes were relatively less prevalent, presenting in 3.03% to 21.21% of COVID-19 convalescents (Figure 5A, Table S5). Notably, we found that a shared clonotype (cl.138) related to SARS-CoV-1 cross-reactive mAb S304 that directs the RBD VII epitope can be detected in 15.15% of COVID-19 convalescents. N11-related cl.208 encoded by IGHV1-24 targeting NTD I is the most prevalent neutralizing clonotype on NTD, presenting in 24.24% of COVID-19 convalescents. The prevalence of other NTD mAbs ranged from 3.03% to 9.09% (Figure 5A). No convergence was observed for the mAbs that target the conserved S2-neutralizing epitopes (Figure 5A).

Figure 5. Abundance and prevalence of the epitope-specific antibody lineages. (A) Heatmap shows the abundance and population prevalence of 251 epitope-defined mAbs and their 160 clonally related published spike-specific mAbs. Representative clonotypes are marked with their related mAbs on the right. (B) Abundance of mAbs approved for emergency use authorization or in clinical trials in SARS-CoV-2-exposed (green) and -unexposed (orange) individuals. (C) Abundance of ACE2-competitive (blue) or non-ACE2-competitive (dark blue) spike-specific IgH lineages in each sample of SARS-CoV-2-exposed and -unexposed individuals. Pie charts on the top panel show the prevalence of ACE2-competitive or non-ACE2-competitive spike-specific IgH lineages in SARS-CoV-2-exposed and -unexposed individuals. (D) Abundance of spike-specific IgH lineages that cross-reactive with SARS-CoV-1 (orange) or not (red) in each sample of SARS-CoV-2-exposed and -unexposed individuals. Pie charts on the top panel show the prevalence of spike-specific IgH lineages that cross-reactive with SARS-CoV-1 in SARS-CoV-2-exposed and -unexposed individuals. The one-sided chi-square test was performed in panels (C) and (D) (*** p < 0.001, * p < 0.5, ns p ≥ 0.5).

We found there is no population prevalence for most of the mAbs approved for emergency use authorization or in clinical trials, except for P2B-1F11 (BRII-196), COV2-2196, and DXP604 (Figure 5B), suggesting that ultrapotent neutralizing antibodies generally do not belong to the multi-donor class. In total, the abundance of ACE2-competitive clonotypes was much higher than none-ACE2-competitive clonotypes (477.86 ± 139.99 versus 326.02 ± 145.71) in the repertoires of COVID-19 convalescents (Figure 5C), indicating preferred ACE2-competitive neutralizing antibody responses elicited by SARS-CoV-2. Not only the abundance of precursors of shared ACE2-competitive and non-ACE2-competitive antibody clonotypes is different in naïve repertoires (29.44 ± 5.05 versus 12.14 ± 9.81) (Figure 5C), but their prevalence is also remarkably distinct (11/24 versus 5/24) (Figure 5C), suggesting that the immunodominance of the neutralizing antibody response against SARS-CoV-2 is likely determined by the occurrence of precursors in human naïve B-cell repertoires. Finally, we also found that natural SARS-CoV-2 infection induced convergent SARS-CoV-1 cross-reactive neutralizing antibody clonotypes, with 75.63 ± 37.46 reads per million per sample in COVID-19 convalescents but nearly undetectable in all unexposed individuals (Figure 5D). This observation suggests that cross-reactive mAbs against SARS-CoV-1 and SARS-CoV-2 in human naïve B cell repertoires are relatively rare.

Long-lived ighv3-53 antibody lineages would be recalled following omicron breakthrough infection and evolve cross-reactivity to SARS-CoV-2 variants

We next determined the dynamics of spike-specific antibody lineages in COVID-19 convalescents throughout a year post-infection. A total of 141, 267, 95, and 28 spike-specific antibody clonotypes were identified at visit 1-4, respectively (Figure 6A). Spike-specific IgH lineages appeared at visit 1, peaked at visit 2, and then decreased at visit 3 and visit 4 (Figure 6A). The clonal size of these spike-specific IgH lineages also peaked at visit 2 and significantly decreased at visit 3 and visit 4 (Figure 6B). However, the SHM levels of these spike-specific IgH clonotypes that present at visit 3 and visit 4 were significantly higher than those present at visit 1 and visit 2, indicating continued antibody evolution during convalescence (Figure 6C), in agreement with recent reports [36, 37]. The 28 spike-specific IgH clonotypes that persist in the one-year IgH repertoires of SARS-CoV-2 convalescents are encoded by 13 different germline genes, 14 of which are encoded by IGHV-3-53 and IGHV3-30 (Figure S7, Table S6). Among the 12 neutralizing epitopes, only RBD I-targeting antibody clonotypes (IGHV3-53-encoded cl.43 and cl.44) were observed (Figure 5A, S7).

Figure 6. Long-term B-cell memory post ancestral SARS-CoV-2 infection and their clonal relationship with breakthrough infection-induced spike-specific mAbs. (A) Abundance of spike-specific clonotypes in COVID-19 patients at different time points. (B) Dynamics of the clonal size of spike-specific clonotypes one year after SARS-CoV-2 infection. (C) Dynamics of SHM levels of the spike-specific clonotypes one year after SARS-CoV-2 infection. (D) Circos plot shows clonal relationship of spike-specific mAbs induced by wild-type SARS-CoV-2 infections and Omicron breakthrough infections. The width of the link is determined by the number of mAbs shared between two groups. (E) Comparison of heavy chain V domain (VH) SHM levels between mAbs isolated from wild-type SARS-CoV-2 infected and Omicron breakthrough-infected donors (F) Heatmap shows the 98 spike-specific IgH lineages that clonally related to the mAbs isolated from Omicron breakthrough-infected donors. (G) Population prevalence and IGHV usage of the 58 shared clonotypes that clonally related to mAbs isolated from Omicron breakthrough-infected donors. The Wilcoxon rank-sum test was performed in panels (B), (C) and (E) (* p < 0.05; ** p < 0.01; *** p < 0.001).

To gain further insight into the evolution of spike-specific antibody lineages during breakthrough infection, particularly the persistent lineages, we determined the clonal relationship between ancestral SARS-CoV-2-induced antibody lineages and Omicron BA.1 breakthrough infection-induced antibody lineages. We curated another 799 mAbs isolated from Omicron BA.1 breakthrough-infected donors (Table S7), these mAbs are experimentally validated to be cross-reactive to ancestral and Omicron BA.1 strains [7, 38, 39]. By sequence clustering, we found that 72 mAbs isolated from Omicron breakthrough-infected donors were derived from the same lineages as those induced by ancestral SARS-CoV-2 infection (Figure 6D). Sequence analysis revealed significantly higher SHMs of breakthrough infection-induced mAbs than those induced by ancestral SARS-CoV-2 infection (Figure 6E), suggesting that cross-reactive antibody lineages undergo further affinity maturation during breakthrough infection. Next, the IgH sequences of the 799 mAbs were compared to our IgH repertoire data to determine the population prevalence of their precursors induced by ancestral SARS-CoV-2 infection. A total of 98 mAbs isolated from Omicron breakthrough-infected donors could be identified to have clonally related lineages in the repertoires of ancestral SARS-CoV-2-infected individuals, which were classified into 58 shared clonotypes (Figure 6F and G, Table S8). Among them, IGHV3-53-encoded cl.43 and cl.44 represented the most prevalent clonotypes that related to the cross-reactive mAbs isolated from Omicron breakthrough-infected donors (Figure 6G). The Omicron breakthrough infection-induced mAbs such as ADI-75733, BD56-044, and Omi-18 are clonally related to the highly prevalent and long-lived clonotypes cl.43 or cl.44 (Figure 6F and G, S7). Sequence analysis found that more SHMs accumulated in mAbs isolated from 3-dose vaccinees or convalescents that experienced BA.1 breakthrough infection than in those isolated from primary infected individuals (Figure S8). Convergent somatic mutations involved in RBD engagement are more frequently observed in mAbs isolated from individuals who have experienced multiple antigen exposures (Figure S8). It has been demonstrated that single somatic mutations of F27I, F27 V, F27L, or Y58F contribute to improving RBD binding affinity and neutralizing potency, as well as resistance to variants [25, 27, 40]. Indeed, cl.43- or cl.44-related mAbs Omi-18, BD56-619, etc (with massive SHMs) isolated from convalescents that experienced BA.1 breakthrough infection broadly neutralize variants from the ancestral strain to Omicron BA.5 [7, 39], suggesting a recall of preexisting cross-reactive B-cell memory and further affinity maturation. Particularly, Omi-18 contains F27I, S53A, and Y58F (Figure S8), and structural analysis revealed that these SHMs generate new interactions with the RBD and are capable of adapting K417N substitution (Figure S9), explaining its broad neutralization by K417N bearing Omicron subvariants. However, most cl.43- or cl.44-related mAbs P5A-3A1, P22-1D1, etc (with limited SHM), isolated from ancestral SARS-CoV-2-infected individuals are sensitive to K417 substitution [41, 42]. Thus, our observations indicate that the IGHV3-53-mediated long-term immunity elicited by ancestral SARS-CoV-2 infection or vaccination has the potential to develop cross-reactivity following booster vaccination or breakthrough infection through the accumulation of somatic hypermutations.

Discussion

Successful future vaccine strategies as well as the effective use of monoclonal antibodies in clinical applications require an in-depth understanding of the protective antibody component in COVID-19. Traditional approaches for studying epitope-specific antibody responses can be limited by the availability of conformational epitopes. Here, we address this limitation by combining IR-seq, structural, and bioinformatic analyses. Antigenic anatomy and unsupervised classification defined 12 neutralizing epitopes on the SARS-CoV-2 spike protein. To our knowledge, this is the most comprehensive analysis of neutralizing epitopes on spikes. Previously, several classifications based on limited structural information were proposed [36, 43, 44]. However, these classifications are artificially supervised and mostly focused on RBD. We collected 251 available structures of human spike-specific mAbs and performed an unsupervised clustering algorithm for the whole spike, improving the resolution of the antigenic anatomy of the spike. In agreement with serological observations [45], our epitope-resolved profiling showed that ACE2-competitive epitopes on the RBD dominated SARS-CoV-2-induced neutralizing antibody responses. Similar to the RBD, dominant antibody lineages targeting the NTD-neutralizing supersite were also identified. In addition, we mapped the germline usage information to corresponding epitopes and revealed epitope-preferred IGHVs. Our epitope-resolved profiling found that IGHV3-53/3-66, IGHV1-58, and IGHV1-69 antibody lineages represented the most prevalent and abundant neutralizing components targeting RBD. IGHV1-24 antibody lineages represented the most prevalent and abundant neutralizing components targeting NTD. Selection pressure from neutralizing antibodies, particularly the highly prevalent antibody lineages, is associated with VOC escape. The emergence of immune escape variants appears to be driven by those belonging to classes that are widely distributed in the human population. Recently, we have proposed evidence to support that IGHV1-69-encoded public antibody lineages contribute to the occurrence of L452R substitution in different SARS-CoV-2 variants [35, 46]. As part of full-length S, the S2 subunit harbors a sizeable proportion of epitopes targeted by antibodies in response to SARS-CoV-2 infection or vaccination [47]. Our analysis reveals a strikingly high fraction of the S2-tagerting antibody response shared among COVID-19 convalescents that target the S2 non-neutralizing epitopes, but no obvious sequence convergence was observed for mAbs targeting S2-neutralizing epitopes. Antibodies targeting the S2 neutralizing epitopes likely represent a minor component of circulating neutralizing antibodies and likely have minimal selective immune pressure on SARS-CoV-2 spike evolution. This is explained by the lower frequency of substitution occurring in the S2 domain. Nevertheless, the complex interaction of virus evolution and the virus-specific antibody response requires in-depth investigation in future work.

Moreover, this study also elevates our understanding of the origin and development of the antibody response against SARS-CoV-2 and its variants. We not only identified the expanded spike-specific antibody lineages in COVID-19 convalescents but also their antibody precursors in SARS-CoV-2-naïve individuals. Consistently, there is evidence that naïve human B cells engage the spike of SARS-CoV-2 [48, 49]. The vast diversity of human naïve B-cell repertoires [50] contributes to shape SARS-CoV-2-specific antibody responses. Our analysis reveals the remarkable and rapid generation of spike-specific clones after SARS-CoV-2 infection and significantly increased clonal sharing compared with naïve individuals. Consistent with the observation of the rapid waning of neutralizing antibodies in both COVID-19 convalescents and vaccinees [51, 52], we found that only approximately 10% of spike-specific antibody lineages that presented at early times persisted for at least 12 months. Among them, the prevalent and persistent IGHV3-53 antibody lineages are highly clonally related with the cross-reactive IGHV3-53 mAbs isolated from Omicron BA.1 breakthrough-infected donors. Several studies have reported the recall of preexisting cross-reactive B-cell memory encoded by IGHV3-53 following Delta or Omicron BA.1 breakthrough infection [38, 39, 53], potentially helping to clear the virus and prevent severe illness. Efficient recall of Omicron-reactive IGHV3-53 B-cell memory was also detected after a third dose of the SARS-CoV-2 mRNA vaccine [54]. Although IGHV3-53 mAbs isolated from early pandemic samples usually show reductions or loss of activity on K417N/T- or N501Y-containing VOCs, somatic maturation adapts a subset of this class of antibodies to recover potency during breakthrough infection [55]. Indeed, we observed more somatic mutations in IGHV3-53 mAbs isolated from breakthrough infection donors than those isolated from ancestral SARS-CoV-2 infection. Recently, several studies showed that the affinity maturated IGHV3-53 mAbs, such as P5S-2B10 (from WT convalescents) [56], 10-5B (from inactivated vaccinees) [57], and TH132 (from BA.1 breakthrough infected donors) [58], are able to protect K18-hACE2 mice from the challenge of Omicron BA.1, BA.2, and BA.5, respectively. Collectively, our analysis provides a holistic view of the kinetics of the spike-specific antibody response and identifies which of them is induced by ancestral SARS-CoV-2 infection or vaccination that will create long-term or/and cross-reactive B-cell memory. The identification of spike-specific antibody precursors in the naïve human population and further tracking of their evolution in COVID-19 convalescents drives the understanding of how to mount long-term protective responses.

Our study has several limitations. The first limitation is the small number of severe cases in our cohort. Therefore, we have not included a comparative analysis to distinguish COVID-19 patients with disparate disease severities and outcomes. Another limitation is the lack of light-chain analysis. Secondly, the healthy group only has one visit. Additional visits or sequence data may be needed to improve the detection of rare clonotypes. In addition, the binding and neutralizing capacity of the identified antibody lineages were not determined experimentally. Instead, their function was inferred by similarity to sequences deposited in previously published literature. While this approach has been proven by us or other groups for efficient identification of antigen-specific antibodies [19, 32, 59], bulk sequencing is limited by the availability of heavy and light chain paired sequence data and results in the difficulty in functional verification. Future work should include single-cell antibody repertoire sequencing to allow the generation of monoclonal antibodies for further characterization.

In summary, we established a pipeline based on IR-seq, structural, and bioinformatic analysis that can quantify the abundance and prevalence of antibody lineages targeting certain neutralizing epitopes on the SARS-CoV-2 spike protein. We tracked the kinetics of spike-specific antibody lineages for 12 months post-symptom onset at the repertoire level and established the most comprehensive map of spike-specific antibodies in COVID-19 patients. Epitope-resolved antibody repertoire analysis is a useful adjunct to that of antigen-specific B-cell responses. Our study provides a unique perspective on SARS-CoV-2-induced spike-specific antibody responses, informing the design and evaluation of next-generation COVID-19 vaccines.

Author contributions

Y.Q., Z.Y., H.R., P.W., L.H., G.X., W.Q., H.X., D.W., Z.X., L.B., H.P., C.Z., H.J., and N.X. performed the experiments and analyzed the data. Y.Q., and Z.Y. conceived and wrote the paper. C.L., N.X., L.S., Z.J., and X.X. supervised the study. All authors reviewed and approved the final version of the manuscript.

Competing interests

Author Wenjing Pan and Congli Tang are employed by Nanjing ARP Biotechnology Co. Ltd, Nanjing, China. Author Jian Han is employed by iRepertoire Inc., Huntsville, USA. The other authors declare no competing interests.

Acknowledgments

This work supported by grants from National Natural Science Foundation of China (82201932, 92269201, 61971187, 32170941). China Postdoctoral Science Foundation (2022M710891), State Key Laboratory of Respiratory Disease (SKLRD-Z-202324). Guangzhou Science and Technology Bureau (Jointly funded by Municipal Schools and Colleges, 2060206-202201020428).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability

The data that support this study are available from the corresponding authors upon reasonable request. Raw IgH sequences generated in this study have been deposited at the National Genomics Data Center (https://bigd.big.ac.cn/) under accession number PRJCA003775.

Additional information

Funding

This work was supported by National Natural Science Foundation of China: [Grant Number 32170941].