Targeted reconstruction of T cell receptor sequence from single cell RNA-seq links CDR3 length to T cell differentiation state

TRAPeS

The TRAPeS algorithm has four main steps, each applied separately to the alpha and beta chains:

1. Identifying putative pairs of V and J segments. In order to recognize the V and J segments of the TCR, TRAPeS takes as input the alignment of the RNA-seq reads to the genome. TRAPeS searches for a paired-end read where one mate maps to a V segment while the other mate is mapped to a J segment, and takes those V–J pairs as putative candidates for the CDR3 reconstruction. In a case where there are no such pairings, TRAPeS takes all possible V-J combinations of V and J segments that have V–C and J–C pairing (i.e. reads where one mate maps to V or J and the other mate maps to the C segments). We note that reads are not successfully aligned to D segments of the beta chain due to their short length. Thus, for the beta chain, reconstruction of the CDR3 includes reconstruction of the D segment sequence. In addition, TRAPeS allows the user to specify the maximum number of reconstructions per chain. If the number of possible V–J pairs exceeds this number, TRAPeS ranks the pairs based on the number of reads initially mapped to them, and only attempts to reconstruct the top pairs.

2. Collecting putative CDR3-originating reads. TRAPeS finds the putative CDR3- originating reads by taking all the unmapped reads whose mates map to the V/J/C segments. In addition, since the first step of the CDR3 reconstruction includes alignment to the genomic V/J sequences (see below), TRAPeS also collects the reads that map to the V and J segments.

3. Reconstructing the CDR3. Using an iterative dynamic programming algorithm, TRAPeS extends the V and J regions. TRAPeS takes only the bases at the ends of the V and J segments closest to the CDR3 (3′ of the V segment and 5′ of the J segment). The number of initial bases is a parameter that can be tuned, set by default to min(length(V), length(J)). If the specified length is longer than the J segment, TRAPeS concatenates the J sequence to the beginning of the C sequence and use this extended segment as the initial J segment. In each iteration, TRAPeS aligns all the reads to the V and J segments separately with the Needleman–Wunsch algorithm, using the following scoring scheme: +1 for a match, –1 for a mismatch, –20 for gap opening and –4 for gap extension. In addition, we don’t penalize for having the read ‘flank’ the V and J toward their 3′ and 5′, respectively.

   Next, TRAPeS takes all the reads that aligned to the V and J segments above a certain score threshold, and build the ‘extended’ V and J sequences based on the reads. For each position, we take the base that appears in most reads as the chosen sequence for this position. This way, we extend the V and J regions in each iteration and also correct for mutations or SNPs in the known genomic V and J segments. TRAPeS repeats this step until the extended V and extended J overlap, or until TRAPeS reaches a number of predefined iterations. If no overlap is found, TRAPeS also offers an optional ‘one-sided’ mode, where it will attempt to determine the productivity (see below) of only the extended V segment. For this work, we used a threshold score of 21 for the alignment of the reads. However, in some cases a lower threshold was required, thus if no sequence was reconstructed we run TRAPeS with a scoring threshold of 15.

4. Separating similar TCRs and determining chain productivity. Since some V and J segments have similar sequence, reads can be mapped to several segments, creating few similar putative V–J pairs. In addition, two alpha or beta chains can be created within a single cell. TRAPeS takes all possible pairing and attempts to reconstruct the CDR3 region for all pairs. After reconstruction, full-length TCR sequences are created by extending the reconstructed region with the known reference sequences. Then, TRAPeS runs RSEM (14) on all reconstructed TCRs and the set of reads used as input (and their mates) in order to rank the TCRs based on the relative abundance. Next, TRAPeS determines if the TCR is productive: V and J segments are in the same reading frame and the CDR3 does not contain a stop codon. TRAPeS outputs a file with a summary of all possible reconstructions (see Supplementary Table S1 for example) for all cells, as well as separate files for each cell with the full-length TCR sequences. For this paper we used the productive chain with the highest expression as the TCR sequence for each cell.

