Allelic variation in class I HLA determines CD8+ T cell repertoire shape and cross-reactive memory responses to SARS-CoV-2.

Effective presentation of antigens by human leukocyte antigen (HLA) class I molecules to CD8+ T cells is required for viral elimination and generation of long-term immunological memory. In this study, we applied a single-cell, multiomic technology to generate a unified ex vivo characterization of the CD8+ T cell response to severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) across four major HLA class I alleles. We found that HLA genotype conditions key features of epitope specificity, TCRα/β sequence diversity, and the utilization of pre-existing SARS-CoV-2-reactive memory T cell pools. Single-cell transcriptomics revealed functionally diverse T cell phenotypes of SARS-CoV-2-reactive T cells, associated with both disease stage and epitope specificity. Our results show that HLA variations notably influence the CD8+ T cell repertoire shape and utilization of immune recall upon SARS-CoV-2 infection.

INTRODUCTION

Elicitation of a robust and durable neutralizing antibody response following immunization of large sections of the population with approved SARS-CoV-2 vaccines is limiting viral transmission and decreasing mortality, providing hope that the global threat from the COVID-19 pandemic is diminishing. However, the appearance of new viral variants warrants continued vigilance. A more complete understanding of the underlying cellular mechanisms that regulate host immunity and contribute to long term protection is required. Infection with SARS-CoV-2 leads to an upper respiratory tract infection, which can be benign or even asymptomatic. If not controlled by the immune response, it can evolve into a lethal pneumonia with immunopathology due to excessive amplification of the innate inflammatory response, complicated by several extra-respiratory manifestations (). While humoral responses play an important role in immunological control of infection, the generation of effective cellular immunity and expansion of cytotoxic CD8+ memory T cells is also required to eliminate virally infected cells as shown from the earlier SARS-CoV-1 epidemic, even in the absence of seroconversion (–).

Several recent studies have focused on the discovery of relevant SARS-CoV-2 epitopes in both CD4+ and CD8+ T cell responses, leveraging in silico predictions, stimulation/expansion with peptide pools (–), tetramer binding (, ), and analysis of presentation in vitro (). Collectively, these studies identified a number of immunodominant epitopes derived from across the viral proteome including structural and non-structural proteins in canonical (–) and non-canonical open reading frames (). Interestingly, some of these specificities were also detected in uninfected individuals, suggesting potential cross-reactivity from endemic human coronaviruses (HCoV) to which the population is routinely exposed (), though a direct connection to pre-existing memory cells has not been established.

The breadth and nature of the cellular immune response to SARS-CoV-2 infection is driven by diversity in both T cell receptor (TCR) repertoire and human leukocyte antigen (HLA) genetics. Mammalian cells express up to six different HLA class I alleles that shape antigen presentation in disease, and allelic diversity has been associated with both disease susceptibility and outcome of viral infections (, ). There are divergent reports regarding HLA polymorphism and COVID-19 incidence and severity, although the major genome-wide association studies clearly show no dominant effect of the locus (–). Together with genetic influences on HLA-associated antigen presentation, the clonal selection of TCRs that compose an individual’s repertoire contributes to the nature and dynamics of the antiviral response, including cellular cytotoxicity and memory formation. Interestingly, despite a potential TCR diversity of 1015 (), several studies have described “public” T cell responses in COVID-19, where complementarity-determining region (CDR) sequences are conserved within and across individuals (, ). The extent to which TCR diversity, especially in the context of epitope specificity restricted to HLA, contributes to response is not well understood.

Here, we leverage an assay technology to elucidate, at single-cell resolution, the connection between T cell specificity, HLA variation, conserved features of paired α/β TCR repertoires, and cellular phenotype observed in CD8+ T cell responses to SARS-CoV-2 infection. We profiled 96,909,416 CD8+ T cells ex vivo across 78 samples from acute, convalescent, or unexposed individuals, and identified T cell specificity to 648 epitopes presented by four HLA alleles across the SARS-CoV-2 proteome, few of which are implicated by the current variants of concern. Estimated frequencies of epitope-specific CD8+ T cells observed in convalescent patients had a mean value of 0.01% and maximum around 1% of the total CD8+ T cell population. We observed that TCR repertoires were surprisingly public in nature, though we found a high degree of pre-existing immunity associated with a clonally diverse response to HLA-B*07:02, which can efficiently present homologous epitopes from SARS-CoV-2 and HCoVs. Transcriptomic analysis and functional validation confirmed a central memory phenotype and TCR cross-reactivity in unexposed individuals with HLA-B*07:02. Our data suggest an association between HLA genotype and the CD8+ T cell response to SARS-CoV-2, which may have important implications for understanding herd immunity and elements of vaccine design that are likely to confer long-term immunity to protect against SARS-CoV-2 variants and related viral pathogens.

RESULTS

Direct ex-vivo detection and decoding of SARS-CoV-2-specific CD8+ T cells

We leveraged single-cell RNA-sequencing with DNA-encoded peptide-HLA tetramers to characterize CD8+ T cell responses to SARS-CoV-2 across multiple Class I alleles in subjects with varying degrees of disease severity. The technology simultaneously determines the specificity of paired α/β TCR sequences for HLA-restricted epitopes, provides transcriptomic phenotype at single-cell resolution, and gives an indication of T cell frequency (). To interrogate TCR specificity in a highly multiplexed fashion, we generated libraries of HLA-conjugated streptavidin tetramers (Fig. S1) using ultraviolet-mediated peptide exchange with subsequent labeling via biotinylated DNA barcodes. The encoded, peptide-exchanged tetramers showed similar staining to conventional tetramers (Fig. S2) while enabling readout via single-cell sequencing.

Fig. 1.

Overview of the experimental approach to decode the CD8+ T cell response to SARS-CoV-2.

(A) Schematic of the method where encoded tetramer libraries, designed independently for each HLA allele to span the entire SARS-CoV-2-proteome, are used to stain enriched CD8+ cells from subject PBMCs, which are then sorted and subjected to single-cell sequencing (left). Using this approach, TCR sequence, peptide/HLA specificity and transcriptomic features are simultaneously acquired for each cell (right). (B) Clonotype specificity detected by HLA allele and epitope across the SARS-CoV-2 proteome in COVID-19 patient samples (n=61). A scheme of the viral ORF structure is shown at the top. Bar colors denote HLA allele. Amino acid sequences of epitopes recognized by the largest number of T cell clonotypes are shown next to the corresponding bar. (C) Single-cell transcriptomic analysis showing global uniform manifold approximation and projection (UMAP) clustering, scoring by functional gene set, and projections onto the transcriptomic UMAP for T cells with specificity toward select epitopes in convalescent individuals (n=33 samples). QYI-A24, PTD-A01, and LLY-A02 correspond to QYIKWPWYI in A*24:02, PTDNYITTY in A*01:01, and LLYDANYFL in A*02:01, respectively.

