Massively parallel single-cell B-cell receptor sequencing enables rapid discovery of diverse antigen-reactive antibodies

Obtaining full-length antibody heavy- and light-chain variable regions from individual B cells at scale remains a challenging problem. Here we use high-throughput single-cell B-cell receptor sequencing (scBCR-seq) to obtain accurately paired full-length variable regions in a massively parallel fashion. We sequenced more than 250,000 B cells from rat, mouse and human repertoires to characterize their lineages and expansion. In addition, we immunized rats with chicken ovalbumin and profiled antigen-reactive B cells from lymph nodes of immunized animals. The scBCR-seq data recovered 81% (n = 56/69) of B-cell lineages identified from hybridomas generated from the same set of B cells subjected to scBCR-seq. Importantly, scBCR-seq identified an additional 710 candidate lineages not recovered as hybridomas. We synthesized, expressed and tested 93 clones from the identified lineages and found that 99% (n = 92/93) of the clones were antigen-reactive. Our results establish scBCR-seq as a powerful tool for antibody discovery.

A ntibody diversity is an important feature of the adaptive immune system. B cells produce a diverse array of antibodies by rearranging variable, diversity, and joining immunoglobulin germline gene segments [1][2][3] . Somatic hypermutation (SHM) and class switching add to antibody diversity. A mature antibody consists of two identical heavy chains linked through disulphide bonds and two identical light chains each linked to one of the heavy chains, generating two identical antigen-binding sites formed by the first immunoglobulin domain of each chain pair 2 . The heavy and light chains are encoded in separate gene loci, and each B cell normally expresses a single functional heavy and light chain sequence.
Next-generation sequencing has been applied to understand the diversity of the variable regions of heavy (VH) and light chains (VL) that determine the antigen specificity of antibodies. Until recently, the majority of high-throughput sequencing approaches produced unpaired VH and VL repertoires, as generating paired information requires obtaining data at the individual cell level 4 . Recently, techniques that isolate individual cells in microwell plates or droplets of an emulsion, followed by physical linking of VH and VL regions through overlap extension RT-PCR, have demonstrated the potential for obtaining VH-VL pairing information in a high-throughput manner [5][6][7] . However, these techniques only infer full-length variable region sequences indirectly, and single-cell information is lost during library construction. High-throughput approaches that yield full-length variable regions for individual B cells at scale would enable routine application of large-scale immune repertoire sequencing to antibody discovery and detailed repertoire characterization.
Here we describe the application of high-throughput single-cell sequencing to obtain the VH and VL sequences for antibodies from individual human, rat, and mouse B cells. We developed a bioinformatics framework to analyze the sequence data and identify accurate VH and VL pairing. Further, we show the utility of the technique for antibody discovery by expressing and testing predicted antigen-reactive antibody sequences. We demonstrate the potential of direct sequencing of individual antigen-reactive B cells to rapidly generate a large and diverse panel of antigenspecific antibody variable regions and thus expand immune repertoire sampling and expedite antibody discovery processes.