TRAPeS is implemented in python. To increase performance, the CDR3 reconstruction using the dynamic programming algorithm is implemented in C++, and uses the seqan package (21). TRAPeS is freely available and can be downloaded in the following link: https://github.com/YosefLab/TRAPeS

TRAPeS can be easily extended to work with single-end data. The reconstruction algorithm only requires the paired-end information for the recognition of V/J segments and CDR3-originating reads, which can be easily done in single-end reads by searching for partial alignment of the read edges to the V/J segments. This feature will be available in the next TRAPeS version.

Single cell sorting

RNA sequencing

Single cell lysates were converted to cDNA following capture with Agencourt RNA Clean beads using the SmartSeq2 protocol as previously described (23). The cDNA was amplified using 22–24 PCR enrichment cycles prior to quantification and dual-index barcoding with the Illumina Nextera XT kit. The libraries were enriched with 12 cycles of PCR, then combined in equal volumes prior to final bead cleanup and sequencing. All libraries were sequenced on an Illumina HiSeq 2500 or NextSeq by either single-end 150 bp reads or short paired-end reads using the following read lengths: mouse samples—30 bp, human donor 1—26 bp for read 1 and 25 bp for read 2, human donor 2—30bp, human donor 3—26 bp. Donor 1 and donor 2 were sequenced using two batches, where every batch had cells from all of the donor's population (i.e. donor 1 batch 1 had both naive and CMV-specific cells, same for batch 2. Donor 2 batch 1 had YFV-specific, naive and effector memory cells, same for batch 2). Donor 3′s entire sample was sequenced on a single batch, and the LCMV samples from both mice were combined and sequenced on a single batch (Supplementary Table S2).

Preprocessing and Normalization of scRNA-seq data

Low quality bases were trimmed with trimmomatic (24) using the following parameters: LEADING:15, TRAILING:15, SLIDINGWINDOW:4:15, MINLEN:16. Trimmed reads were then aligned to the genome (hg38 or mm10 for human or mouse samples, respectively) with TopHat2 (25) for TCR reconstruction, and aligned to the transcriptome with RSEM (14) for transcriptome quantification.

For transcriptome analysis of the human CMV and YFV donors (donors 1 and 2), low quality cells were filtered out prior to normalization. Cells were filtered out if their read depth was less than 1 million pairs or if the cell expressed less than 20% of all expressed transcript, where a transcript was considered expressed if it had a transcripts per million (TPM) value of >10 in at least 10% of cells, leaving 353 out of 378 cells for further analysis.