We designed peptide-HLA tetramer libraries to ensure comprehensive coverage of SARS-CoV-2 and related human coronaviruses across four class I HLA alleles prevalent in North America (A*02:01, B*07:02, A*01:01, and A*24:02, hereafter A*02, B*07, A*01, and A*24). Library inclusion was determined computationally using predicted HLA binding (NetMHC-4.0 ()) of candidate peptides from a set of all possible 9-mers from the SARS-CoV-2 proteome (40% from structural, 60% from non-structural proteins), potentially immunogenic neopeptides from known SARS-CoV-2 variants, and immunogenic epitopes from SARS-CoV-1. A total of 1,355 SARS-CoV-2-related epitopes were included in the libraries in addition to well-characterized epitopes from common endemic viruses (cytomegalovirus (CMV), Epstein-Barr virus (EBV), and influenza).

The peptide-HLA tetramer libraries were used to interrogate peripheral blood mononuclear cells (PBMCs) from individuals who had been infected with SARS-CoV-2 (n = 28 convalescent, n = 27 with acute disease that required hospitalization), or who were unexposed (n = 23) (summarized in Data file S2). For each sample, CD8+ cells were isolated from PBMCs using magnetic separation, incubated with HLA-matched tetramer libraries, and sorted by flow cytometry to enrich viable, tetramer positive cells (Fig. S3). Sorted single cells were encapsulated with DNA-encoded hydrogel beads to provide cell-specific barcodes and unique molecular identifiers (UMIs) that could be used to unify reads across independent sequencing libraries for TCR, peptide-HLA tetramer, and messenger RNA (mRNA) (). We determined the specificity of TCRs using a classification method that identified UMI counts for TCR-peptide-HLA interactions that were outliers when Z-score transformed within and across cells for each sample. The resulting classifier was evaluated against functional assay data for each allele by a receiver-operator curve (ROC) analysis to identify thresholds, which were then used for normalization. The normalized classifier evaluated by ROC analysis provided an area under the curve (AUC) of 0.82 (Fig. S4), and at a threshold of 1, yielded a true positive rate of 93% and a false positive rate of 32%.

From the 52,728,647 CD8+ cells interrogated from acute and convalescent COVID-19 patients, we identified high-confidence TCR-peptide-HLA interactions across 434 immunogenic SARS-CoV-2-derived epitopes and 1,163 independent α/β TCR clonotypes (, Data file S3). The immunodominant epitopes we identified ex vivo were consistent with those measured by other means (–), but we also identified many epitopes with less dominant representation (yet observed with two or more reactive clonotypes), 178 of which had not been previously reported as minimal epitopes (Data file S4). Importantly, CD8+ T cell reactivity to SARS-CoV-2 epitopes was observed across the entire proteome, generally distributed in a manner consistent with protein lengths (Table S1). Of relevance, 85 of these epitopes were derived from the Spike protein currently used in vaccines, but only five of them (a total of 35 CD8+ T cell clonotypes in our study) would be affected by the delta (B.1.617.2) SARS-CoV-2 variant (Data file S5). Interestingly, the epitope NYNYLYRLF in A*24 (NYN-A24) is one of those epitopes, which was detected in 71% of convalescent patients with A*24, but at a modest mean frequency of 0.025%.

Dimensionality reduced projections of mRNA expression for 224,780 CD8+ T cells revealed the broad phenotypic variance observed within this study spread across 8 clusters (). We defined the phenotypic features of clusters using gene signatures generally associated with various CD8+ T cell states, including those with naïve, memory, effector, and proliferative status (). In this space, cells from convalescent patients that recognized different dominant epitopes were commonly associated with divergent phenotypes, as shown for representative epitopes (). For example, T cells specific for QYIKWPWYI in A*24 (QYI-A24, Spike) were clustered in regions with high effector scores while those specific for PTDNYITTY in A*01 (PTD-A01, Orf1ab) and LLYDANYFL in A*02 (LLY-A02, Orf3a) resided at opposite ends of memory-rich regions. Thus, and as will be further detailed below, the different immunoreactive epitopes of SARS-CoV-2 elicit distinct CD8+ T cell phenotypes.

Evolution of immunoreactivity through COVID-19 disease progression

Having established a broad landscape of SARS-CoV-2-reactive CD8+ T cells, we asked how TCR repertoires evolve over the course of infection and recovery. As our approach does not require cell expansion to determine TCR specificity, we were able to estimate the frequency of epitope-specific cells in the CD8+ T cell pools of convalescent, acute, and unexposed individuals (see Methods). While these calculated frequencies are not precise measurements, they are indicative of T cell frequency in parental samples, reflecting robustness of response observed by subject, HLA, and epitope. We show the frequency, for each subject, of T cells reactive to the top five epitopes detected across each of the four HLA variants analyzed (). Notably, we detected markedly fewer SARS-CoV-2-specific T cells in patients with acute disease compared to those in convalescence, which may have been impacted by the lower recovery of T cells that could subsequently be profiled (Fig. S5). The striking reduction also applied to memory T cells from prior antiviral responses in these patients, including influenza and EBV, but potentially less to the CMV-specific pool in multiple acute subjects (Fig. S6). The paucity of virus-reactive T cells is consistent with the T cell lymphopenia that has been reported to occur in patients with acute COVID-19 (, ).

Fig. 2.

CD8+ T cell specificity for major SARS-CoV-2 epitopes across HLA, cohort, and subject.

