Recent advancements in microfluidics and high-throughput sequencing technologies have enabled recovery of paired H and L chains of Igs and VDJ and VJ chains of TCRs from thousands of single cells simultaneously in humans and mice. Despite rhesus macaques being one of the most well-studied model organisms for the human adaptive immune response, high-throughput single-cell immune repertoire sequencing assays are not yet available due to the complexity of these polyclonal receptors. We used custom primers that capture all known rhesus macaque Ig and TCR isotypes and chains that are fully compatible with a commercial solution for single-cell immune repertoire profiling. Using these rhesus-specific assays, we sequenced Ig and TCR repertoires in >60,000 cells from cryopreserved rhesus PBMCs, splenocytes, and FACS-sorted B and T cells. We were able to recover every Ig isotype and TCR chain, measure clonal expansion in proliferating T cells, and pair Ig and TCR repertoires with gene expression profiles of the same single cells. Our results establish the ability to perform high-throughput immune repertoire analysis in rhesus macaques at the single-cell level.

This article is featured in Top Reads, p.

Immunoglobulin and TCR repertoire analysis plays a key role in understanding the development of host immunity. These receptor molecules are responsible for recognizing a myriad of foreign Ags from infectious agents. The ability of T and B lymphocytes to give rise to such a diversity of receptor molecules with affinity to these potential Ags is, in part, because of their generation and structure. Igs are tetrameric proteins typically composed of two identical L chains (IgL or IgK) and two identical H chains (IgH) (1). TCRs are heterodimeric proteins composed of paired β (TCRβ) and α (TCRα) or γ (TCRγ) and δ (TCRδ) chains, respectively (2). The IgH chain, as well as TRB and TRD chains, consist of V, D, J, and C region gene segments. Ig L chains, as well as TRA and TRG chains, do not possess D gene segments. Germline V(D)J gene segments exist at large loci within the genome and are somatically rearranged to produce functional and diversified mRNA transcripts and proteins. It is estimated that, within a single individual, the number of possible TCR and Ig V region domains is on the order of 1013 and 1017, respectfully (3).

Increasingly, high-throughput single-cell–based sequencing techniques are being used to profile Ig and TCR repertoires. Single-cell immune repertoire sequencing (scIRS) is a promising new sequencing approach that allows for paired V(D)J repertoire analysis of thousands of cells simultaneously. For example, it has been used to identify multiple neutralizing Abs against SARS-Cov-2 infection in humans (4, 5). However, scIRS assays are species specific, and the development of scIRS assays relies on the complete Ig and TCR reference sequences for the species of interest. Ig and TCR loci are characterized by high levels of repetitive sequences and allelic variation, making targeted sequencing and assembly difficult technical challenges (6). To our knowledge, current commercial scIRS assays are available only for human and mouse.

Due to the close phylogenetic relationship and highly similar physiology to humans, rhesus macaques (Macaca mulatta) have been one of the most popular and well-studied nonhuman primates for modeling immune responses in humans (7, 8). For example, rhesus macaques have been used to model the adaptive immune response and progression of infectious diseases from such agents as varicella zoster (9), HIV (1012), and SARS-Cov-2 (13), as well as many other immune-related studies and diseases, such as allograft rejection (14) and graft-versus-host disease (15). Recently, using long-read transcriptome sequencing, we generated the first complete reference set of C regions of all known isotypes and chain types of rhesus Ig and TCR repertoires (16). We also designed in silico rhesus-specific scIRS assays that remove the need for primers conventionally targeting V regions.

In this study, we sought to experimentally validate and optimize rhesus-specific scIRS assays that are fully compatible with commercial solutions for single-cell immune repertoire profiling. Based on the complete rhesus macaque C region reference set of Ig and TCR isotypes and chains (16), we designed and validated primers that target these C regions in mRNA transcripts. We further adopted these rhesus-specific primers into the human single-cell immune profiling workflow provided by 10x Genomics. These rhesus-specific scIRS assays were validated using cryopreserved PBMCs and splenocytes, as well as FACS-sorted B and T cells from various rhesus animals. We were able to recover every known Ig and TCR isotype and pair Ig/TCR repertoire analysis with transcriptome profiles from the same single cells. We also observed clonal expansion in proliferating (P) versus nonproliferating (NP) rhesus T cells. These results establish the ability to perform high-throughput scIRS analysis in rhesus macaques with comparable performances to commercially available platforms.

PBMCs were obtained from peripheral blood of rhesus macaques using CPT tubes with sodium citrate (BD Biosciences, San Jose, CA). Cells from spleens of rhesus macaques were obtained via manual maceration of the spleen and lysis of RBCs using high-yield lysing solution (Life Technologies, Carlsbad, CA). Cells were cryopreserved in 10% dimethyl sulfoxide (Sigma-Aldrich, Raleigh, NC) and 90% FBS (Corning, Corning, NY) using temperature-controlled freezing containers. To prepare the MLRs, we obtained cryopreserved cells, thawed them, and counted them. CD14+ monocytes were isolated from splenocytes using magnetically labeled Abs (Miltenyi, Auburn, CA). These cells were cultured in R10 Media (consisting of RPMI 1640 supplemented with penicillin-streptomycin at 100 U/ml, l-glutamine at 2 mM [all Life Technologies, Gaithersburg, MD] and 10% FCS) with IL-4 (Miltenyi Biotech, Auburn, CA) and GM-CSF (Miltenyi Biotech) for 7 total days, with the addition of TNF-α (Novus Biologicals, Centennial, CO) on the sixth day to induce activated dendritic cells (DCs). PBMCs were thawed and labeled using Violet Proliferation Dye-450 (BD Biosciences, San Jose, CA) and combined with DCs at a 1:4 ratio. These cells were cocultured in R10 for 5 d. At this time, the culture was stained with CD3-allophycocyanin-Cy7 (557757; BD Biosciences, San Jose, CA) and sent for cell sorting on a BD DIVA (BD Biosciences, San Jose, CA). Our gating strategy is shown in Supplemental Fig. 1.

Multiple sample types from Indian-origin rhesus macaques were used to optimize and validate rhesus-specific scIRS assays: cryopreserved PBMCs and splenocytes, FACS-sorted stimulated B cells, and FACS-sorted stimulated NP and P T cells. Cells were washed in resuspension buffer (RPMI, 10% FBS), collected by centrifugation at 700 × g, and resuspended at an appropriate concentration, between 700 and 1200 cells/μl. All cells were counted by hemocytometer and assayed for viability using a Countess II automated cell counter (Thermo Fisher Scientific, Waltham, MA). From each sample, a maximum volume containing 17,000 cells, having ∼90% of viable cells, were loaded onto the 10× Chromium (10x Genomics, Pleasanton, CA) controller for a targeted cell recovery of 10,000 cells (Chromium Next GEM Single Cell V(D)J Reagents Kits v1.1_RevE; 10x Genomics). Two replicates of the same PBMC sample were established (PMBC1 and PMBC2), whereas only one sample was prepared for each of the remaining cell types. cDNA amplification was carried out at 13 cycles of amplification for all samples and assayed for quality and concentration using a Bioanalyzer 2100 and High Sensitivity DNA Kit (Agilent, Santa Clara, CA).