Normalization of TPM values was done with our newly developed normalization framework SCONE (https://niryosef.wordpress.com/tools/scone/). SCONE considers a large number of unsupervised normalization pipelines (i.e. without using any prior biological information about samples’ origin), applying different ways to scale the data (e.g. full quantile, upper quantile) and perform factor analysis to eliminate unwanted variation. SCONE then uses a number of quality metrics to choose the best normalization, which reduces technical variation and maintain prior biological knowledge. In our study, the chosen normalization first scaled each sample with the DEseq (26) scaling factor to account for differences in sequencing depth. Then, we ran RUVg (27) with k = 1. In order to run RUVg, a list of genes that are constant across conditions should be provided. To find constant genes across the specific conditions that were tested in this paper, we also sequenced bulk populations of naive CD8⁺ T cells from donor 1 and CMV-specific effector memory CD8⁺ T cells, as well as populations of 50 cells of naive CD8⁺ T cells from donor 2 and YFV-specific effector memory CD8⁺ T cells. We ran DESeq2 (28) on those samples and defined the set of constant genes as the genes that showed no change (FDR-adjusted P-value > 0.98 and absolute log fold change < 0.2) across all pairwise comparisons (naive vs. all effector memory cells, naive versus CMV-specific effector memory, naive vs. YFV-specific effector memory and CMV-specific effector memory versus YFV-specific effector memory), resulting in a total of 373 genes.

Dimensionality reduction with PCA on samples from each donor after normalization revealed that the normalization process maintain biological information, while reducing the correlation between the data and technical variables such as batch, number of expressed genes in each cell, and the values of the first PC of the quality matrix (where the quality matrix includes for each cell technical information as previously described (29) (Supplementary Figure S10)).

Reconstructing TCR sequence from long reads

Detection of CDR3 sequence using long (150 bp) reads was performed similar to Venturi et al. (7). In short, reads were aligned against the set of known V and J segment using blastn (30). Reads with V and J segments aligning to their edges were selected, extracting the CDR3 sequence in each read. In case where more than one productive CDR3 sequence was discovered in a cell, the sequence with the highest number of supporting reads was selected.

Reconstructing TCR sequence from short paired-end reads using Trinity

Trinity (20) was run on each cell with the following parameters: –max_memory 10G, –min_contig_length 50. In addition, using the –KMER_SIZE parameter Trinity was run with four different k-mer sizes—13, 15, 17 and 19. For each k-mer size we ran Trinity twice: once in single-end mode, using the set of reads used by TRAPeS for CDR3 reconstruction, and once in paired-end mode, taking all the mapped and unmapped reads along with their pairs. Then, for each k-mer we combined the final Trinity output from both runs (paired-end and single-end) for each cell. To determine whether or not a transcript is productive and to annotate the CDR3 sequence, all possible reconstructed transcripts were run through IMGT/HighV-QUEST (31,32). We considered each productive chain output by IMGT as a successful reconstruction.

Comparing TRAPeS to TraCeR

TraCeR was run using default parameters. To compare TRAPeS to TraCeR on the benchmark data used by TraCeR (15), raw single cell RNA-seq data was downloaded as fastq files from ArrayExpress (accession number E-MTAB-3857). While the original data consisted of 100 bp paired-end reads, we converted it to that equivalent of short-read sequences by trimming each fragment to leave only the outer 25 or 30 bp of each read. We also ran TRAPeS on the original 100 bp paired-end data with the following parameters: -score 80 -bases 150 -top 15 -byExp -oneSide.

Comparing TRAPeS to scTCRseq and VDJPuzzle

VDJPuzzle was run using the default parameters. For scTCRseq, since running the software with the default parameters resulted in no alignments for human TRBV segments, we ran the software using the parameters -e 1e-7 -c 2. In addition, since scTCRseq does not summarize the data, we collected the fasta sequences of scTCRseq final results (*.gapfilled.final.vdj.fa files) and ran them through IMGT to annotate the junction sequence in each cell, taking only productive CDR3 with a complete reconstruction (no missing amino acids) as successful reconstructions. To compare TRAPeS and scTCRseq on the benchmark data used by scTCRseq (33), raw single cell RNA-seq data was downloaded as fastq files from ArrayExpress (accession number E-MTAB-2512) and trimmed from 75 bp paired-end into 25 or 30 bp paired-end. We also ran TRAPeS on the original 75 bp paired-end data with the following parameters: -score 65 -bases 150 -top 10 -bases 100

Gini coefficient calculation:

For each population, cells were considered from the same clone if they had identical CDR3 sequences of both alpha and beta chains. Cells with only one reconstructed chain were excluded from this analysis. The number of cells for each clone was counted and the Gini coefficient was calculated by using the Gini command in R from the ‘ineq’ package.

Inference of cell clusters, visualization and differential expression analysis

For cluster inference in the YFV + CMV human data, we defined an expression matrix consisting of normalized TPM values of 353 cells by 10827 transcripts (expressed at a level of ≥5 TPM in at least 1% of cells; Supplementary Table S11). We applied the SC3 software (34) for clustering the cells in this matrix using default parameters.

To visualize the data, we first used the jackStraw package (35) to reduce the dimensionality of the data and retain only principal components (PC) that are statistically significant (P-value < 10⁻⁴) in terms of the respective percent of explained variance. This analysis retained the first three PCs. We then applied t-SNE (36) with default parameters and 2000 iterations to these significant PCs, further reducing the data for visualization in two dimensions.

We used the DESeq2 package (28) to identify genes that are differentially expressed (DE) between the different clusters. In this application, each cluster was compared to the other two clusters, looking for genes that are differentially expressed. Genes were called as differentially expressed using an FDR-adjusted P-value cutoff of 0.05. The heatmap in Figure 3B was populated with log₂(TPM) values for genes identified as uniquely up- or down- regulated in each of three major phenotypic groups. We also see similar results of DE genes using the scRNA-seq analysis package Seurat (37,38). Enrichment of DE genes with respect to immunological pathways was determined using a Fisher exact test (FDR-adjusted P-value < 10⁻³) quantifying the significance of overlap between differential genes and signatures from the ImmuneSigDB database (39).

Gene enrichment by signature analysis

We used FastProject (40) together with large collection of transcriptional signatures from ImmuneSigDB (39) to characterize the phenotype of our single cells. In short, each transcriptional signature is comprised of genes that are either over-expressed or under-expressed between two cell states of interest (e.g. using published bulk RNA-seq data from naive versus memory cells). For each single cell, the signature score is computed as:

where is the signature, is the cell, for genes under-expressed in this signature and +1 for over-expressed genes, is the standardized (Z-normalized across all cells) log expression level of gene in cell , and is the estimated false-negative weight for gene in cell . To identify transcriptional signatures that are associated with an scRNA-seq data set of interest, FastProject looks for consistency between signatures and low-dimensional projections of the data. To this end, FastProject first computes a wide range of 2D projections (e.g. PCA, ICA, spectral embedding, tSNE), each capturing (possibly different) key axes of variation in the data. For each transcriptional signature and each projection it then computes a consistency score, which reflects the extent to which cells that have a similar signature score reside close to each other in the projection (thus extending our previous work (30) and facilitating the analysis of non-linear projections). The significance of the consistency score is evaluated by random shuffling.

To include only relevant signatures, we analyzed only signatures with a significant consistency score (FDR-adjusted P-value < 0.05) in at least one projection. In addition, only signatures that include ‘CD8’ in their name were used for further analysis, leaving a total of 95 signatures for the YFV + CMV human data and 154 signatures for the YFV-specific analysis.

Characterization of TCR properties of YFV-specific cells

TRAPeS reconstructs TCR sequences using short (25–30 bp) scRNA-seq

TRAPeS starts by recognizing putative pairs of V and J segments that flank the CDR3 region, using genome alignment (46) (Figure 1, top; see Materials and Methods for a complete description of the algorithm). It then identifies the set of unaligned reads that may have originated from the CDR3 region, taking the unmapped mates of reads aligned to the putative V-J segments or to the constant (C) segment (Figure 1, middle). Next, it uses an iterative dynamic programming scheme to piece together the putative CDR3 reads, gradually extending the CDR3 reconstruction on both ends (V and J) until convergence (Figure 1, bottom). Finally, after the TCR chain has been reconstructed, TRAPeS determines whether it is productive (i.e., has an in-frame CDR3 without a stop codon) and determines its exact CDR3 sequence, based on the criteria established by the international ImMunoGeneTics information system (IMGT) (47). For each cell, TRAPeS outputs a set of reconstructed TCR transcripts (from both chains), along with their complete sequence, an indication of whether or not they are productive, and the number of reads mapped to them. In some cases multiple reconstructions can be generated for the same cell. This may happen when more than one chain is produced in the cell (a phenomenon that have been previously reported (15,17,18)), or when sequence similarity between some V or J segments results in several possible V-J pairs with an identical CDR3 reconstruction. In such cases, we report all V-J pairs, while ranking the putative TCR transcripts in accordance to their estimated expression levels (Supplementary Table S1). The average running time of TRAPeS on a Human single cell library with an average two million reads per cell is less than two minutes per cell on a standard machine (Supplementary Figure S1).

TRAPeS is accurate and more sensitive than previous methods using short reads and comparable to previous methods using long reads

We applied and tested TRAPeS to scRNA-seq data from a range of CD8⁺ T cell responses (Methods, Figure 2A). These data sets were selected to include both mouse and human CD8⁺ T cells as well as those expected to have a range of TCR complexities (Supplementary Figure S2). In mice, we used the lymphocytic choriomeningitis virus (LCMV) infection model, and profiled CD8⁺ T cells responding to either acute or chronic infection (using the Armstrong and Clone 13 strains of LCMV, respectively). In healthy human subjects we profiled naive CD8⁺ T cells, effector memory CD8⁺ T cells, and antigen-specific CD8⁺ T cells elicited by CMV infection; vaccination with the live attenuated yellow fever virus infection (YFV-17D) (48); or by vaccination with adenoviral and modified vaccinia Ankara vectors encoding HCV proteins (6,22). We sorted up to 128 single CD8⁺ T cells from each dataset to a total of 565 cells, and generated scRNA-seq libraries with short (25–30 bp) paired-end reads as previously described (23,49) and observed good quality metrics using previously used measures (29) (Supplementary Table S2, Materials and Methods). To test TRAPeS, we applied cell quality filtering scheme similar to the criteria used by others (15), removing samples with <2000 genes or with >10% of reads mapping to mitochondrial genes, resulting in a total of 513 high quality cells (Figure 2A). Importantly, our results below remain consistent also when cell filtering is not applied (Supplementary Figure S3).

To evaluate the accuracy of TRAPeS, we compared its output with that of directly sequencing the TCR sequence using long reads (in which reconstruction is not required, Materials and Methods). To that end, we sequenced libraries of epitope-specific cells for Clone 13, Armstrong and CMV, and naive T cells from the CMV donor with both short (25–30 bp) paired-end and 150 bp single-end sequence reads (Figure 2A). TCR sequences identified by TRAPeS were almost perfectly consistent with those produced based on the long read data (Methods; Figure 2B and C), indicating a high level of specificity.

We compared TRAPeS to previously published methods for TCR reconstruction in single cells. First, we compared TRAPeS to TraCeR (15)—a TCR reconstruction software that is built upon Trinity (20), a de-novo transcriptome assembly tool. We found that the sensitivity of TRAPeS was markedly higher (Figure 2A–C). On average (across all data sets), TRAPeS successfully reconstructed productive alpha chains from 66% of the cells and productive beta chains from 80% of the cells, using the short (25–30 bp) libraries. In contrast, TraCeR resulted in no reconstruction for the 25 bp paired-end libraries, and was able, for the 30 bp libraries, to reconstruct CDR3 regions in an average of 43% and 15% of the cells for alpha and beta chains respectively.

Next, we considered two additional recently published methods—VDJPuzzle (18) and scTCRseq (17), both based on de-novo assembly algorithms (Trinity and GapFiller (19), respectively). As above, we observe substantially higher sensitivity with TRAPeS (Figure 2A–C, Materials and Methods). VDJPuzzle was also unable to reconstruct any productive chains in the 25 bp data and, for the 30 bp libraries, reconstructed CDR3 regions in an average of 40% and 63% of the cells for alpha and beta chains, respectively. scTCRseq, which is build upon GapFiller (19), managed to successfully reconstruct CDR3 regions in an average of 50% and 60% of the cells for alpha and beta chains, respectively. While scTCRseq achieves better results compared with Trinity-based methods, TRAPeS clearly outperforms all methods in terms of specificity and sensitivity (Figure 2A–C).

The low success rate of Trinity-based methods TraCeR and VDJPuzzle is likely due to its requirement for seed k-mer length (25nt) that is unsuitable for short reads. Thus, we also directly ran Trinity on our set of CDR3-originating reads, using a k-mer length of 13 (Materials and Methods). This resulted in an increased sensitivity for the 30 bp libraries compared to TraCeR and VDJPuzzle, but did not improve the reconstruction rates for 25 bp libraries (Figure 2A–C). Running Trinity with several other k-mer lengths (15, 17 and 19) did not significantly change the results (Supplementary Figure S4).

Notably, the average rate of successful reconstruction of TRAPeS in our mouse libraries is 93.7% (with 30 bp reads), which is higher than that achieved by TraCeR with the mouse libraries used by Stubbington et al. (86.3% with 100 bp reads) (15). To further substantiate this result, we applied TRAPeS and TraCeR on a trimmed version of this published data. We found that discarding 70–75% of the information (i.e. taking only 25 or 30 bp out of each 100 bp read) substantially hurts the performance of TraCeR, while TRAPeS is able to maintain rates of successful TCR reconstructions that are similar to those achieved in the original paper (15) (Supplementary Figure S5). Running TRAPeS on the original long read data is also comparable to the success rates obtained by TraCeR, demonstrating the ability of TRAPeS to be applied on long reads as well (Supplementary Figure S5). In addition, running TRAPeS on short or long reads is comparable to running scTCRseq using long reads, as evident by running TRAPeS on the original and a trimmed version of the data used to benchmark scTCRseq (17,33) (Supplementary Figure S6).

TRAPeS captures various clonality levels

We investigated the clonality of the TCR repertoire measured by TRAPeS among the human CD8⁺ T cells (Figure 2D, Supplementary Table S3), using the Gini Index, a clonality measure (50) ranging from zero (i.e. no two cells share the same TCR) to one (i.e. all cells are from the same clone; Methods). As expected, the naïve population had a Gini index of zero, indicating that each naive CD8⁺ T cell expressed a unique TCR. The CMV-specific CD8⁺ T cell population had a high Gini index (with 83% of CMV-specific CD8⁺ T cells with reconstructed alpha and beta chains originated from a single clone), indicating a high degree of oligoclonality as previously described (51, 52). In contrast, CD8⁺ T cells elicited by YFV or HCV vaccines showed much greater heterogeneity in TCR repertoire, consistent with a more limited, rather than persistent, exposure to antigen (6,53–56). This demonstrates the ability of TRAPeS to capture cells from the same clone even with relatively small number of antigen-specific cells, assuming a clonal response.

Single-cell transcriptome analysis detects subpopulations of YFV cells

In order to determine the relationship between TCR use and CD8⁺ T cell state, we focused on CD8⁺ T cells from two healthy donors (YFV and CMV peptide-specific, as well as naive and effector memory cells without sorting for peptide specificity; Methods) to avoid introducing additional complexity from chronic infection. To identify groups of cells with similar expression profiles, we used SC3 (34), a robust clustering method for sparse datasets, to identify subpopulations of cells (Supplementary Figure S7, Supplementary Table S4, Materials and Methods) which we then visualized using t-SNE (36) (Figure 3A). We found three clusters of cells: one that contained all CMV-specific cells (Figure 3A, purple symbols); one that contained all effector memory cells (blue symbols); and one that contained all naive CD8⁺ T cells (green symbols). In contrast to these discrete groupings, we observed that YFV-specific CD8⁺ T cells were split between two clusters: one containing effector memory CD8⁺ T cells and one containing naive CD8⁺ T cells.

Differential gene expression analysis between cell clusters revealed transcripts consistent with the known patterns of gene expression in antigen-experienced or naive CD8⁺ T cells (Supplementary Tables S5–S7, Figure 3B, Materials and Methods). CMV-specific CD8⁺ T cells expressed effector molecules and transcription factors characteristic of antigen experienced cells (e.g., Granzyme B, PRDM1), which were not detected in naive cells. Naive CD8⁺ T cells expressed canonical markers of the naive state (CCR7, SATB1, LEF1) that were absent in CMV-specific and effector memory CD8⁺ T cells. The expression of these genes in YFV-specific CD8⁺ T cells was consistent with the cluster in which they were associated, with those in the naive cluster expressing minimal Granzyme B or PRDM1, but showing robust expression of CCR7, SATB1, and LEF1 (Supplementary Figure S8).

To identify broader patterns of transcriptional signatures, we applied FastProject (40)—a software tool that enables the expression of gene sets of interest to be quantified in transcriptional profiles of single cells (Materials and Methods). We surveyed the enrichment of a collection of gene sets, from the C7 (ImmuneSigDB) (39) collection of MSigDB (57) corresponding to cell states and perturbations of CD8⁺ T cells. We found significant up-regulation of multiple gene sets corresponding to naive CD8⁺ T cells (K–S test FDR- adjusted P-value<0.01) in the naive cluster (cluster 3) compared to the other two clusters. Consistent with this, we found significantly greater up-regulation of effector signatures in clusters 1 and 2 compared with the other clusters (FDR-adjusted P- value < 0.01; Figure 3C and Supplementary Table S8).

To confirm these patterns of transcript abundance at the protein level, we compared flow cytometry data for a set of surface markers acquired at the time of sorting (Methods) with transcript abundance in the same cell (Figure 3D). Consistent with the gene expression profiles, we observed that YFV-specific CD8⁺ T cells in the naive-like cluster (open symbols) showed higher protein levels of CCR7 and CD45RA than those in the effector memory cluster (purple symbols). Thus, single-cell analysis shows that CD8⁺ T cells specific for the same peptide epitope from YFV are heterogeneous and includes both effector-memory and naive-like gene expression profiles, as has been reported previously for cells analyzed at the bulk level (4,5,58).

Combined TCR-transcriptome analysis reveals longer CDR3 regions for naive-like YFV-specific cells

We reasoned that differences in TCR might contribute to the heterogeneous differentiation of CD8⁺ T cells following YFV vaccination. To that end, we evaluated a number of properties to characterize each reconstructed TCR—CDR3 specific properties such as length, hydrophobicity and germline score as well as TCR expression. In addition, we measured the normalized tetramer staining intensity per cell (Supplementary Table S9, Methods). We then asked whether any of those properties differed between naive-like and effector memory-like YFV-specific CD8⁺ T cells. Naive-like and effector memory-like YFV-specific CD8⁺ T cells were indistinguishable (P-value>0.05, FDR-adjusted P-value > 0.1) in terms of TCR transcript expression, hydrophobicity of the CDR3 region and normalized tetramer staining intensity (Materials and Methods). However, we found that the CDR3 sequence was significantly longer in YFV-specific CD8⁺ T cells with a naive-like state compared with those with an effector memory profile for both alpha and beta chains (Figure 4A, K–S test P-value 0.038 and 0.027 for alpha and beta chains respectively, FDR-adjusted P-value = 0.084 for both alpha and beta chains).

We next evaluated the germline score of CDR3 regions in YFV-specific CD8⁺ T cells, a measure of the contribution of germline nucleotides to the CDR3 region. The germline score is defined as the ratio between the number of nucleotides in the CDR3 that originate from the germline (V, D, J segments) to the total number of nucleotides in the CDR3 (16) (Materials and Methods). Consistent with the differences in the CDR3 length, we found that naive-like YFV-specific CD8⁺ T cells had a significantly lower germline score in both alpha and beta chains than did effector memory-like cells (Figure 4B, K–S test P-value of 0.034 and 0.029 for alpha and beta chains respectively, FDR-adjusted P-value 0.084 for both alpha and beta chains), suggesting that generating the CDR3 region of these TCRs involved a greater degree of nucleotide addition/subtraction.

To further characterize the relationship between CDR3 length and cellular state in YFV-specific CD8⁺ T cells, we identified CD8⁺ transcriptional signatures (extracted from ImmuneSigDB (39) and scored with FastProject (40), as above) that correlated with CDR3 length across all YFV-specific CD8⁺ T cells (Supplementary Table S10, Materials and Methods). Of all signatures evaluated, we found that only naive CD8⁺ T cell signatures showed a significant positive correlation with CDR3 length (FDR-adjusted P-value < 0.1; Figures 4C and D). Previous work has suggested that YFV-specific CD8⁺ T cells with a naive-like phenotype include those with a stem-cell memory (Tstem-memory) differentiation state. We found that Tstem-memory signatures were more enriched in naive-like YFV-specific CD8⁺ T cells than in effector memory YFV-specific CD8⁺ T cells (Supplementary Figure S9). However, the enrichment for these signatures was equivalent between naive-like YFV-specific and phenotypically naive CD8⁺ T cells, making it difficult to discern whether these cells manifest a specific stem-cell-like state. Our results, however, show that heterogeneity in the differentiation state of CD8⁺ T cells responding to a single epitope of YFV is strongly associated with the CDR3 length.