(A) Estimated frequency of T cell response in parent CD8+ T cell samples by subject and cohort. Heat maps show mean frequency of epitope specificity calculated for the top five epitopes (x-axis, ranked by cumulative frequency in convalescent patients) by subject (y-axis) across each cohort. Sample counts for convalescent, acute, and unexposed cohorts: A*02:01 (n=18,18,7), A*24:02 (n=7,6,5), A*01:01 (n=3,1,5), B*07:02 (n=5,3,9). (B) Prevalence of T cell specificity observed in unexposed versus convalescent cohorts represented as percentage of cohort with any detectable frequency of T cell specificity against each epitope. Dot sizes represent the mean combined estimated frequency of an epitope specificity across convalescent and unexposed subjects. Sequences of select epitopes from each HLA are annotated. (C) Sequence alignment between select SARS-CoV-2 epitopes and related common cold coronaviruses (HCoV) epitopes. Mismatches are represented in red and HLA anchor residues with a grey background. Green arrows indicate sequences where anchor and all internal residues are conserved between SARS-CoV-2 and HCoV species.

We also observed that the frequencies of SARS-CoV-2-specific T cells in unexposed individuals varied markedly with the HLA allele (). While the top five dominant epitopes in HLA-A*02, A*24, and A*01 were generally associated with measurable responses across convalescent subjects (), the frequency of responses was significantly lower in unexposed compared to convalescent individuals (p=1.5e-5, 8.5e-4, 1.7e-6 by Wilcoxon rank-sum, respectively). In stark contrast, there was no discernible difference in response frequency detected across the most immunodominant epitopes in B*07:02 individuals (p=0.23). In fact, CD8+ T cells recognizing nucleocapsid-derived SPRWYFYYL in B*07 (SPR-B07) were found in almost 80% of unexposed subjects (n=9 samples) with a mean frequency of 0.01% (), presaging the immunodominance of this epitope in convalescent COVID-19 patients, where reactivity was detected in 100% of the samples (n=5) at a mean frequency of 0.29%.

The broad presence of SARS-CoV-2-specific T cells in unexposed B*07 subjects could originate from fortuitous cross-reactivity of a public specificity, or from priming via previous exposure to a highly related endemic HCoV. Indeed, SPR-B07 shows marked homology to the corresponding segments of the nucleocapsid proteins from multiple prevalent HCoVs, including HKU1 and OC43, with only a single amino acid residue mismatched at the N terminus (). The nature of the homology preserves internal TCR-contact residues as well as the P and L anchors for HLA binding in peptide positions 2 and 9. Accordingly, the HCoV epitope (LPRWYFYYL, LPR-B07) is predicted to bind with high affinity to HLA-B*07 and could reasonably be expected to cross-react with SPR-B07-specific TCRs. Broader sequence alignment with HCoVs revealed very little homology to the immunodominant epitopes of A*02 and A*01 but did identify a complete match to VYIGDPAQL for A*24 (VYI-A24, Orf1ab). Surprisingly, T cell specificity to VYI-A24 was not detected in a single unexposed subject (n=5 samples). This likely reflects the lower frequency of response elicited by this epitope or an insufficient commitment to memory following exposure to HCoVs. Overall, we found that the response to SARS-CoV-2 is sharply distinguished by HLA genotype, as can be seen clearly in the case of A*02 and B*07, where it appears that highly specific CD8+ responses are either generated de novo or amplified from an abundant pre-existing pool, respectively.

Functional reactivity and cross-reactivity of SARS-CoV2-specific clonotypes

To confirm the specificity and functionality of TCR-peptide-HLA interactions identified in this study, we cloned several of the observed α/β TCRs clonotypes and expressed them in the TCR-null Jurkat J76 cell line (). Activation of these transductants upon stimulation by SARS-CoV-2 peptides, presented by an HLA-matched lymphoblast cell line, was evaluated by measuring the induction of surface CD69 (). Altogether, we validated 28 interactions for epitopes derived from Orf1ab, Spike, Nucleocapsid, Membrane, and Orf3a proteins, spanning high confidence interactions observed across multiple cells as well as interactions observed exclusively in single cells (Data file S6). Dose-response curves for a subset of interactions in A*02 and B*07 are shown in . The effective concentrations (EC50s) measured for these interactions ranged from 1 to 100 nM, with no particular relationship to epitope immunodominance or clonotype frequency measured ex vivo from the respective subject. These values are consistent with interactions measured for CMV-specific epitopes in A*02 using the same system. We next used these recombinant TCR expressing cell lines to compare the functional reactivity elicited by homologous epitopes from HCoVs (). Activation was insignificant for the closest homologs of Orf3a-derived LLY-A02 and Orf1ab-derived ALWEIQQVV in A*02 (ALW-A02, Orf1ab), all of which actually originated from HCoV spike proteins. In contrast, HKU1 and OC43 homologs of nucleocapsid-derived SPR-B07 and KPRQKRTAT in B*07 (KPR-B07, Nucleocapsid) epitopes drove substantial T cell activation ().

Fig. 3.

Functional confirmation of identified epitope/HLA interactions with clonotypic TCRs and comparison of SARS-CoV-2 epitopes with common cold coronaviruses.

(A) Schematic showing lentiviral transduction of clonotypic recombinant TCRs (rTCRs) into J76 cells, loading of antigen presenting cells (APC) with synthetic peptide, and quantification of activated J76 cells expressing surface CD69. (B) Dose-response curves for TCR-pHLA interactions observed across several canonical epitopes in A*02:01 and B*07:02 backgrounds. Shown are fractions of CD69+ cells after a 16-hour stimulation. Error bars represent standard deviation across replicate measurements. (C) Functional activation of TCRs by canonical and homologous epitopes, represented as fraction of CD69+ cells after 16-hour stimulation with 10 μM peptide. (D) Dose-response curves for several rTCRs from COVID-19 patients (left) or unexposed subjects (right) stimulated with peptides from SARS-CoV-2 or HCoV HKU1/OC43. Error bars represent standard deviation across replicate measurements.

We further assessed the sensitivity of B*07 interactions, comparing the reactivity of SPR-B07-specific clonotypes identified from COVID-19 patients or unexposed subjects to SARS-CoV-2-derived SPR-B07 or HCoV-derived LPR-B07 (). The three TCRs identified from COVID-19 individuals yielded EC50s that were essentially identical for the two epitopes, all falling between 50-100 nM (, left). Two of the TCRs from unexposed individuals yielded EC50s in the same range, again comparable for the HCoV and SARS-CoV-2 variants, while a third showed a >10-fold preference for the HCoV epitope (even though it was originally detected as binding to the SARS-CoV-2 peptide). Aside from providing validation that the specificities detected in our barcoded tetramer technology indeed correspond to antigen-reactive T cells, these findings support that the homologies between SARS-CoV-2 and HCoV epitopes are functionally relevant, and that pre-existing cellular reactivity to SARS-CoV-2 in B*07 subjects likely result from previous exposure to HCoVs like HKU1 or OC43.