Results
High-throughput single-cell B-cell receptor sequencing. We analyzed >250,000 individual IgG pos B cells from three human donors and two mice, and IgM neg B cells from two rats using emulsion-based encapsulation, cDNA generation and sequencing. Briefly, we generated 5′ barcoded cDNA from thousands of individual B cells in parallel, and amplified the VH and VL regions using custom primers while retaining the cell barcode ( Fig. 1, Supplementary Figs. 1, 2, and "Methods" section). The 5′ barcoded VH and VL domain-encoding cDNAs were sheared and converted into sequencing-ready libraries by addition of appropriate adapter oligonucleotides ( Supplementary Fig. 1). The library construction method involves 3′ cDNA shearing after amplification to create a set of fragments with variable 3′ end, while retaining the 5′ end for all fragments. This resulted in sequencing reads with constant 5′ sequence and variable 3′ sequence, allowing de novo assembly of full-length VH and VL sequences from short-read data (2 × 150 bp). We devised a computational pipeline for cell detection, de novo contig assembly, variable domain annotation, and pairing of full-length VH and VL sequences (Fig. 2a, Supplementary Fig. 3). Assembled VH and VL sequences were parsed for framework regions and complementarity-determining regions (CDR) and selected for open reading frames encoding the entire variable region. Cells with complete VH and VL domains were filtered by requiring a minimum number of reads (10 and 100 for VH and VL, respectively) and requiring one dominant VH and VL contig (≥80% read support; Fig. 2b, Supplementary Fig. 4). As expected, filtered contigs had highest read coverage at the constant 5′ end and more variable coverage at the 3′ end (Fig. 2c, Supplementary  Fig. 5). The median number of reads supporting filtered contigs was at least 175, 213, and 55 throughout the third CDR of the light chain (CDR-L3), and at least 3, 6, and 11 throughout the third CDR of the heavy chain (CDR-H3) for rat, mouse, and human B-cell repertoires, respectively (Fig. 2c, Supplementary  Fig. 5). Pairing efficiency, defined as the percentage of cells with at least one detected VH and VL, was 48-94%, depending on sample type, and 48-74% of cells with pairing information passed quality filters (Fig. 2d). Variation in pairing and filtering efficiency may be due, in part, to differences in heavy-and light-chain transcript levels, primer efficiency, and primer coverage. Overall, a total of 261,539 sequenced B cells yielded high-quality VH-VL pairing information for 116,115 cells from rat (n = 30,380), mouse (n = 9459) and human (n = 76,276) B-cell repertoires ( VH-VL pairing accuracy. To assess the accuracy of assembled VH and VL sequences and VH-VL pairing, we used hybridomas obtained after immunizing rats with chicken ovalbumin (OVA) as a test sample. We analyzed individual clones by a standard sequencing approach to obtain a reference set of paired VH-VL sequences for 127 unique clones (see "Methods"). The same clones were pooled, grown in culture, and subjected to scBCRseq. The scBCR-seq approach yielded high-quality VH-VL sequence pairs for 1989 cells, representing 120 unique VH-VL sequence pairs. More than 90% of cells (1,801/1,989) showed a perfect match with one of the reference VH-VL pairs and more than 99% (1,972/1,989) matched one of the reference VH-VL pairs when allowing for up to two nucleotide mismatches (Fig. 3a).
Next, we assessed VH-VL pairing accuracy based on scBCRseq data from human, mouse, and rat B-cell repertoires. Cells within the same B-cell lineage are expected to have light chains with concordant light-chain variable germline gene segment (V L ) in the majority of cases 8 . We therefore assumed cells with identical heavy-chain variable germline gene segment (V H ) and CDR-H3 belonged to clonally related B cells, and calculated the percentage of cells with concordant paired V L . Concordance in the subset of lineages with more than one B cell was >98% for human (based on 754-2131 cells in 351-894 lineages per sample), >96% for mouse (445 and 979 cells in 148 and 225 lineages, respectively, per sample), and >93% for rat (209 and 598 cells in 97 and 269 lineages, respectively, per sample) (Fig. 3b). These estimates are likely lower bounds for pairing accuracy, due to coincidental V H and CDR-H3 matches in some B-cell lineages, in particular for short CDR-H3 sequences more prevalent in rodents 9 .
Rat and mouse B-cell repertoires. We profiled IgM neg B-cell repertoires from two nonimmunized rats (n = 30,380) and IgG pos B-cell repertoires from two mice (n = 9459) (Fig. 2d, Supplementary Fig. 6a, b). Data at single-cell resolution allowed us to characterize unique B-cell lineages and quantify their expansion based on the number of cells observed for each lineage. We defined B-cell lineages by grouping cells with identical V H and V L germline gene segments and requiring at least 80% nucleotide identity in the CDR-H3 region. For the two rat samples, only 7% (1,081/15,338) and 3% (431/15,042) of cells belonged to clonally expanded lineages ( Supplementary Fig. 7). For the two mouse samples, 13% (699/5,448) and 34% (1,353/4,011) of cells belonged to clonally expanded lineages, suggesting some level of antigenic stimulation in these mice ( Supplementary Fig. 7). We asked whether the identified B-cell lineages showed preferential usage of particular V H genes, V L genes, or V H -V L gene pairings. We observed that some germline gene segments were consistently used more frequently than others across replicates (Fig. 4). This was particularly noticeable for mouse V H and V L genes (ρ = 0.89, ρ = 0.84, Spearman correlation coefficient) and to a lesser extent for rat V H genes (ρ = 0.52) ( Supplementary Fig. 8). For example, IGHV3-2 was the most commonly used mouse V H gene in both animals, present in 6.2% and 4.4% of lineages, respectively (Fig.  4c, d). After correcting for variation in V H and V L gene usage, some individual V H -V L gene pairings showed higher than expected frequencies across replicates ( Supplementary Fig. 8). For example, both mouse samples showed increased frequencies for lineages with IGHV8-12:IGKV3-7 (n = 9, n = 7) and IGHV3-2: IGKV5-43 (n = 45, n = 16). Higher than expected pairings may be due to immune responses to common antigens. Consistent with this interpretation, IGHV8-12:IGKV3-7 lineages showed evidence for expansion in both animals in 3/9 and 5/7 lineages, respectively.
Human B-cell repertoires. Next, we analyzed human IgG pos Bcell repertoires from three donors, each profiled at two different time points (n = 76,276) (Fig. 2d, Supplementary Fig. 6c). All samples showed evidence for lineage expansion, with 14-26% of cells belonging to expanded lineages ( Supplementary Fig. 7). V H and V L gene usage was highly reproducible between replicates from the same individual (ρ > 0.98, Supplementary Fig. 9). Overall V H and V L gene usage was similar among donors (ρ > 0.65 and ρ > 0.87 for V H and V L , respectively, Supplementary Fig. 9) and in general agreement with the known distribution of germline usage in human repertoires 10,11 . Some donors lacked antibodies for a subset of germline gene segments (e.g., IGHV3-9 in Donor 1 or IGHV4-38-2 in Donor 2, Fig. 5), most likely due to genotype differences in the germline repertoires 11 . Interestingly all samples showed higher than expected pairing frequencies for IGHV3-7: IGKV2-30 (19-62 lineages per sample) (Fig. 5, Supplementary  Fig. 9). To rule out a technical artifact due to the profiling method, we reanalyzed published V H -V L pairing information for naive and antigen-experienced human B cells that were obtained by overlap extension RT-PCR and independent computational methods 12 . Interestingly, the published data showed strongest enrichment for IGHV3-7:IGKV2-30 among all variable germline gene segment pairings for antigen-experienced, but not naive B cells, suggesting this pairing may be the result of stereotypical immune responses ( Supplementary Fig. 10).
Rapid discovery of antigen-reactive antibody candidates. To assess the potential of scBCR-seq for antibody discovery, we immunized rats with chicken OVA and subjected IgM neg /OVA pos lymph node B cells from three immunized animals to scBCR-seq ( Supplementary Figs. 6d and 11). After quality filtering we obtained    (Fig. 6b). Chain pairing accuracy assessed by light-chain germline concordance was 99%, consistent with results obtained for naive rats. Somatic mutation load in the V H and V L -derived regions (i.e., excluding regions containing CDR-H3, and the joining gene segments in both chains) was higher in anti-OVA cells than in IgM neg B cells from naive rats (Fig. 6c). In addition to directly sequencing B cells from OVA-immunized animals, we also generated and sequenced OVA-specific hybridomas derived from a fraction of the IgM neg B cells from the same rats. In this dataset we identified 69 unique B-cell lineages, 56 of which were shared with those identified by direct B-cell scBCR-seq (  OVA pos B-cell scBCR-seq data for functional analysis. We selected the BCR sequence pairs from a range of lineage sizes observed in the scBCR-seq data (Supplementary Data 13). A random B-cell clone was selected from lineages represented by three or more individual B cells. For lineages with 1-2 cells, which may be enriched for antigen-negative clones due to imperfect cell sorting, we prioritized clones with high read count based on the assumption that these may be antibody-secreting cells induced by an active immune response. Of the 96 selected B cell receptor (BCR) sequence pairs, 77 were from lineages with three or more cells and 19 were from lineages with 1-2 cells. Of the 96 selected clones, 93 were from lineages not found in the hybridoma panel and 3 had nonidentical but clonally related OVA-specific monoclonal antibodies. Synthesized VH and VL DNA sequences were cloned into human IgG1/κ expression vectors and rat/ human chimeric IgG was expressed in mammalian cells (see "Methods" section). We found that 93 of the 96 clones showed robust protein expression with a median yield of 85 µg of purified IgG per ml of culture (Supplementary Data 13). All but 4 of the 93 clones (96%) specifically bound OVA in an enzyme-linked immunosorbent assay (ELISA) ( Supplementary Fig. 12a). To assess the importance of correct chain pairing for antigen binding in this panel, we shuffled each of the heavy and light chain expression constructs with a noncognate chain to generate 96 additional antibodies with nonnative light and heavy chain pairs. Eighty-six of the 96 shuffled clones expressed as IgG and only five of these showed OVA reactivity ( Supplementary Fig. 12b). While some level of chain pairing promiscuity is expected 13 , the small number of shuffled clones with OVA reactivity confirms the relevance of accurate chain pairing for identification of antigenbinding antibodies. We further determined the affinity of the 93 expressed antibodies for OVA by surface plasmon resonance (SPR). Eighty-nine clones, including three clones with weak binding in ELISA, bound to OVA in SPR, with monovalent equilibrium dissociation constants (K D ) ranging from the limit of detection (10 pM) to 30 nM (geometric mean K D , 2.4 nM; median K D , 3.5 nM) (Fig. 6e, Supplementary Fig. 12a, Supplementary Data 13). The observed affinities were similar to those of antibody panels with comparable size 14 . Interestingly, binding affinity did not correlate with read count, lineage size, or SHM load (Supplementary Fig. 13, Supplementary Data 13). However, as might be expected, B-cell clones with lower SHM loads (≤17 nt mutations in VH and VL combined) showed lower binding affinities (K D ≥ 1.9 nM) in the tested panel ( Supplementary Fig. 13d).
Overall, 86 of the 93 expressed antibodies (92%) bound OVA in both ELISA and SPR, and 92 antibodies (99%) bound OVA in one of the two assays. These findings show scBCR-seq to be a robust tool for rapid discovery of a large panel of antigen-reactive antibodies.

Discussion
Antibody discovery methods based on molecular cloning of antibody repertoires are limited by the ability to process large numbers of single cells in parallel. A number of next-generation sequencing approaches that capture paired antibody heavy-and light-chain information have been described. Techniques that employ manual sorting of cells into microwell plates, followed by deep single-cell RNA sequencing can yield high-quality fulllength VH-VL sequence information, but are inherently lowthroughput [15][16][17]  Low read coverage at the variable region 3′ end can result in partially truncated FR4 sequences, in particular for VH. However, unlike CDR regions that are directly involved in antigen binding, short missing FR4 regions can be reconstructed from germline sequences with minimal impact on antigen-binding properties. In future studies, it may be possible to improve read coverage through optimization of the experimental protocol. For example, reanalysis of a dataset generated by 10x Genomics, using a protocol with reduced fragmentation time, showed higher coverage in the third CDR region (Supplementary Fig. 14).
The scBCR-seq approach described here was performed using commercially available instruments and reagents (10x Genomics) with custom reagents limited to species-specific primers. Importantly, scBCR-seq does not require complex primer sets that cover the highly diverse V H and V L germline gene segments of human and rodent repertoires, minimizing biases and blind spots introduced in the amplification process. Thus scBCR-seq should be applicable with minimal modification to other species, including those for which only constant region sequences are known. Here we used custom primers for rat and mouse, and 10x Genomics primers for human. More recently, 10x Genomics released mouse primers, allowing a direct comparison between custom and commercial primers. Reanalysis of a mouse dataset generated by 10x Genomics (Supplementary Fig. 14) showed largely similar characteristics as mouse datasets generated with custom primers ( Supplementary Figs. 3-5).
Here we only tested a subset of the diversity yielded by scBCRseq. Clones were selected for validation to assess the overall sequence quality in the panel, including pairing reliability. Within each unique lineage, clones were selected randomly. Still, 30% of the selected clones bound with monovalent K D of 1 nM or better.
Although not attempted here, sequence information, including lineage size, SHM load, and read counts per cell, as possible proxies for memory B cells and plasmablasts, or a combination of these, may be useful in selecting clones with higher affinities. In our dataset we did not find clear correlations between these factors and clone affinity. In the case of SHM, similar findings have been previously reported 22 . However, clones tested here were from different lineages, which may have differing ranges of affinities imposed by epitope chemistry and structural properties. Therefore, SHM load may not necessarily correlate with affinities across lineages. Establishing the predictive value of each of these factors for selection of higher affinity clones requires additional testing. However, even in the absence of rational methods for selecting the highest affinity clones, we expect that, once antibodies with the desired epitope specificity are identified, followup analyses of clones within expanded antigen-specific B-cell lineages in the scBCR-seq panels will allow easy and rapid empirical identification of higher affinity variants 23 . The availability of datasets with thousands of paired-chain, full-length variable region sequences from individual antigen-specific B cells will allow detailed large-scale characterization of the dynamics of immune responses, and the functional impact of germline usage, clonal expansion, and somatic mutation.
Contig assembly. We trimmed the first 39 bases for read 1 covering the 16 nt cell barcode, 10 nt unique molecular identifier (UMI) and 13 nt switch oligo. Barcode and UMI sequences were retained for each read. Reads were demultiplexed based on perfect matches to one of the expected barcode sequences. Subsequent processing was done independently for each barcode. If reads for a given barcode exceeded 100,000, they were downsampled to 100,000. Reads were used as input for de novo assembly with SSAKE (3.8.5) 25 with options '-p 1 -c 1 -w 1 -e 1.5' and expected mean insert size 600. We defined concordant read pairs as those with both reads fully embedded in the assembly and paired in the expected orientation. Contigs without concordant read pairs were discarded. Contigs were trimmed to regions supported by concordant read pairs. In cases with more than one contiguous supported region, the region with highest number of concordant read pairs was selected.
Contig annotation. Sequences were annotated with a custom bioinformatics pipeline for variable domain analysis (https://github.com/Genentech/Absolve). Sequences were aligned using HMMER (http://www.hmmer.org/) to Hidden Markov Models that were trained on heavy and light germline amino acid sequences from the International Immunogenetics (IMGT) database 26 in order to determine chain identity, framework and CDR boundaries and residue Kabat numbering 27 . Sequence germline assignments were determined by aligning to IMGT heavy-and light-chain variable and joining germline gene segment sequences 28 using the SSW library 29 and selecting the highest scoring germline. In order to assess the fidelity of germline classifications, Absolve was benchmarked with a set of 2000 simulated VH sequences, including random combinations of 252 variable, 44 diversity and 13 joining human germline gene segments and alleles and random trimming and addition of nucleotides. Each sequence was mutated up to 40 times, preferentially in CDR regions, to yield a set of 80,000 in-frame VH sequences with a variable number of mutations without stop codons. Processing of the simulated dataset with Absolve and IgBlast 30 yielded comparable V H and J H germline call error rates (<0.5% in V H at the highest mutational load). Reference germlines were from the IGMT database except for 52 additional rat germlines mined from Rattus norvegicus whole genome contigs from GenBank and germline sequences deduced from Sprague Dawley rat repertoire sequencing data (I. Hötzel, unpublished) (Supplementary Information).
Heavy-and light-chain pairing. We considered contigs identified as either heavyor light-chain variable domains with HMM score ≥30. Identified variable domains were considered complete if they contained all of FR1 and the first four positions of FR4. The latter requirement ensures that a minimum of four codons match germline sequences without frameshifts immediately following the third CDR in both chains, providing enough sequence information to determine CDR-H3 and CDR-L3 boundaries with high confidence. Amino acid sequences downstream of FR4 position 4 were ignored. For each cell we reported the complete heavy-and light-chain variable domain with highest number of concordant read pairs. For each reported heavy-and light-chain variable domain we calculated a certainty score, defined as the number of read pairs supporting the top contig divided by the total number of read pairs for all identified heavy-and light-chain contigs, respectively. Filtered high-quality pairings are those with (1) top heavy chain supported by ≥10 read pairs (2) top light chain supported by ≥100 read pairs, and (3) certainty ≥80% for both heavy and light chain.
Lineage definition. Filtered cells were grouped into lineage clusters. Two cells were grouped together if they had identical germline V H and V L genes, ignoring allele number, identical CDR-H3 length and ≥80% identical CDR-H3 sequence. J H and J L genes were ignored for lineage definition due to the lower sequence coverage in this region.
Antibody expression and binding assays. Selected B-cell clones from the OVA immunization dataset were produced by DNA synthesis and cloned in mammalian expression vectors 31 as chimeric human IgG1/kappa antibodies. Sequences with clearly incorrect FR4 terminal sequences were corrected based on the assigned joining germline gene segment information for that clone but no additional coding changes were introduced in sequences. Antibodies were transfected into Expi293F TM cells (Thermo Fisher Scientific, Cat. A14635) at 1 ml scale and purified in batch mode by protein A chromatography 32,33 . Antibodies with mispaired chains were made by combining each of the 96 heavy chain clones with a light chain clone from a different anti-OVA antibody. Purified IgGs were tested for binding to OVA in ELISA and in an SPR assay using a BIAcore T200 apparatus (GE Life Sciences, Piscataway, NJ) in a protein A capture format. Purified IgG was immobilized by binding to protein A on chips and soluble monomeric OVA was used as analyte in SPR, using the single cycle kinetics method at 25°C.
Statistics and reproducibility. Replicates refer to independent biological samples. Rat and mouse repertoire datasets each include two samples from two distinct animals. The human repertoire dataset includes samples from three different donors, with two samples obtained at different time points for each donor, six samples in total. The antigen-positive dataset includes one sample pooled from three animals. Reproducibility was verified using replicate samples and publicly available datasets, as indicated. The number of single-cell libraries generated for each sample are listed in Supplementary Data 1.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Rat and mouse datasets generated during the current study are available at the NCBI Sequence Read Archive under project number PRJNA544118. Human datasets generated during the current study have been deposited at the European Genome-phenome Archive under accession code EGAS00001003663. Source data underlying Figs. 2-6 are included in Supplementary Data 1-15.