Construction of V(D)J libraries required the implementation of rhesus-specific primers (16) for two subsequent enrichments of Ig and TCR transcripts as is typical for the 10× V(D)J workflow (Supplemental Table I). Quantitative PCR results of various primer pools, which varied primer concentrations to enrich for lower abundant transcripts, showed no impact on final sequencing data. Therefore, an equimolar ratio of gene-specific primers was used to construct each pool, whereby Ig assays used a final concentration of 0.5 μM of each gene-specific primer paired with 1 μM 10× forward primer, and TCR assays used a final concentration of 1 μM of each gene-specific primer paired with 2 μM of the 10× forward primer. Due to initial sequencing results, later Ig and TCR primer pools were split into VDJ and VJ pools to improve transcript capture and chain pairing efficiency. Primer pools were constructed in volumes of 10 μl, requiring the nuclease-free water volume added to the cDNA at the beginning of target enrichment to be reduced to 28 μl from 33 μl. Other than the needed changes to accommodate novel primers to the 10× VDJ workflow, VDJ and 5′ GEX libraries were prepared according to the manufacturer’s instructions (Chromium Next GEM Single Cell VDJ Reagents Kits v1.1; 10x Genomics).

Final V(D)J and 5′ GEX libraries were run on a Bioanalyzer 2100 with a High Sensitivity DNA kit (Agilent) to assess library size and concentration. Further analysis of library concentration was performed using a Qubit 3 fluorometer (Thermo Fisher Scientific) coupled with the sizing data from the Bioanalyzer to determine appropriate loading concentration for each library. Sequencing was performed using a NextSeq 500 sequencer (Illumina, San Diego, CA) using 150 cycle kits in a paired-end fashion. Runs were programmed to generate a 26-bp read 1 sequence consisting of the 16-bp 10× barcode and 10-bp 10× unique molecular identifier (UMI), an 8-bp index read, and a 133-bp read 2 sequence of the cDNA insert to exhaust the number of cycles available in the kit. Sequencing depth was targeted at a minimum of 5000 reads per cell for V(D)J libraries and 20,000 reads per cell for 5′ GEX libraries as recommended by the 10× workflow (Chromium Next GEM Single Cell VDJ Reagents Kits v1.1; 10x Genomics).