HLA Restricted SARS-CoV-2 Epitopes Impact V(D)J Gene Usage

Given the comprehensive landscape of TCR specificity determined with our approach, we sought to elucidate the extent to which TCR usage is shared within and across subjects. We examined the linkage between paired TCR α/β sequences and their epitope specificity to determine if any features are implicated in the CD8+ T cell response to SARS-CoV-2. We used TCRs from 2,469 SARS-CoV-2-specific T cells to perform network mapping of epitope-specific subsets across several immunodominant epitopes identified (). Importantly, because it is known that during development, a TCR β-chain can be paired with many different α-chains, the network analysis allowed clonotype linkages by α or β complementarity-determining region 3 (CDR3) sequences (indicated by edges), identifying conserved motifs based on physicochemical similarity (via BLOSUM matrices) within in the epitope specific T cell population (). T cells from COVID-19 patients that recognize the most dominant A*02-, A*24-, and A*01-restricted epitopes, which have no counterpart in unexposed repertoires, showed a high degree of motif sharing with the exception of KLWAQCVQL in A*02 (KLW-A02, Orf1ab) (). Interestingly, all of these epitopes, including KLW-A02, show dominant usage of a single TCR alpha variable (TRAV) region, and in the cases of QYI-A24 and PTD-A01, dominant usage of both TRAV and TCR beta variable (TRBV) regions (). In marked contrast, SPR-B07-specific T cells, including those that also recognize homologs from HCoV, were far more diverse in CDR3 across subjects (), using 8 TRAV and 3 TRBV regions to cover 50% of the clonotypes represented. We observed two instances of CDR3 homology shared across cohorts, as indicated by the presence of nodes with unconnected edges, which are represented in both network maps.

Fig. 4.

Analysis of TCR sequences from cells specific for the most immunodominant epitopes for each HLA allele tested.

(A) Network plots showing biochemical similarity of TCRα or TCRβ CDR3 regions in unexposed subjects (left) or COVID-19 patients (right). Unique subjects are identified by node color. Each node is a unique clonotype within a subject, and the size of the node represents the relative frequency of the response detected. Edges drawn between nodes represent CDR3 homology, and the size of each node represents relative cell frequency. (B) TCRα or TCRβ (V) gene usage across all sequences represented in (A) with the most frequently used gene labeled. (C) Distributions of CDR3 lengths. (D) Paired CDR3 motifs for the most interconnected nodes identified in the network analysis (A).

These comparisons show that the reactivities that appear during SARS-CoV-2 infection may stem from both the amplification of highly related TCRs, or from the usage of diverse pre-existing T cell populations. This conclusion extended to CDR3 lengths (), which were tightly distributed for α− and/or β- chains in T cells reactive to the top epitopes in A*02, A*24, and A*01, but substantially less so for SPR-B07. To further elucidate the extent of the public nature of paired α/β TCR usage in COVID-19, we generated consensus sequences from select interconnected network clusters (). This representation provides insight into α/β linkage in the context of public responses that cannot be afforded by bulk sequencing approaches. Most motifs were represented by multiple sequences and shared by at least 40% of the subjects studied, with the exception of KLW-A02 that was shared across only 22%, and SPR-B07 that was shared across only 14%, notably with identical α/β sequences (). Thus, we have observed divergent TCR repertoire utilization, conditioned by HLA and the presence of diverse, pre-existing reactivity resulting from prior viral exposure.

CD8+ Memory T cell Phenotypes vary with recognition of SARS-CoV-2 epitopes

To examine how CD8+ T cell phenotype varied in relation to disease status, HLA/epitope specificity, and TCR diversity, we performed a more detailed analysis of the single-cell transcriptomic data. We leveraged, as an internal reference, the transcriptomic phenotype of T cells reactive to common acute and latent infections, including influenza, EBV, and CMV. To relate these data to existing knowledge on differentially expressed genes that delineate CD8+ T subsets, we used supervised partition clustering based on imputed expression (see Methods) of a set of 51 curated transcripts (Table S2) characteristic of naïve, memory, effector, or chronically-activated/exhausted populations (). This resulted in the identification of seven distinct cell clusters. Some were easily assigned (naïve cells in C1, central memory in C2, and fully activated cytotoxic effectors in C7). Other memory/effector intermediates were more tentatively labeled, as they did not easily fit into existing categorizations (–). These included a population (C3, here “CD127+ Memory”), which expresses markers of naïve, memory and effector cells, and 3 other clusters with characteristics of memory or chronically activated cells (C4-6).

Fig. 5.

Transcriptomic clustering of T cells based on function-specific gene sets.

(A) Single cell gene expression heat map of single CD8+ T cells specific for SARS-CoV-2, CMV, EBV, Influenza, or with no observed (N.O.) specificity. Units are log transformed, normalized counts per 10,000 UMIs per cell plus one, ln(CP10K+1). Kmeans clustering identified seven distinct clusters showing gene expression consistent with a range of functional states. (B) Ticks indicate the location and cohort assignment of individual cells with specificity for epitopes from SARS-CoV-2, CMV, EBV, influenza or no observed (N.O.) specificity. (C) Gene expression of single cells with individual epitope specificities indicated in the heat map. In cases where specificity was detected in the unexposed cohort, pie charts are shown to indicate the fraction of cells corresponding to each cluster identified in (A).

SARS-CoV-2 specific T cells were found in all clusters (, bottom), but at proportions that varied with stage of disease and epitope specificity (Data file S7). Cells from acute patients predominantly showed effector phenotypes, but also paradoxically naïve types. In convalescent donors, T cells from several epitope specificities were broadly distributed, consistent with the resolution of an infection. Several epitope-specific T cell pools were predominantly found in central memory (C2), including PTD-A01 (49%) and LLY-A02 (42%), while others predominantly resided in the cytotoxic terminal effector cluster (C7), including TLMNVITLV in A*02 (TLM-A02, Orf1ab, 80%), and LLLDRLNQL in A*02 (LLL-A02, Nucleocapsid, 61%) (). In most other reactivities, including SPR-B07, transcriptional profiles in convalescent patients were fairly broadly distributed across all clusters. In contrast, the reactivity in unexposed subjects was dominated by the central memory pool, confirming that the CD8+ cells likely result from long-term exposure to cross-reactive antigens. This was especially clear in the case of B*07, where epitope-specific T cells for SPR-B07, QAWQPGVAM in B*07 (QPG-B07, Orf1ab), and SIIAYTMSL in B*07 (SII-B07, Spike) were represented in central memory (C2) at proportions of 88%, 75%, and 67%, respectively. Other notable reactivities associated with central memory include TSQWLTNIF in A*24 (TSQ-A24, Orf1ab, 70%) and NSSTCMMCY in A*01 (NSS-A01, Orf1ab, 68%), though the source of these memory cells, like QPG-B07 and SII-B07, does not appear to be from HCoV exposure based on a lack of homology. Overall, this analysis provides further evidence that SPR-B07 responses to SARS-CoV-2 are likely drawn from a cross-reactive memory pool and that commitment to different cell fate is dependent on epitope specificity.

We also observed an interesting dynamic between SARS-CoV-2 infection and existing T cell pools specific for common viral infections, with differentiated outcomes likely shaped by exposure history (). Influenza-specific CD8+ T cells, which result from vaccination or past infections, mapped primarily to the central memory (C1) and effector memory (C3) compartments in unexposed individuals. Proportions were stable across epitope specificities in COVID-19 patients with the exception of GILGFVFTL in A*02 (GIL-A02, Influenza A Matrix protein 1), where the proportion of effector memory cells decreased from 50% to 0% and a naïve population representing 30% of the cells paradoxically emerged. CMV- and EBV-specific T cells, likely subject to more chronic stimulation from low-level re-activation of these integrated herpesviruses, mapped to more activated pools in unexposed subjects, as has been described by others (). After SARS-CoV-2 infection, EBV-specific cells shifted markedly from central memory (C2) and chronically stimulated compartments (C5) into the CD127+ memory cluster (C3). These changes may reflect either bystander activation, perhaps as a result of the high cytokine release in COVID-19 patients, or from changes in homing or recirculation patterns that bring into the blood cells normally sequestered in tissues. These observations suggest that, in addition to inducing lymphopenia, COVID-19 strongly reshuffles third-party antiviral T cell pools, the extent of which may be associated with exposure history and, at least to some degree, epitope specificity.

DISCUSSION

Here we presented a unified description of the CD8+ T cell response to SARS-CoV-2, highlighting the importance of HLA genetics, TCR repertoire diversity, and epitope-specific navigation through a complex transcriptomic phenotype at various stages of disease. In building a comprehensive map of immunodominant, HLA-restricted epitopes broadly derived from proteins across the entire SARS-CoV-2 proteome, we highlight how only some HLA haplotypes are associated with the existence of a pre-existing CD8+ T cell memory pool in unexposed individuals. We further show how HLA variation plays an important role in shaping the diversity of CD8+ T cell repertoires upon exposure to SARS-CoV-2, and that cellular phenotype and commitment to memory can be associated with epitope-specificity in the context of both SARS-CoV-2 and latent EBV infections.

The presence of SARS-CoV-2 reactive CD8+ T cells has been linked to milder disease (, , ), although the precise link between cellular immunity and host protection still remains to be further understood (, , ). We found that individuals carrying HLA-B*07 show a CD8+ T cell response that is dominated by pre-existing memory pools reactive to multiple SARS-CoV-2 epitopes, especially SPR-B07, which is likely induced by previous exposures to benign HCoVs. In contrast, the immunodominant responses in A*02 individuals (e.g., to YLQPRTFLL in A*02 (YLQ-A02, Spike) and LLY-A02) appear to be driven largely by the expansion of antigen-inexperienced SARS-CoV-2-specific T cells. It is interesting to note that CD8+ T cell cross-reactivity may be less widespread in unexposed individuals than for CD4+ T cross-reactivity, for which ~50% of unexposed individuals exhibited CD4+ T cell memory (). Our data provides a basis for this limited representation of the CD8+ T cell repertoire in that only a subpopulation of individuals carrying a specific HLA allele would have these cross-reactive memory CD8+ T cells. The extent to which pre-existing memory specific to SPR-B07 contributes to protection would need to be explored with longitudinal studies spanning SARS-CoV-2 exposure.

The interplay between HLA-restricted epitope presentation and available TCR repertoire shapes the cellular response to SARS-CoV-2. There are few limited studies suggesting an influence of HLA genotype on COVID-19 severity (, –). Large-scale studies evaluating T cell responses across a comprehensive HLA coverage per patient may help identify or deconvolute relationships between HLA genotype, like B*07 in this study, and protection against severe disease, ideally uncovering mechanism. Here, we observed an interesting connection between TCR repertoire diversity and HLA restriction. Responses seen in A*02, A*24, and A*01 were more often associated with “public” CDR3 motifs and consistent V gene segment usage in the α− and/or β− chains. In contrast, the dominant immune response in B*07 leveraged a significantly more diverse TCR repertoire. Several contributors to public TCR responses have been proposed, focusing on the physicochemical features of HLA-restricted peptides (e.g., “featureless” peptide-HLAs may drive a public response) and convergent recombination of TCR sequences (). The method described in this work provides an ideal system to address this question. Perhaps counterintuitively, our results show that in the case of COVID-19, the largest pool of potentially protective, pre-existing cellular immunity is derived from one of the least public epitope-specific repertoires, possibly reflecting the influence of repeated acute infections with HCoVs throughout the life of the individuals.

Beyond the comprehensive deciphering of TCR specificity reported here, we also provided a detailed picture of the complex and dynamic transcriptional landscape of the CD8+ T response to SARS-CoV-2. Importantly, we were able to demonstrate that the pre-existing SPR-B07 reactivity, observed in ~80% of unexposed subjects with HLA-B*07, was predominantly associated with a central memory-like transcriptional profile (88% of SPR-B07-reactive T cells), confirming that it originates from prior exposures. In convalescent patients, we observed a much broader distribution of SPR-B07-reactive T cells spanning every functional state at proportions ranging from 5-29% (Data file S7). This is consistent with late contraction/early memory formation described for SARS-CoV-2 in a recent study (), where cells spanned naïve, central memory, various classifications of effector memory, and terminally differentiated effector memory expressing RA (TEMRA). There was no evidence for a particularly frequent “exhausted” state among SARS-CoV-2-specific CD8+ T cells, as suggested elsewhere (, ) (acknowledging that the phenotypic state is a proxy for true reactivity testing, and that blood T cells may not fully reflect what happens in the lung). We also did not find evidence of “antigenic sin” resulting from HCoV pre-exposure () that would stifle an effective response to SARS-CoV-2-unexposed B*07 individuals. Whether HLA haplotype plays a role in the durability of the CD8+ T cell responses, especially to SARS-CoV-2 vaccines, may have impact for long-term protection across different ethnic groups and geographic regions.