Raw sequencing data were demultiplexed and converted to fastq files also using the Cellranger (v4.0.0) command mkfastq. Next, raw reads were de novo assembled into contigs using the Cellranger vdj command. A de novo, rather than reference-based, contig assembly was performed because of incomplete annotation of the rhesus germline V(D)J sequences. To annotate the C regions of the assembled V(D)J transcripts, we searched assembled contigs against inner-enrichment primers and a reference of C region amplification sequences using the ublast (e value ≤ 1e−4) and usearch_global (sequence ID threshold = 0.9) commands, respectively, from the USEARCH (v10.0.240) (17) suite of sequence analysis tools. Subsequently, to annotate the V regions of the assembled V(D)J transcripts, we aligned the assembled contigs to a custom IgBLAST (v1.8.0) (18) database of human and rhesus germline V(D)J sequences (available at https://github.com/ncsu-penglab/RhesusIgTCR). An e value cutoff of 0.01 and default parameters were used in IgBLAST queries. After annotation, we systematically filtered assembled contigs from our analysis as follows: contigs representing misassembled chimeras (e.g., IgL primer alignment and IgM C region alignment) or off-target transcripts (e.g., C region primer match at the 3′ end but without a corresponding C region alignment) were removed. Further, assembled contigs without a productive and complete V(D)J sequence or that had a premature stop codon were removed. Only contigs that passed all of these quality-control (QC) criteria were considered in downstream analyses, including the assessment of pairing efficiency.

We assessed read coverage of filtered contigs using the all_contig.bam file outputted by cellranger vdj, where sequencing reads were aligned to contigs after assembly. Reads were only aligned to contigs of their respective cell barcode. For each dataset, the samtools depth command was run on the bam file to generate the read coverage per contig base position. To facilitate the comparison of read coverages across the contigs of different lengths, we normalized the read coverages as follows: for each assembled contig, we first labeled each base position as part of the 5′ untranslated region, V region, or C region, based on the start positions of the V gene segments and the end positions of the J gene segments from IgBLAST results. Base positions within each labeled region were placed into 100 bins of equal length, and the mean relative coverage per bin was calculated.

VDJ versus VJ chain pairing efficiency was assessed using B and T cells defined by the 5′ GEX data collected from the same cell samples described later. After removing low-quality cells and identifying the B and T cell clusters, the remaining cell barcodes were used as a ground truth to assess chain pairing efficiency in V(D)J sequence data because the same cell barcodes were already independently validated as a productive B or T lymphocyte based on the filtered V(D)J transcript contigs. The V(D)J sequences of each cell were integrated into a Seurat object as metadata for gene expression and clonotype analysis.

Expanded clones were defined by more than one cell sharing the same V germline gene segments and with at least 85% identical CDR3 nucleotides sequences, for both VDJ and VJ contigs. Cells that had unpaired contigs were not considered in the clonotyping analysis. CDR3 sequence clustering was performed using the CD-HIT (v4.6.6) (19) sequence analysis tool.

Raw sequencing data were demultiplexed, aligned to the RheMac10 genome annotation, and UMI-collapsed using the 10× Genomics cellranger (v4.0.0) commands mkfastq and count, respectively. Raw gene expression matrices were normalized and scaled using the SCTransform (20) method with the Seurat R package (v3.1.5). QC was performed on each dataset to remove poor-quality cells. For each sample, cells that had >5% of their UMIs map to mitochondrial genes were removed from the analysis.

Principal-component analysis was performed using the normalized and scaled expression levels of the 3000 most variable genes in each dataset. Based on Seurat’s recommendations, the first 30 principal components were used as input to uniform manifold approximation and projection (UMAP) dimension reduction and K-nearest neighbor cell clustering. Canonical marker gene expression and differential expression testing were used to determine the cell types present in the tissue culture samples. Differential expression analysis was performed using a Wilcoxon rank sum test within the FindMarkers Seurat function.

The potential for CD3+ B cells to be technical B+T cell doublets was assessed by simulating B+T cell doublets by combining UMI counts from randomly selected pairs of B and T cells. The counts for gene j in doublet i with parent cells a and b was defined as yi,j = xa,j + xb,j. The new simulated doublets were added to the unnormalized gene expression matrix, whereas the original pair of singlets was removed. Next, the standard Seurat analysis workflow was performed as described earlier. Once cells were embedded with individual UMAP coordinates, the Euclidian distance between CD3+ B cells and simulated doublets was measured. This analysis was performed using between 10 and 100 simulated doublets by an increase of 10 cells. The simulation and distance calculations were repeated 10 times for each number of simulated doublets. Wilcoxon rank sum tests were performed to test for a difference in distance between CD3+ B cells to each other and to simulated B+T cell doublets.

As shown in (Fig. 1, we adopted rhesus-specific V(D)J primers into the human single-cell immune profiling assays commercially available from 10× Genomics and sequenced in total >80,000 single-cell barcodes from cryopreserved rhesus PBMCs, splenocytes, FACS-sorted memory B cells (CD20+CD27+), and stimulated T cells (CD3+; sorted into P and NP fractions by dye dilution). De novo assembly of raw sequencing reads resulted in >150,000 unfiltered, unannotated rhesus V(D)J transcript contigs. We devised a custom computational pipeline for contig isotype and germline annotation, as well as QC filtering standards (see Materials and Methods). Assembled VDJ and VJ transcript sequences were searched against a custom IgBLAST database of rhesus and human germline segments and selected for sequences encoding complete, productive variable domains. Filtering using these criteria yielded the final set of 63,623 (78.8%) cell barcodes with a total of 97,932 (61.8%) V(D)J sequences from six rhesus samples (two PBMCs, splenocytes, sorted memory B cells, and P and NP T cells) (Supplemental Table I).

FIGURE 1.

Overview of assay development. (A) Droplet-based single-cell RNA sequencing was used to generate the initial barcoded cDNA libraries. To maximize the rate at which cells with paired VDJ and VJ contigs were recovered, we generated sequencing libraries using pooled and separate primers for VDJ and VJ chain types. After sequencing and data analysis, the chain pairing rate was assessed and compared between the two configurations using filtered B and T cell barcodes recovered from 5′ gene expression analysis. (B) Targeted enrichment performed on each cDNA sample. The PBMC1 sample is labeled with an asterisk (*) to indicate that a gene expression library was not constructed for this dataset, so a comparison between split and pooled library preps was not performed. (C) Bar charts indicate the percentage of cell barcodes with either no V(D)J contigs that pass QC filtering, an unpaired contig(s) that does pass, or paired contigs that do pass. In each sample, the fraction of paired contigs that pass QC increases when VDJ and VJ primers are split in separate library preps.

FIGURE 1.

Overview of assay development. (A) Droplet-based single-cell RNA sequencing was used to generate the initial barcoded cDNA libraries. To maximize the rate at which cells with paired VDJ and VJ contigs were recovered, we generated sequencing libraries using pooled and separate primers for VDJ and VJ chain types. After sequencing and data analysis, the chain pairing rate was assessed and compared between the two configurations using filtered B and T cell barcodes recovered from 5′ gene expression analysis. (B) Targeted enrichment performed on each cDNA sample. The PBMC1 sample is labeled with an asterisk (*) to indicate that a gene expression library was not constructed for this dataset, so a comparison between split and pooled library preps was not performed. (C) Bar charts indicate the percentage of cell barcodes with either no V(D)J contigs that pass QC filtering, an unpaired contig(s) that does pass, or paired contigs that do pass. In each sample, the fraction of paired contigs that pass QC increases when VDJ and VJ primers are split in separate library preps.

Close modal

We observed normally distributed read coverages across the 5′ untranslated regions and V regions of the V(D)J contigs (Supplemental Fig. 2A). C region reads were uniformly distributed and did not contain bias as a result of the 10× capture method. We reason that our lack of bias is the product of an increased read 2 length (133 bp) being of similar length to the 5′ portion of C regions (78–108 bp) captured by our assay (Supplemental Fig. 2B). Therefore, those sequencing reads corresponding to the unfragmented cDNA molecules would cover the entire amplified C region. The average numbers of reads supporting the final filtered contigs in Ig enrichment libraries were between 117 and 2434, and the numbers of UMIs were between 1 and 9. The median numbers of reads and UMIs supporting the final filtered contigs in TCR enrichment libraries were comparable, ranging from 763 to 2282 for reads and from 1 to 4 for UMIs (Supplemental Table I).

To improve the recovery of paired VDJ and VJ sequences from the same cells (e.g., TRB and TRA or IgH and IgL chains), we devised a sequencing library construction strategy in which VDJ and VJ transcripts were enriched in both separate and pooled sequencing libraries (Fig. 1A, 1B). Pairing efficiency, defined as the percentage of cells with at least one final filtered VDJ and one final filtered VJ sequence, may be affected by differences in VDJ and VJ chain transcript abundances, as well as individual primer amplification efficiencies. To facilitate the comparisons of pairing efficiency between VDJ and VJ split and pooled libraries, we additionally sequenced and generated gene expression profiles of the same single cells from five rhesus samples. B and T cell barcodes that passed gene expression QC filtering served as an independent “ground-truth” of viable B and T cells for assessing pairing efficiency. When VDJ and VJ primers were used for separate library preparation, an increase in the fraction of cells with paired VDJ-VJ repertoires was observed across all five samples and both Ig and TCR enrichments (Fig. 1C). Between 60 and 74% of filtered B cell barcodes from gene expression libraries had paired H and L chains (Supplemental Table I). We observed between 30 and 68% of T cell barcodes that had paired VDJ and VJ sequences (Supplemental Table I). Both Ig and TCR pairing efficiency ranges were comparable with other scIRS assays in the split library configurations (21, 22). We noted a relatively lower pairing efficiency rate observed in T cells from the one splenocyte sample. Overall, we recovered a total of 7539 B cells and 9906 T cells with paired VDJ and VJ contigs, as well as a filtered gene expression profile (Supplemental Table I).

We analyzed Ig repertoires from two PBMC replicate samples (the total number of B cells sequenced: n = 3418 and n = 3035), one splenocyte sample (n = 15,682), and one sorted B cell sample (n = 9746) (Supplemental Table I). The number of cell barcodes with at least one filtered V(D)J sequence will be higher in the absence of paired GEX data because of the inability to filter out doublets, nonviable cells, or ambient transcripts assigned to cell barcodes. Overall, we detected every known Ig isotype (Fig. 2A) and were also able to recover every known IgA allotype and IgG subclass identified in our previous full-length isoform sequencing analysis, observing preferential usage of IgA*01 and IgA*02, as well as IgG1, respectively (16). No IgE contigs were detected in the PBMC2 and splenocyte samples; however, IgE is known to be detected at the lowest levels among any IgH isotype (22). Ig L chains (IgK and IgL) were more frequently detected than IgH contigs across all four samples. We also investigated the ratio of cells with IgA, IgE, and IgG contigs to those with IgM and IgD. IgM and IgD isotypes are known to be expressed in mature naive B cells, while IgA, IgE, and IgG are expressed by activated (memory) B cells after undergoing an IgH class-switch recombination process (3). Overall, this IgH ratio ranged between 0.39 and 0.60 for our four Ig datasets, suggesting that we might have captured about twice as many naive relative to activated B cells, on average (Supplemental Fig. 2D). As expected, the largest ratio of activated to naive B cells appeared in the FACS B cell dataset, which were sorted using cell surface markers of B cell memory (CD20 and CD27) (23).

FIGURE 2.

Characterization of isotypes Ig enrichment libraries. (A) Recovery of 3′ of enriched V(D)J transcripts. Tiled plots indicate the number of contigs with respective primer and C region alignments. (B) Distribution of contig length by Ig isotype across the Ig enrichment datasets. IGHE was omitted because of the low number of contigs detected. A different IgM inner-enrichment primer was used for the PBMC1 library construction, which was 59 bp 5′ of the primer used for all other datasets. This is consistent with the differences in the length distributions. (C) Distribution of the distances between C region and primer alignments and 3′ end of contigs.

FIGURE 2.

Characterization of isotypes Ig enrichment libraries. (A) Recovery of 3′ of enriched V(D)J transcripts. Tiled plots indicate the number of contigs with respective primer and C region alignments. (B) Distribution of contig length by Ig isotype across the Ig enrichment datasets. IGHE was omitted because of the low number of contigs detected. A different IgM inner-enrichment primer was used for the PBMC1 library construction, which was 59 bp 5′ of the primer used for all other datasets. This is consistent with the differences in the length distributions. (C) Distribution of the distances between C region and primer alignments and 3′ end of contigs.

Close modal

The obtained Ig V(D)J contigs exhibited a high rate of complete assembly of the C regions at the 3′ end, and 97.5–98.4% mapped to an inner-enrichment primer or expected amplified C region sequence (Fig. 2A). Ig V(D)J contigs that mapped to the inner-enrichment primer or amplified C region sequences also exhibited higher sequencing depth in terms of both reads and UMIs compared with those that did not map in the same fashion (Supplemental Fig. 1C). 3′ C regions of Ig V(D)J sequences mapped to the 3′ end of Ig V(D)J contigs to ensure that contigs were correctly assembled (Fig. 2C). The distribution of isotype contig lengths was also largely consistent across the different datasets, providing confidence that our assembled contigs represent true BCR mRNA molecules (Fig. 2B).

Similarly, we analyzed TCR repertoires from two PBMC replicate samples (the total number of T cells sequenced: n = 3988 and n = 3543), one splenocyte sample (n = 3162), and two samples derived from a MLR experiment representing P (n = 12,964) and NP (n = 8085) FACS-sorted T cells (Supplemental Table I). Like our B cell repertoire analysis, TCR V(D)J contigs were searched against the inner-enrichment primers used to build TCR libraries in addition to a reference of expected amplified C region sequences. Importantly, in addition to the more commonly targeted TRA and TRB transcripts, we also included primers capturing TRD and TRG transcripts. Notably, primers that target TRD and TRG chains are not currently available in the human and mouse single-cell immune profiling assays from 10x Genomics, resulting in repertoire analyses solely covering αβ T cells (24, 25). γδ T cells have been reported to make up only 4% of T cells, on average (26). As expected, we detected TRD and TRG chains at a lower frequency than the predominantly expressed TRA and TRB chains (Fig. 3A). However, we observed TRA+ cells to be only as much as 1.59 times more abundant than TRG+ cells, and TRB+ cells were ∼15.8 times more abundant than TRD+ cells in mixed cell datasets (i.e., PBMCs and splenocytes), suggesting that leaving out TRG and TRD sequences could potentially miss relevant biological insights. Interestingly, overall TCR VDJ contigs were 1.38× more abundant than VJ contigs in αβ T cells, but 4.52× less abundant in γδ T cells (Fig. 3A).

FIGURE 3.

Characterization of isotypes in TCR enrichment libraries. (A) Recovery of 3′ of enriched V(D)J transcripts. Tiled plots indicate the number of contigs with respective primer and C region alignments. (B) Distribution of contig length by TCR isotype across the TCR enrichment datasets. (C) Distribution of the distances between C region and primer alignments and 3′ end of contigs.

FIGURE 3.

Characterization of isotypes in TCR enrichment libraries. (A) Recovery of 3′ of enriched V(D)J transcripts. Tiled plots indicate the number of contigs with respective primer and C region alignments. (B) Distribution of contig length by TCR isotype across the TCR enrichment datasets. (C) Distribution of the distances between C region and primer alignments and 3′ end of contigs.

Close modal

Similar to our Ig repertoire analysis, we also observed a high rate of complete C region assembly in our TCR enrichment datasets with 97.8–98.5% of TCR V(D)J contigs mapping to an inner-enrichment primer or amplified C region sequence (Fig. 3A). As was observed with Ig V(D)J contigs, TCR V(D)J contigs that mapped to the inner-enrichment primer sequences or C region sequences also exhibited higher sequencing depth in terms of reads and UMIs than those that did map accordingly (Supplemental Fig. 2C). Correct C region assembly was confirmed by observing C region amplification sequences and inner-enrichment primers mapped directly to the 3′ end of the TCR V(D)J contigs (Fig. 3C). The distributions of TCR V(D)J contig lengths were also largely consistent across different samples, supporting that our assembled TCR contigs represented the full-length TCR V(D)J transcripts (Fig. 3B).

To assess the potential of our rhesus scIRS assays to detect clonally expanded lineages, we compared the clonal expansion rates in our P and NP T cell datasets. These cells were isolated from an MLR, in which responding PBMCs were labeled with a proliferation dye and incubated with MHC mismatched activated DCs for 5 days using a previously described protocol (Supplemental Fig. 1) (27). We then sorted these cells based on CD3 positivity and whether they had undergone proliferation as indicated by a loss of proliferation dye. P cells were enriched for those cells responding to alloantigens from the MHC mismatched DCs. We defined T cell lineages by grouping cells with identical VDJ and VJ germline variable gene segments and required at least 85% nucleotide identity in both CDR3 regions, although the number of clones did not change significantly when this sequence identity threshold was increased (Fig. 4A).

FIGURE 4.

Detection of clonally expanded lineages in P T cells. (A) Number of cells in expanded lineages in NP and P T cells as a function of the CDR3 sequence identity threshold. (B) Lineage expansions among P sorted T cells. Pie charts indicate percentage of cells belonging to expanded lineages. (C) Same as (B), but for NP T cells. (D) Tiled plots indicate the number and percentage of cells with the respective TCR chain types across the two TCR enrichment datasets.

FIGURE 4.

Detection of clonally expanded lineages in P T cells. (A) Number of cells in expanded lineages in NP and P T cells as a function of the CDR3 sequence identity threshold. (B) Lineage expansions among P sorted T cells. Pie charts indicate percentage of cells belonging to expanded lineages. (C) Same as (B), but for NP T cells. (D) Tiled plots indicate the number and percentage of cells with the respective TCR chain types across the two TCR enrichment datasets.

Close modal

As expected, we observed a substantial clonal expansion in the P T cell dataset with 42.5% of cells belonging to clonally expanded lineages, compared with <2% in the NP T cell dataset (Fig. 4A–C). Despite the high clonal expansion rate in P T cells, we did not observe any lineages that dominated the clonally expanded cells. Out of the 2166 cells of expanded lineages in the P T cell dataset, 530 (24.5%) belonged to a lineage of size two. The largest lineage size consisted of only 24 cells. Hierarchical clustering of TCR VDJ and VJ variable germline gene segments in clonally expanded P T cells did not reveal any modules of highly overlapping VJ and VDJ combinations. MLRs typically result in a generalized activation of some T cells predisposed for proliferation. There are a multitude of potential Ags that can activate T cells, there is a propensity of stimulatory cytokines to induce bystander proliferation in previously primed cells (e.g., those with a memory phenotype), and heterologous cross-reactive stimulation can be expected (28). Nevertheless, the products of an MLR are substantially enriched for allospecific cells.

Interestingly, among P T cells, the most frequent VJ germline gene segment detected in expanded lineages was TRGV10 (Supplemental Fig. 3A). This observation correlated with the codetection of TRG molecules in many T cells with TRB molecules (Fig. 4D). In addition, TRGV10 showed the largest increase in usage between expanded P T cells and NP T cells among any VJ or VDJ germline segment (Supplemental Fig. 3B). This observation was intriguing, given the high correlation of germline usage rate between expanded P and NP T cells (Supplemental Fig. 3C). In some of the clonotypes, no additional TRB contig was detected, but the frequency of αβ T cells that also expressed a TRG chain showed the largest increase from NP to expanded P T cells (Fig. 4D). Similar to what others have reported, we observed a diverse number of TRA and TRB chains and a limited number of TRG chains in our γ-αβ T cells in both NP and P + expanded T cell datasets (29). In both NP and P + expanded γ-αβ T cells, the TRGV10 variable gene was used in >80% of lineages (singleton or expanded). We did not observe one TRAV and TRBV gene used in >12% of lineages.

Some studies have reported T cells that coexpress αβ and γδ receptors (29, 30). Others have characterized δ/αβ T cells that possess chimeric contigs that encompass a δ variable gene but α J and C region domains (31). This chimeric scenario is unlikely to be the case in these data because chimeric contigs as the one described in δ/αβ T cells would be removed from our analyses (see Materials and Methods). The frequency of γ-αβ T cells was between 9.32 and 12.46% of the T cells in PBMC and splenocyte samples, while γ-β T cells were between 11.77 and 25.39% (Supplemental Fig. 3D). It shows that the observed coexpression was not unique to these FACS-sorted T cell datasets. Because typical single-cell–based TCR repertoire sequencing analyses have focused on αβ T cells, it is unclear whether these αβ T cells coexpressing TRG or TRD tend to be undetected or are atypically deleted at some point postproduction. In addition, the presence of TRG transcripts in αβ T cells does not imply γδ TCRs are on the surface of the same T cells. Obviously future analyses of more samples are needed. However, this complete chain pairing coverage highlights the unbiased potential of our rhesus scIRS assays, which could also enable future exploratory analysis of TRA/TRB and TRG/TRD coexpressions.

Parallel 5′ gene expression profiling analysis allowed us to cluster the same single cells into groups of cell types and to overlay individual cells with their Ig and TCR repertoire data (Fig. 5A, 5B). For example, in the PBMC2 sample (Fig. 5A), we identified 2100 B cells across three cell clusters. We recovered 1344 (64.1%) B cells with paired H and L chains, 636 (30.3%) with only an L chain, and 62 (3%) with only a H chain. Among the 3071 T cells in the same PBMC2 sample, we recovered 1495 (48.7%) cells with paired TCR VDJ and VJ chains, 651 (21.2%) cells with only a VDJ, and 582 (19%) cells with only a VJ chain. The rates of paired TCR chains were similar across the different T cell subtypes. Among the 3531 B cells identified in the splenocyte sample (Fig. 5B), we recovered 2118 (60%) B cells with paired H and L chains, 817 (23%) with only an L chain, and 328 (9.3%) with only an H chain. In the same splenocyte sample, we captured 1886 T cells. Of those, we recovered 575 (30.5%) cells with paired TCR VDJ and VJ chains, 472 (25%) cells with only a VDJ, and 347 (18.4%) cells with only a VJ chain. In addition, we observed 492 (26.1%) T cells without any TCR V(D)J contigs (Supplemental Table I).

FIGURE 5.

Integration of gene expression and V(D)J profiles of single cells. (A) UMAP plots of the (left to right) assignment of cell type, Ig chains, TCR chains, and dual Ig and TCR detection for each cell in PBMC 2 sample. (B) Same as (A), but the splenocyte sample. (C) Normalized expression of Ig and TCR signal transduction genes, as well as CD3+ B cell gene expression signatures across CD3 B cells (green) and randomly sampled B (gold) and T cells (blue). (D) Distribution of UMIs and unique features in CD3+ B and B cells.

FIGURE 5.

Integration of gene expression and V(D)J profiles of single cells. (A) UMAP plots of the (left to right) assignment of cell type, Ig chains, TCR chains, and dual Ig and TCR detection for each cell in PBMC 2 sample. (B) Same as (A), but the splenocyte sample. (C) Normalized expression of Ig and TCR signal transduction genes, as well as CD3+ B cell gene expression signatures across CD3 B cells (green) and randomly sampled B (gold) and T cells (blue). (D) Distribution of UMIs and unique features in CD3+ B and B cells.

Close modal

Intriguingly, among the cells sequenced in the PBMC2 dataset, we observed a small cluster of 23 B cells clustered separately from the others (Fig. 5A). On investigation, we obtained both Ig and TCR V(D)J contigs from a high fraction of these cells. Of these 23 cells, Ig and TCR dual expression was detected in 22 (>95%) cells. Dual expression of Ig and TCR molecules was also observed in 229 additional cells, but these cells appeared to be randomly distributed in UMAP space. Therefore, we focused our exploratory Ig/TCR dual-expression analysis on this small cluster of 23 cells.

Dual expression of Ig and TCRs or other B and T cell markers has been reported by others, but the role that these hybrid lymphocytes play is largely unknown (3234). The possibility of doublets, where two cells are given the same cellular barcode, does arise when using droplet-based single-cell technologies, such as the 10x Genomics used in this study. To assess the potential for these cells to be B+T cell doublets, we used differential expression analysis to identify transcriptomic signatures that were absent in B and T cells. This small B cell cluster not only expressed both Ig and TCR molecules but also the genes involved in B and TCR signal transduction. CD79A, CD79B, CD3E, CD3G, and CD3D were all upregulated in this small cluster of B cells. As such, we labeled this cluster as CD3+ B cells. The top differentially expressed genes appeared to be constitutively expressed across these CD3+ B cells but absent in regular B and T cells (Fig. 5C). Furthermore, CD3+ B cells did not have a greater sequencing depth than other B cells in terms of number of UMIs or unique features expressed (Fig. 5D). To further demonstrate the transcriptional differences of CD3+ B cells from technical B+T cell doublets, we simulated B+T cell doublets and added their expression profiles to the gene expression matrix to rerun dimension reduction and clustering (see Materials and Methods). To assess the transcriptional similarity of CD3+ B cells and simulated B+T cell doublets, we measured the Euclidean distance between CD3+ B cells to each other and to simulated doublets. Consistently, CD3+ B cells were closer to each other (3.06 ± 0.397) (i.e., more transcriptionally similar) than to simulated doublets (5.14 ± 0.570). These results were additionally found to be statistically significant by Wilcoxon rank sum test (adjusted p < 0.05). Taken together, these data suggest that these CD3+ B cells were less likely to be B and T cell doublets and more likely to be expressing both Ig and TCR transcripts.

In this study, using extensive experimental optimization and validation, we established rhesus-specific single-cell high-throughput assays for Ig and TCR repertoire sequencing analysis that recover every known rhesus Ig and TCR isotype and chain type. Although rhesus macaque is one of the most widely used nonhuman primate models, it is our understanding that this is the first such assay developed for rhesus macaques that are also fully compatible with a commercially available platform. We also devised a custom computational pipeline to process and analyze the generated rhesus single-cell V(D)J sequencing data.

The establishment of these single-cell–based Ig/TCR repertoire sequencing assays will improve our understanding of immune responses in rhesus macaques, which will in turn allow for applications to further understanding the immune response in humans. For example, we were able to detect clonally expanded lineages of T cells in a P T cell sample, using the full-length productive VDJ and VJ sequences paired within single cells. The overlay of cellular gene expression profiles allowed the identification of the underlying cell types and subsets. We also observed rare cells expressing both Ig and TCR V(D)J transcripts that could be investigated in the future. Compared with the human and mouse assays commercially available from 10× Genomics, our assays also capture TRG and TRD chains. Although αβ T cells are more frequent in general, our results showed that the frequencies of γδ T cells could be extremely variable, and sometimes their frequency could be as high as 26% of V(D)J contigs detected in a PBMC sample. This suggests that the additional coverage of TRG and TRD TCRs may enable the discovery of relevant biological insights that are overlooked by current assays.

Unexpectedly, we detected γ-αβ T cells in all four types of rhesus samples that we analyzed in this study: PBMC, splenocyte, FACS-sorted NP, and clonally expanded P T cell samples. Interestingly, we also observed a large increase in the frequencies of γ-αβ T cells in the clonally expanded P T cell sample compared with the NP T cell sample (Fig. 4D). Generally, a mutually exclusive pattern of TRG/TRD and TRA/TRB expression is expected (35). However, reports have demonstrated that there are exceptions. Using FACS in combination with quantitative PCR and single-cell RNA sequencing, research studies have observed and validated the presence of dual-expressing γδ-αβ T cells at the mRNA and cell surface receptor level in multiple tissues and developmental stages of human and mice (29, 30). Furthermore, transgenic mice constructed with a pair of productively rearranged TRG and TRD genes have been shown to produce a normal number of αβ T cells (36), in addition to other reported unconventional expression patterns of TCR molecules (37, 38). Considering potential technical confounding factors, such as ambient transcripts and technical doublets, future independent analyses will be needed to experimentally investigate the rhesus γ-αβ T cells observed in this study.

Analysis of the adaptive immune response has coevolved with the advancement of high-throughput immune repertoire sequencing (IRS) strategies. For RNA-based analysis, 5′ rapid amplification of cDNA ends has been used to reliably generate unbiased amplification of Ig and TCR cDNA molecules and can be done using a simple primer set (3941). The major limitations of 5′ rapid amplification of cDNA ends and similar strategies have been extracting the pairing of V region information required for the determination of Ag/epitope specificity. Therefore, single-cell high-throughput sequencing strategies started to emerge to recover paired V region sequences. Emulsion PCR, where TCRβ and TCRα RNA molecules from single cells are fused together (42, 43), and pairSeq, where single cells are isolated and mRNA reverse transcription reactions occur individually within a 96-well plate, became commonplace (44). Despite the advancement in chain pairing of these IRS strategies, these methods are generally low throughput in terms of the number of cells sequenced in a single experiment and are not applicable to many Ig and TCR chain types. The recent advent of droplet-based single-cell sequencing technologies has greatly improved the throughput of IRS analysis. The human and mouse immune profiling assays from 10× Genomics allow for the recovery of paired, full-length V(D)J sequences, as well as gene expression profiles of thousands of cells simultaneously.

However, the development of similar scIRS assays for other species is challenging. For a species of interest, it requires the complete Ig and TCR reference sequences, which are often lacking due to the extreme complexity of Ig and TCR loci. Recently, we obtained a complete reference set for the C regions of rhesus Ig and TCR isotypes and chain types using long-read transcriptome sequencing (16), which bypasses the need of rhesus-specific germline reference sequences. This reference set enabled us to design and validate primers targeting rhesus V(D)J C region sequences to achieve complete and unbiased repertoire analysis, whereas primer pools targeting the highly diverse V regions can miss as much as 50% of repertoire sequences (15). In this study, we used de novo assembly to reconstruct captured Ig/TCR transcripts with full-length productive V(D)J regions instead of using germline reference sequence-based approaches. Still, the availability of comprehensive rhesus germline reference sequences could greatly facilitate the analysis and the biological interpretation of the obtained Ig/TCR repertoire data. Building rhesus germline reference sequence collections would require future developments.

In summary, these results demonstrate that scIRS assays we established in this study for rhesus macaques can be used to recover every known Ig and TCR isotype and chain type in an unbiased fashion, pair VDJ and VJ transcripts in the same cells, detect clonally expanded lineages and dual-expressing cels, and integrate V(D)J repertoires with cellular gene expression profiles. These scIRS assays will greatly facilitate the analysis of the adaptive immune responses in rhesus macaques.

Cryopreserved splenocyte samples were provided by Dr. Sallie R. Permar at Duke University. The rhesus B-lymphoblastoid cell line used in this study was provided by the National Institutes of Health Nonhuman Primate Reagent Resource (U24 AI126683). Rhesus T cell lines RH447 and RH444 were provided by Dr. Vanessa M. Hirsch at National Institute of Allergy and Infectious Diseases, National Institutes of Health.

This work was supported by the U.S. Department of Health and Human Services (HHS), National Institutes of Health (NIH), National Institute of Allergy and Infectious Diseases (NIAID) (Contract HHSN272201800008C) and by the Washington National Primate Research Center Core (Grant P51 OD010425 to M.G. and Grant R21AI120713 to X.P. from the HHS, NIH, NIAID, Office of the Director). The nonhuman primate work has also been funded by multiple awards from the HHS, NIH, NIAID, including Grants R38AI40297 (to B.I.S.) and U19AI131471 (to A.D.K.).

The 5′ gene expression sequencing data presented in this article have been submitted to the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE179722) under accession number GSE179722. Datasets generated during this study have been submitted to the NCBI Sequence Read Archive (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA746267) under project number PRJNA746267.

The online version of this article contains supplemental material.

Abbreviations used in this article

DC

dendritic cell

IRS

immune repertoire sequencing

NP

nonproliferating

P

proliferating

QC

quality-control

scIRS

single-cell immune repertoire sequencing

UMAP, uniform manifold approximation and projection; UMI

unique molecular identifier

UTR

untranslated region

1.
Schroeder
H. W.
 Jr.
,
L.
Cavacini
.
2010
.
Structure and function of immunoglobulins.
J. Allergy Clin. Immunol.
125
(
Suppl. 2
):
S41
S52
.
2.
Davis
M. M.
,
P. J.
Bjorkman
.
1988
.
T-cell antigen receptor genes and T-cell recognition. [Published erratum appears in 1988 Nature 335: 744.]
Nature
334
:
395
402
.
3.
Murphy
K.
,
P.
Travers
,
M.
Walport
.
2017
.
The generation of lymphocyte antigen receptors.
In
Janeway’s Immunobiology
, 7th Ed.
K.
Murphy
,
P.
Travers
,
M.
Walport
.
Garland Science
,
New York
, p.
188
189
.
4.
Li
F.
,
M.
Luo
,
W.
Zhou
,
J.
Li
,
X.
Jin
,
Z.
Xu
,
L.
Juan
,
Z.
Zhang
,
Y.
Li
,
R.
Liu
, et al
2021
.
Single cell RNA and immune repertoire profiling of COVID-19 patients reveal novel neutralizing antibody.
Protein Cell
12
:
751
755
.
5.
Mor
M.
,
M.
Werbner
,
J.
Alter
,
M.
Safra
,
E.
Chomsky
,
J. C.
Lee
,
S.
Hada-Neeman
,
K.
Polonsky
,
C. J.
Nowell
,
A. E.
Clark
, et al
2021
.
Multi-clonal SARS-CoV-2 neutralization by antibodies isolated from severe COVID-19 convalescent donors.
PLoS Pathog.
17
:
e1009165
.
6.
Watson
C. T.
,
J.
Glanville
,
W. A.
Marasco
.
2017
.
The individual and population genetics of antibody immunity.
Trends Immunol.
38
:
459
470
.
7.
Williams
K.
,
A.
Lackner
,
J.
Mallard
.
2016
.
Non-human primate models of SIV infection and CNS neuropathology.
Curr. Opin. Virol.
19
:
92
98
.
8.
Flynn
J. L.
,
H. P.
Gideon
,
J. T.
Mattila
,
P. L.
Lin
.
2015
.
Immunology studies in non-human primate models of tuberculosis.
Immunol. Rev.
264
:
60
73
.
9.
Meyer
C.
,
A.
Kerns
,
K.
Haberthur
,
J.
Dewane
,
J.
Walker
,
W.
Gray
,
I.
Messaoudi
.
2013
.
Attenuation of the adaptive immune response in rhesus macaques infected with simian varicella virus lacking open reading frame 61.
J. Virol.
87
:
2151
2163
.
10.
Lackner
A. A.
,
R. S.
Veazey
.
2007
.
Current concepts in AIDS pathogenesis: insights from the SIV/macaque model.
Annu. Rev. Med.
58
:
461
476
.
11.
Hansen
S. G.
,
J. C.
Ford
,
M. S.
Lewis
,
A. B.
Ventura
,
C. M.
Hughes
,
L.
Coyne-Johnson
,
N.
Whizin
,
K.
Oswald
,
R.
Shoemaker
,
T.
Swanson
, et al
2011
.
Profound early control of highly pathogenic SIV by an effector memory T-cell vaccine.
Nature
473
:
523
527
.
12.
Hansen
S. G.
,
M.
Piatak
Jr.
,
A. B.
Ventura
,
C. M.
Hughes
,
R. M.
Gilbride
,
J. C.
Ford
,
K.
Oswald
,
R.
Shoemaker
,
Y.
Li
,
M. S.
Lewis
, et al
2013
.
Immune clearance of highly pathogenic SIV infection. [Published errata appear in 2014 Nature 514: 654 and 2017 Nature 547: 123–124.]
Nature
502
:
100
104
.
13.
Munster
V. J.
,
F.
Feldmann
,
B. N.
Williamson
,
N.
van Doremalen
,
L.
Pérez-Pérez
,
J.
Schulz
,
K.
Meade-White
,
A.
Okumura
,
J.
Callison
,
B.
Brumbaugh
, et al
2020
.
Respiratory disease in rhesus macaques inoculated with SARS-CoV-2.
Nature
585
:
268
272
.
14.
Knechtle
S. J.
,
J. M.
Shaw
,
B. J.
Hering
,
K.
Kraemer
,
J. C.
Madsen
.
2019
.
Translational impact of NIH-funded nonhuman primate research in transplantation.
Sci. Transl. Med.
11
:
eaau0143
.
15.
Miller
W. P.
,
S.
Srinivasan
,
A.
Panoskaltsis-Mortari
,
K.
Singh
,
S.
Sen
,
K.
Hamby
,
T.
Deane
,
L.
Stempora
,
J.
Beus
,
A.
Turner
, et al
2010
.
GVHD after haploidentical transplantation: a novel, MHC-defined rhesus macaque model identifies CD28- CD8+ T cells as a reservoir of breakthrough T-cell proliferation during costimulation blockade and sirolimus-based immunosuppression.
Blood
116
:
5403
5418
.
16.
Brochu
H. N.
,
E.
Tseng
,
E.
Smith
,
M. J.
Thomas
,
A. M.
Jones
,
K. R.
Diveley
,
L.
Law
,
S. G.
Hansen
,
L. J.
Picker
,
M.
Gale
Jr.
,
X.
Peng
.
2020
.
Systematic profiling of full-length Ig and TCR repertoire diversity in rhesus macaque through long read transcriptome sequencing.
J. Immunol.
204
:
3434
3444
.
17.
Edgar
R. C.
2010
.
Search and clustering orders of magnitude faster than BLAST.
Bioinformatics
26
:
2460
2461
.
18.
Ye
J.
,
N.
Ma
,
T. L.
Madden
,
J. M.
Ostell
.
2013
.
IgBLAST: an immunoglobulin variable domain sequence analysis tool.
Nucleic Acids Res.
41
:
W34
W40
.
19.
Fu
L.
,
B.
Niu
,
Z.
Zhu
,
S.
Wu
,
W.
Li
.
2012
.
CD-HIT: accelerated for clustering the next-generation sequencing data.
Bioinformatics
28
:
3150
3152
.
20.
Hafemeister
C.
,
R.
Satija
.
2019
.
Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression.
Genome Biol.
20
:
296
.
21.
Goldstein
L. D.
,
Y. J.
Chen
,
J.
Wu
,
S.
Chaudhuri
,
Y. C.
Hsiao
,
K.
Schneider
,
K. H.
Hoi
,
Z.
Lin
,
S.
Guerrero
,
B. S.
Jaiswal
, et al
2019
.
Massively parallel single-cell B-cell receptor sequencing enables rapid discovery of diverse antigen-reactive antibodies.
Commun. Biol.
2
:
304
.
22.
Singh
M.
,
G.
Al-Eryani
,
S.
Carswell
,
J. M.
Ferguson
,
J.
Blackburn
,
K.
Barton
,
D.
Roden
,
F.
Luciani
,
T.
Giang Phan
,
S.
Junankar
, et al
2019
.
High-throughput targeted long-read single cell sequencing reveals the clonal and transcriptional landscape of lymphocytes.
Nat. Commun.
10
:
3120
.
23.
Tong
P.
,
A.
Granato
,
T.
Zuo
,
N.
Chaudhary
,
A.
Zuiani
,
S. S.
Han
,
R.
Donthula
,
A.
Shrestha
,
D.
Sen
,
J. M.
Magee
, et al
2017
.
IgH isotype-specific B cell receptor expression influences B cell fate. [Published erratum appears in 2017 Proc. Natl Acad. Sci. USA 114: E9750–E9751.]
Proc. Natl. Acad. Sci. USA
114
:
E8411
E8420
.
24.
De Simone
M.
,
G.
Rossetti
,
M.
Pagani
.
2018
.
Single cell T cell receptor sequencing: techniques and future challenges.
Front. Immunol.
9
:
1638
.
25.
Carter
J. A.
,
J. B.
Preall
,
K.
Grigaityte
,
S. J.
Goldfless
,
E.
Jeffery
,
A. W.
Briggs
,
F.
Vigneault
,
G. S.
Atwal
.
2019
.
Single T cell sequencing demonstrates the functional role of αβ TCR pairing in cell lineage and antigen specificity.
Front. Immunol.
10
:
1516
.
26.
Chien
Y. H.
,
C.
Meyer
,
M.
Bonneville
.
2014
.
γδ T cells: first line of defense and beyond.
Annu. Rev. Immunol.
32
:
121
155
.
27.
Espinosa
J.
,
F.
Herr
,
G.
Tharp
,
S.
Bosinger
,
M.
Song
,
A. B.
Farris
III
,
R.
George
,
J.
Cheeseman
,
L.
Stempora
,
R.
Townsend
, et al
2016
.
CD57(+) CD4 T cells underlie belatacept-resistant allograft rejection.
Am. J. Transplant.
16
:
1102
1112
.
28.
DeWolf
S.
,
B.
Grinshpun
,
T.
Savage
,
S. P.
Lau
,
A.
Obradovic
,
B.
Shonts
,
S.
Yang
,
H.
Morris
,
J.
Zuber
,
R.
Winchester
, et al
2018
.
Quantifying size and diversity of the human T cell alloresponse.
JCI Insight
3
:
e121256
.
29.
Edwards
S. C.
,
C. E.
Sutton
,
K.
Ladell
,
E. J.
Grant
,
J. E.
McLaren
,
F.
Roche
,
P.
Dash
,
N.
Apiwattanakul
,
W.
Awad
,
K. L.
Miners
, et al
2020
.
A population of proinflammatory T cells coexpresses αβ and γδ T cell receptors in mice and humans.
J. Exp. Med.
217
:
e20190834
.
30.
Reitermaier
R.
,
T.
Krausgruber
,
N.
Fortelny
,
T.
Ayub
,
P. A.
Vieyra-Garcia
,
P.
Kienzl
,
P.
Wolf
,
A.
Scharrer
,
C.
Fiala
,
M.
Kölz
, et al
2021
.
αβγδ T cells play a vital role in fetal human skin development and immunity.
J. Exp. Med.
218
:
e20201189
.
31.
Pellicci
D. G.
,
A. P.
Uldrich
,
J.
Le Nours
,
F.
Ross
,
E.
Chabrol
,
S. B. G.
Eckle
,
R.
de Boer
,
R. T.
Lim
,
K.
McPherson
,
G.
Besra
, et al
2014
.
The molecular bases of δ/αβ T cell-mediated antigen recognition.
J. Exp. Med.
211
:
2599
2615
.
32.
Japp
A. S.
,
W.
Meng
,
A. M.
Rosenfeld
,
D. J.
Perry
,
P.
Thirawatananond
,
R. L.
Bacher
,
C.
Liu
,
J. S.
Gardner
;
HPAP Consortium
, 
M. A.
Atkinson
,
K. H.
Kaestner
, et al
.
2021
.
TCR+/BCR+ dual-expressing cells and their associated public BCR clonotype are not enriched in type 1 diabetes.
Cell
184
:
827
839.e14
.
33.
Ahmed
R.
,
Z.
Omidian
,
A.
Giwa
,
B.
Cornwell
,
N.
Majety
,
D. R.
Bell
,
S.
Lee
,
H.
Zhang
,
A.
Michels
,
S.
Desiderio
, et al
2019
.
A public BCR present in a unique dual-receptor-expressing lymphocyte from type 1 diabetes patients encodes a potent T cell autoantigen.
Cell
177
:
1583
1599.e16
.
34.
Liu
Y.
,
S.
Ye
,
X.
Guo
,
W.
Li
,
Y.
Xia
,
X.
Wen
,
J.
Yu
,
Y.
Jia
,
X.
Liu
,
Y.
Guo
,
Y.
Zhao
.
2020
.
Discovery and characteristics of B cell-like T cells: a potential novel tumor immune marker?
Immunol. Lett.
220
:
44
50
.
35.
Hayes
S. M.
,
L.
Li
,
P. E.
Love
.
2005
.
TCR signal strength influences alphabeta/gammadelta lineage fate.
Immunity
22
:
583
593
.
36.
Ishida
I.
,
S.
Verbeek
,
M.
Bonneville
,
S.
Itohara
,
A.
Berns
,
S.
Tonegawa
.
1990
.
T-cell receptor gamma delta and gamma transgenic mice suggest a role of a gamma gene silencer in the generation of alpha beta T cells.
Proc. Natl. Acad. Sci. USA
87
:
3067
3071
.
37.
Bowen
S.
,
P.
Sun
,
F.
Livak
,
S.
Sharrow
,
R. J.
Hodes
.
2014
.
A novel T cell subset with trans-rearranged Vγ-Cβ TCRs shows Vβ expression is dispensable for lineage choice and MHC restriction.
J. Immunol.
192
:
169
177
.
38.
Hochstenbach
F.
,
M. B.
Brenner
.
1989
.
T-cell receptor delta-chain can substitute for alpha to form a beta delta heterodimer.
Nature
340
:
562
565
.
39.
Ichinohe
T.
,
T.
Miyama
,
T.
Kawase
,
Y.
Honjo
,
K.
Kitaura
,
H.
Sato
,
T.
Shin-I
,
R.
Suzuki
.
2018
.
Next-generation immune repertoire sequencing as a clue to elucidate the landscape of immune modulation by host-gut microbiome interactions.
Front. Immunol.
9
:
668
.
40.
Gao
F.
,
K.
Wang
.
2015
.
Ligation-anchored PCR unveils immune repertoire of TCR-beta from whole blood.
BMC Biotechnol.
15
:
39
.
41.
Heather
J. M.
,
K.
Best
,
T.
Oakes
,
E. R.
Gray
,
J. K.
Roe
,
N.
Thomas
,
N.
Friedman
,
M.
Noursadeghi
,
B.
Chain
.
2016
.
Dynamic perturbations of the T-cell receptor repertoire in chronic HIV infection and following antiretroviral therapy.
Front. Immunol.
6
:
644
.
42.
Turchaninova
M. A.
,
O. V.
Britanova
,
D. A.
Bolotin
,
M.
Shugay
,
E. V.
Putintseva
,
D. B.
Staroverov
,
G.
Sharonov
,
D.
Shcherbo
,
I. V.
Zvyagin
,
I. Z.
Mamedov
, et al
2013
.
Pairing of T-cell receptor chains via emulsion PCR.
Eur. J. Immunol.
43
:
2507
2515
.
43.
Liu
X.
,
J.
Wu
.
2018
.
History, applications, and challenges of immune repertoire research.
Cell Biol. Toxicol.
34
:
441
457
.
44.
Howie
B.
,
A. M.
Sherwood
,
A. D.
Berkebile
,
J.
Berka
,
R. O.
Emerson
,
D. W.
Williamson
,
I.
Kirsch
,
M.
Vignali
,
M. J.
Rieder
,
C. S.
Carlson
,
H. S.
Robins
.
2015
.
High-throughput pairing of T cell receptor α and β sequences.
Sci. Transl. Med.
7
:
301ra131
.

The authors have no financial conflicts of interest.

Supplementary data