Another interesting observation from this work, as noted by others (), is that even at the height of infection or shortly after viral clearance, the cumulative anti-SARS-CoV-2 CD8+ T cell response barely reached the frequency of anti-influenza memory responses and was well below the frequencies that could be achieved by CMV-specific cells in the same individuals (Fig. S6). This was notably evident in the acutely infected individuals, at a time where the contribution of cytotoxic CD8+ T cells would have been most important. We acknowledge the caveat that peripheral frequencies were measured, and some degree of sequestration in viral target tissues, such as the lung, is likely to occur in acute patients. Yet, the response seems much more muted than the robust response observed in some other viral infections (). This meager outcome was seen both for the cross-reactive “secondary responses” by memory T cells pre-primed by endemic HCoVs, as well as for the primary responses of truly SARS-CoV-2 species-specific CD8+ T cells amplified de novo. This suggests that the paucity likely does not result from a blocking of primary activation, but from a dampening of all specific CD8+ T cells. Consistent with this notion, the detection of influenza/EBV/CMV reactive cells were also lower in acute COVID-19 patients, compared to SARS-CoV-2 “naïve” individuals. It has been proposed that the lethal cytokine storm in severe COVID-19 stems from innate immune functions overcompensating for adaptive immune system failures ().

Given the widespread lymphopenia observed in acute COVID-19, we considered the possibility of latent virus reactivation with the loss of protective CMV- and EBV-specific T memory pools. While we have no direct evidence of impact on disease outcome, we do observe a significant alteration of cell state within these subsets. While CMV-reactive cells remained within, though somewhat shuffled, the same effector/memory transcriptional phenotypes between unexposed and COVID-19 cohorts (including chronic stimulation, cytotoxic terminal effector, and terminal effector memory), we observed a striking shift of EBV-specific cells from chronic stimulation and central memory into the “CD127+ memory” state in COVID-19-exposed individuals. These cells expressed moderate to high levels of many naïve (IL7R, SELL, CCR7), memory (GZMK), and effector-associated genes (NKG7, CST7, GZMA), along with markers of activation/exhaustion (TIGIT, LAG3), making them particularly interesting and difficult to ascribe to conventional phenotype labels. Recently, two transcriptionally distinct stem-like CD8+ T cell memory states were described, one of which was functionally committed to a dysfunctional lineage (). As these cell states were differentiated by many of the same markers observed in our “CD127+ memory” compartment, it would be interesting to determine to what extent these “CD127+ memory” cells, dominated by EBV-reactive pools, experience similar fates of dysfunction. We speculate that this phenotype may be a consequence of the particular inflammatory milieu of COVID-19 patients.

There are several limitations to this study. While we have investigated the CD8+ T cell response to epitopes predicted to be presented with high affinity in four common HLA alleles, the selection of HLAs and epitopes was not exhaustive. We assessed predominantly 9-mer epitopes from canonical open reading frames of a single SARS-CoV-2 variant. Subsequent studies may include a more comprehensive set of epitopes, broader coverage of HLAs, the exploration of non-canonical open reading frames, and inclusion of several SARS-CoV-2 variants. Another limitation of this study is the small sample size for specific HLA alleles and limited cell recovery for samples from acute patients. Our findings on response prevalence, public features of T cell repertoires, and T cell phenotype could be further substantiated or broadened with deeper sampling across genetically diverse populations and larger cell inputs. Related to this, there are limitations in the interpretation of response frequencies calculated in this work, especially in the cases of low cell input. The frequencies calculated are intended to provide a qualitative assessment of T cell response, allowing for comparisons across subjects, HLAs, and epitopes.

In conclusion, we leveraged a powerful single-cell technology to better elucidate the roles of HLA variation, TCR diversity, and cellular phenotypes in establishing pre-existing immunity to SARS-CoV-2. We observed the presence of a diverse and immuno-dominant nucleocapsid epitope-specific memory pool in subjects with HLA-B*07 but saw little evidence of similar reactivity in individuals with other HLA alleles. Outside of the HLA-B*07, the epitope-specific TCR repertoires observed were largely public in nature. We measured a diverse landscape of T cell phenotypes associated with SARS-CoV-2 infection, and also observed an influence on T cell repertoires reactive to persistent and latent infections with other viruses. Overall, this work provides a framework for the unified characterization of the cellular response to novel viral infections. The ability to understand the basis of cellular immunity to SARS-CoV-2 and other pathogens will provide insight for the continued assessment of immune surveillance, health security, and long-term protection from future respiratory pathogens.

MATERIALS AND METHODS

Study design. The aim of this study was to identify features of CD8+ T cell responses to SARS-CoV-2 associated with disease state and HLA genetics, including immunodominant T cell epitopes, evidence of immune recall, and shared TCR sequence motifs. We used libraries of peptide-HLA tetramers with epitopes derived from across the SARS-CoV-2 proteome presented in four HLAs with high prevalence in North America. Samples from acute and convalescent patients, with HLAs matching the tetramer libraries, were acquired as they became available and screened in several batches alongside samples from unexposed subjects. A total of 27 acute, 28 convalescent, and 23 unexposed subjects were screened providing HLA-matched analysis for 43 A*02:01, 18 A*24:02, 17 B*07:02, and 9 A*01:01 samples.

Antigen library design. Antigenic peptide libraries were designed by scoring all possible 9mer peptides derived from the entire SARS-CoV-2 proteome (NC_045512.2) using netMHC-4.0 () in the HLA-A*02:01, HLA-A*01:01, HLA-A*24:02 or HLA-B*07:02 alleles. SARS-CoV-1 peptides that had evidence of T cell positive assays, obtained from the Immune Epitope Database (www.iedb.org; ()), and that were highly homologous to their SARS-CoV2 counterparts within hamming-distance of 2 were converted to 9-mers. Additionally, SARS-CoV-2 peptides predicted to raise immunogenic responses by others were also included (, ). Finally, libraries included a set of well-defined viral epitopes from Cytomegalovirus, Epstein-Barr virus, and Influenza viruses (CEF peptide pool) that elicit T cell responses in the population at large. Antigenic peptides with 500 nM affinity or lower were then selected for inclusion (Data file S8).

Production of tetramer library pools. HLA-A*01:01, -A*02:01, -A*24:02 and HLA-B*07:02 extracellular domains were expressed in E. coli and refolded along with beta-2-microglobulin and ultraviolet (UV)-labile place-holder peptides STAPGJLEY, KILGFVFJV, VYGJVRACL and AARGJTLAM, respectively (). A C-terminal sortase recognition sequence on the HLA was modified by sortase transpeptidation (, ) with a synthetic alkynylated linker peptide, featuring an N-terminal triglycine connected to propargylglycine via a PEG linker (Genscript, Piscataway, NJ). The modified HLA monomer was then purified by size exclusion chromatography (SEC). Full-length streptavidin with an N-terminal Flag tag and a C-terminal sortase recognition sequence and 6xHisTag was prepared by expression and purification from E. coli using immobilized metal affinity chromatography and SEC. Streptavidin was modified by sortase transpeptidation with a synthetic azidylated linker peptide, featuring an N-terminal triglycine connected to picolyl azide via a PEG linker (Click Chemistry Tools, Scottsdale, AZ). HLA tetramers were produced by mixing alkynylated HLA monomers and azidylated streptavidin in 0.5 mM copper sulfate, 2.5 mM BTTAA (2-(4-((Bis((1-(tert-butyl)-1H-1,2,3-triazol-4-yl)methyl)amino)methyl)-1H-1,2,3-triazol-1-yl)acetic acid) and 5 mM ascorbic acid for up to 4 hours on ice, followed by purification of highly multimeric fractions by SEC. Individual peptide exchange reactions containing 500 nM HLA tetramer and 60 uM peptide were exposed to long-wave UV (366 nm) at a distance of 2-5 cm for 30 min at 4°C, followed by 30 min incubation at 30°C. A biotinylated oligonucleotide barcode (Integrated DNA Technologies) was added to each individual reaction followed by 30 min incubation at 4°C. Individual tetramer reactions were then pooled and concentrated using 30 kDa molecular weight cut-off centrifugal filter units (Amicon). Tetramer production was quality controlled using SEC (Fig. S1a), sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) (Fig. S1b), and UV-mediated peptide exchange by assessing binding to peptide-expanded cell lines (Fig. S2).

Patient Samples. Peripheral blood mononuclear cells (PBMCs) from COVID-19 positive donors or unexposed donors were obtained from Precision 4 Medicine (USA), the Massachusetts Consortium on Pathogen Readiness (MassCPR, Boston, USA), or CTL (USA), all under appropriate informed consent. Patients were defined COVID-19 positive based on positive SARS-CoV-2 real-time reverse transcriptase–polymerase-chain-reaction (RT-PCR) using nasopharyngeal swabs. Patient samples were characterized as “acute” if collected while the patient was hospitalized and as “convalescent” if collected after recovery or when presenting mild disease. Samples from unexposed subjects were collected prior to December, 2019. A summary of patient samples used in this study are presented in Data file S2.

Cell Staining. PBMCs were thawed, and CD8+ T cells were enriched by magnetic-activated cell sorting (MACS) using a CD8+ T Cell Isolation Kit (Miltenyi) following the manufacturers protocol. The CD8+ T cells were then stained with tetramer libraries (Data file S8), matched to subject HLAs (Data file S2) and at 1nM final concentration for each member, in the presence of 2 mg/mL salmon sperm DNA in PBS with 0.5% BSA solution for 20 min. Cells were then labeled with anti-TCR antibody-derived tag (ADT, clone IP26, Biolegend, CA, USA) for 15 min followed by washing. Tetramer bound cells were then labeled with phycoerythrin (PE) conjugated anti-DKDDDDK-Flag antibody (BioLegend, CA, USA) followed by dead cell discrimination using 7-amino-actinomycin D (7-AAD). The live, tetramer positive cells were sorted (Fig. S3) using a Sony MA900 Sorter (Sony). When necessary, sorting gates were set liberally to enable sufficient cell recovery for single-cell sequencing.

Sample multiplexing. To ensure sufficient cell loading and subsequent cDNA production in single-cell sequencing, we used sample multiplexing for several experiments. When applied, samples were independently stained with tetramer libraries, labeled using custom anti-TCR ADTs with unique 15 base pair DNA barcodes (clone IP26, BioLegend, CA, USA), and sorted. ADT-labeled, sorted samples were combined prior to encapsulation and single-cell sequencing. In several cases, an expanded T cell line (Cellero Anti-MART-1, MA, USA) was labeled with a BV785 anti-CD8 antibody (BioLegend, CA, USA), stained using a tetramer for ELAGIGILTV in A*02:01, and subsequently mixed and co-sorted alongside samples interrogated for this study. This provided confirmation of tetramer staining, guidance for gating, and verification of the multiplexing strategy (Fig. S3). The anti-MART-1 T cells (TCR sequences provided in Data file S9) were excluded from any subsequent analyses.

Single-cell Sequencing. Tetramer positive cells were counted by Nexcelom Cellometer (Lawrence, MA, USA) using AOPI stain following manufacturer’s recommended conditions. When possible, 15,000 cells were targeted for encapsulation. Single-cell encapsulations were generated utilizing 5′ v1 Gem beads from 10x Genomics (Pleasanton, CA, USA) on a 10x Chromium controller and downstream TCR, Gene Expression, and Surface marker libraries were made following manufacturer recommended conditions. All libraries were quantified on a BioRad CFX 384 (Hercules, CA, USA) using Kapa Biosystems (Wilmington, MA, USA) library quantified kits and pooled at an equimolar ratio. TCRs, Gene Expression, surface markers, and tetramer generated libraries were sequenced on Illumina (San Diego, CA, USA) NextSeq550 instruments. Sequencing data were processed using the Cell Ranger Software Suite (Version 3). Samples were demultiplexed and unique molecular identifier (UMI) counts were quantified for TCRs, tetramers, and gene expression.

Single-cell Transcriptomic Analysis. Hydrogel-based RNA-seq data were analyzed using the Cell Ranger package from 10X Genomics (v3.1.0) with the GRCh38 human expression reference (v3.0.0). Except where noted, Scanpy (v1.6.0 ()) was used to perform the subsequent single cell analyses. Any exogenous control cells identified by TCR clonotype were removed before further gene expression processing. Hydrogels that contain UMIs for less than 300 genes were excluded. Genes that were detected in less than 3 cells were also excluded from further analysis. Several additional quality control thresholds were also enforced. To remove data generated from cells likely to be damaged, upper thresholds were set for percent UMIs arising from mitochondrial genes (13%). To exclude data likely arising from multiple cells captured in a single drop, upper thresholds were set for total UMI counts based on individual distributions from each encapsulation (from 1500 to 3000 UMIs). A lower threshold of 10% was set for UMIs arising from ribosomal protein genes. Finally, an upper threshold of 5% of UMIs was set for the MALAT1 gene. Any hydrogel outside of any of the thresholds was omitted from further analysis. A total of 15,683 hydrogels were carried forward. Gene expression data were normalized to counts per 10,000 UMIs per cell (CP10K) followed by log1p transformation: ln(CP10K + 1).

Highly variable genes were identified (1,567) and scaled to have a mean of zero and unit variance. They were then provided to scanorama (v1.7, ()) to perform batch integration and dimension reduction. The data were used to generate the nearest neighbor graph which was in turn used to generate a UMAP representation that was used for Leiden clustering. The hydrogel data (not scaled to mean zero, unit variance, and before extraction of highly variable genes) were labeled with cluster membership and provided to SingleR (v1.4.0, ()) using the following references from Celldex (v1.0.0, ()): Monaco Immune Data, Database Immune Cell Expression Data, and Blueprint Encode Data. SingleR was used to annotate the clusters with their best-fit match from the cell types in the references. Clusters that yielded cell types other than types of the T Cell lineage were removed from consideration and the process was repeated starting from the batch integration step. The best-fit annotations from SingleR after the second round of clustering and the annotation was assigned as putative labels for each Leiden cluster. Further clustering of transcriptomic data was performed across the genes shown in Fig. 5 using KMeans in sklearn (v0.24) with n_clusters set to 8. As the method has a preference to assign like-sized clusters, further consolidation of two central memory clusters was performed.

In order to provide corroboration for the SingleR best-fit annotations and further evidence as to the phenotype of the clusters, gene panels representing functional categories (Naïve, Effector, Memory, Exhaustion, Proliferation) were used to score each hydrogel’s expression profiles using scanpy’s “score_genes” function () which compares the mean expression values of the target gene set against a larger set of randomly chosen genes that represent background expression levels. The gene panels for each class were: Naïve - TCF7, LEF1, CCR7; Effector - GZMB, PRF1, GNLY; Memory - AQP3, CD69, GZMK; Exhaustion - PDCD1, TIGIT, LAG3; Proliferation - MKI67, TYMS. The gene expression matrix for all hydrogels were first imputed using the MAGIC algorithm (v2.0.4, (). These functional scores were the only data generated from imputed expression values.

Scoring peptide-HLA-TCR interactions. Tetramer data analysis was performed using built-in methods of pandas (v1.2.5) and numpy (v1.20.3) in Python (v3.7.3). For each single-cell encapsulation, tetramer UMI counts (columns) were matrixed by cell (rows) and log-transformed. Duplicates of this matrix were independently Z-score transformed by row or column, and subsequently median-centered by the opposite axis (column or row), respectively (Fig. S7). For each peptide-HLA-cell interaction, this provided two scores - inter-tetramer () and inter-cell (), which were used to calculate a classifier for unique CDR3 a/b clonotypes across cells as . Classifier thresholds for positive interactions were set at 40, 36, 50, and 65 for A*02:01, B*07:02, A*24:02, and A*01:01, respectively.

Frequency Calculation. The frequency of reactive T cells in parent CD8+ T cell populations was estimated using a calculation of compounded frequency by taking the product of the fraction of reactive cells in the sorted population and the fraction of cells sorted (Fig. S8). When sample multiplexing was applied, care was taken to include only de-multiplexed cells from the corresponding sample to determine reactive cell fraction.

TCR Network Analysis. TCR motif analysis was performed using scirpy (v0.6.1) with receptor_arms = “any,” metric = “alignment,” and default cutoff of 10. Once clusters were identified, sequence alignment was performed using the pairwise2 module in Biopython (v1.78) and visualized using logomaker (v0.8).

Recombinant TCR validation. Recombinant TCRs identified from patient samples were ordered from TWIST Biosciences in the pLVX-EF1a lentiviral backbone (Takara) as a bicistronic TCRb-T2A-TCRa vector. Viral supernatants from transfected HEK 293T cells were collected 48 and 72 hours after transfection and added to the parental TCRab−/− Jurkat J76 cell like () expressing CD8 and a nuclear factor of activated T cells (NFAT)-green fluorescent protein (GFP) reporter, referred to as J76-CD8-NFAT-GFP. Recombinant TCR surface expression was confirmed through flow cytometry by staining transduced J76-CD8-NFAT-GFP cells with anti-CD3-PE (Clone UCHT1) and anti-TCRab-allophycocyanin (APC) antibodies (Clone IP26).

To assess functional activity of recombinant TCRs, J76-CD8-NFAT-GFP expressing recombinant TCRs were incubated at a 1:1 ratio with the HLA-A*02:01+and HLA-B*07:02+ HCC 1428 BL (ATCC CRL-2327) lymphoblastic cell line, with a final concentration of 0.5% dimethylsulfoxide (DMSO, vehicle) or 50 uM of cognate peptide (New England Peptide, >95% pure). Cell mixtures were incubated in the Sartorius IncuCyte at 37°C, 5% CO2 overnight and analyzed for NFAT-GFP expression measured as total integrated intensity (GCU x mm2/image) at 12 hours after assay setup. At 16 hours, cells were removed from the IncuCyte and subsequently washed and blocked with staining buffer (BD 554656), stained with anti-CD3-PE-Cy7 (Clone UCHT1) and anti-CD69-APC (Clone FN50) antibodies, and analyzed using the Intellicyt iQue Screener Plus and FlowJo v10. CD69 activity was measured as percent positive of CD3+ cells.