Identification of TCR repertoires in asymptomatic COVID-19 patients by single-cell T-cell receptor sequencing

The T cell-mediated immune responses associated with asymptomatic infection (AS) of SARS-CoV-2 remain largely unknown. The diversity of T-cell receptor (TCR) repertoire is essential for generating effective immunity against viral infections in T cell response. Here, we performed the single-cell TCR sequencing of the PBMC samples from five AS subjects, 33 symptomatic COVID-19 patients and eleven healthy controls to investigate the size and the diversity of TCR repertoire. We subsequently analyzed the TCR repertoire diversity, the V and J gene segment deference, and the dominant combination of αβ VJ gene pairing among these three study groups. Notably, we revealed significant TCR preference in the AS group, including the skewed usage of TRAV1-2-J33-TRBV6-4-J2-2 and TRAV1-2-J33-TRBV6-1-J2-3. Our findings may shed new light on understanding the immunopathogenesis of COVID-19 and help identify optimal TCRs for development of novel therapeutic strategies against SARS-CoV-2 infection.


Introduction
The ongoing pandemic of coronavirus disease 2019 (COVID-19) has posed a serious public health threat to the international communities. As of May 27, 2022, there have been 525,467,084 confirmed cases of COVID-19, including 6,285,171 deaths worldwide reported to the World Health Organization [1]. The clinical presentations of SARS-CoV-2 infections are highly variable, ranging from asymptomatic to fatal conditions [2,3]. As the number of asymptomatic individuals continues to increase worldwide, the asymptomatic infection has raised concerns among global public health authorities [4]. It is conceivable that the host immune responses may play important roles in shaping the disease outcomes. Indeed, the asymptomatic subjects have been suggested to generate timely and well-coordinated protective immune responses against SARS-CoV-2 [5,6]. Sekine et al. reported that SARS-CoV-2specific T cells were detectable in antibody-seronegative exposed family members and convalescent individuals with a history of asymptomatic and mild COVID-19 [7]. Altmann et al. observed that T cell memory is potentially more durable across asymptomatic, mild, and severe COVID-19 than B cell memory [8]. Mazzoni et al. found that asymptomatic patients displayed a significantly lower magnitude of both cell-mediated and humoral immune responses to the virus, as compared to symptomatic individuals [9]. While Reynolds et al. reported that neutralizing antibody titers were maintained irrespective of symptoms [10]. However, the immune responses, especially the T cellmediated responses and the underlying mechanisms in the asymptomatic COVID-19 patients remain to be further characterized.
T cells are important components of cellular immunity, killing virusinfected cells through T cell receptor (TCR)-mediated recognition of viral antigens [11,12]. Approximately 95 % of the TCRs are heterodimers of alpha (α) and beta (β) subunits in humans, which are shaped through the recombination (rearrangement) of variable (V), diversity (D), and joining (J) gene segments for the specificity of given antigen [13]. T cells undergo clonal expansion to generate cells with identical TCRs after activation by pathogen/antigen [14]. Therefore, unraveling the size and the diversity of TCR repertoire in the asymptomatic subjects may help understand the T cell-mediated immune responses against SARS-CoV-2 infection. Previous studies reported the TCR bias in COVID-19 patients [15][16][17]. In addition, Minervina A. et al. revealed the dynamics of memory T cells in patients with mild COVID-19 infection by Tcell receptor (TCR) sequencing [18]. Snyder T. et al. developed a method for diagnosis of SARS-CoV-2 infection based on TCR sequencing of blood samples [19]. Moreover, Wang P. et al. designed a novel T cell-based vaccine for SARS-CoV-2 variants and sarbecoviruses [20]. However, the TCR repertoires in the asymptomatic subjects (AS) COVID-19 patients remain unknown.
In this study, we performed single-cell TCR sequencing (scTCR-seq) of PBMCs from five asymptomatic individuals (AS), 33 symptomatic COVID-19 patients (SM) and eleven healthy controls (HC), and revealed the TCR repertoire diversity, the V and J gene segment preference, and the dominant combination of αβ VJ gene pairing across different disease states. Our findings may help understand the T-cell mediated immunity and identify optimal TCRs for immunotherapy against SASR-CoV-2 infection.

Study cohort and blood sample collection
This study was approved by the Ethics Committee of The First affiliated Hospital of Xi'an Jiaotong University (XJTU1AF2020LSK-015) and The Renmin Hospital of Wuhan University (WDRY2020-K130). All participants enrolled in this study offered the written informed consent by themselves or their surrogates. The case definition and classification of the study subjects followed the Guidelines of the World Health Organization (WHO) and the "Guidelines on the Diagnosis and Treatment of the Novel Coronavirus Infected Pneumonia" developed by the National Health Commission of People's Republic of China (Table S1) [21][22][23]. The study cohort included eleven healthy controls (HC, n = 11) and 38 SARS-CoV-2-infected individuals (COV, n = 38) consisting of five asymptomatic subjects (AS, n = 5) and 33 symptomatic patients (SM, n = 33) ( Table 1, Table S1). The peripheral blood was collected into the standard EDTA-K2 Vacuette Blood Collection Tubes (Jiangsu Yuli Medical Equipment Co., Ltd., China; Cat.Y09012282) and the peripheral blood mononuclear cells (PBMCs) were isolated by Ficoll density gradient centrifugation method as described [24]. PBMCs were stored in the liquid nitrogen storage tank (− 196 • C) until being used for the studies. Of note, two of the patients (C19 and C26) had blood samples collected at two different time points during this study. All the experimental procedures were completed inside a biosafety level 2 (BSL-2) laboratory at the Department of Clinical Diagnostic Laboratories, Renmin Hospital of Wuhan University. All the clinical laboratory tests were conducted by the Clinical Laboratories at the Renmin Hospital of Wuhan University, including the  routine blood tests, cytokines (IL-2, IL-4, IL-6, IL-10, AU5400, Beckman  Coulter, Inc. USA), flow cytometry analysis of lymphocytes (CD3, CD4,  CD8, CD19, CytoFLEX LX, Beckman Coulter, Inc. USA), and the SARS-CoV-2 specific IgM (Cat. 20203400769, YHLO Biotech) and SARS-CoV-2 specific IgG (Cat. 20203400770, YHLO Biotech) antibodies.

Preparation of single-cell suspensions
The frozen PBMC samples were retrieved from the liquid nitrogen storage tank, which were mixed with 10 mL of washing medium (90 % DMEM+10 % FBS) in a 15-mL polypropylene tube and centrifuged at 500g for 20 min. The supernatant was then removed (repeated this step twice). The cell pellets were resuspended with 500 μL of 1× PBS with 0.04 % BSA (A1933-25G, SIGMA), followed by adding 5 mL of 1× Red blood cell (RBC) lysis buffer (130-094-183, 10× stock solution, Miltenyi Biotech) and incubated at room temperature for 10 min to lyse the remaining red blood cells. After incubation, the cell suspension was centrifuged at 500g for 20 min at room temperature and then resus- The cell pellet was resuspended in 50 μL of 1× PBS with 0.04 % BSA. The cell viability was measured by the trypan blue exclusion method (15250-061, Invitrogen) and confirmed to be 85 % or higher. The cell number of the single cell suspension was counted using a Countess II Automated Cell Counter and the final concentration was adjusted to 700-1200 cells/μL.

Chromium 10× Genomics library construction and TCR sequencing
Cells were captured using the Chromium Single-Cell 5 ′ kit (V1) according to the manufacturer's instructions (10× Genomics), followed by cDNA amplification and library construction according to the standard protocols. The libraries were sequenced on an Illumina NovaSeq 6000 sequencing platform (paired-end multiplexing run, 150 bp) by LC-Bio Technology (Hangzhou, China) at a minimum depth of 5000 reads per cell. The sequencing results were demultiplexed and converted to FASTQ format using Illumina bcl2fastq software. To visualize the data, we further reduced the dimensionality of all cells by Seurat and used tdistributed stochastic neighbor embedding (t-SNE) to project the cells into 2D space. The steps included: (1) Using the LogNormalize method of the "Normalization" function of the Seurat software to calculate the expression level of genes; (2) the principal component analysis (PCA) was performed using the normalized expression level, within all the PCs, the top 10 PCs were used to conduct the clustering and t-SNE analysis; (3) using weighted shared nearest neighbor (SNN) graph-based clustering method to find clusters. The marker genes for each cluster were identified with the "bimod" (likelihood-ratio test) with default parameters via the FindAllMarkers function in Seurat. This selects marker genes that expressed in >10 % of the cells in a cluster and the average log (fold change) of >0. 26. To further avoid interference of putative multiplets (where more than one cell was loaded into a given well on an array), cells in a defined cluster that had high expression of more than one cell type canonical marker gene were filtered to ensure the data quality. As a result, a total of 119,799 cells were used for the final analysis in this study. The nine cell types were integrated for further subclustering. After integration, genes were scaled to unit variance. Scaling, principal component analysis and clustering were performed as described above.
Sample demultiplexing, barcode processing and single-cell 5 ′ gene counting were completed by using the Cell Ranger pipeline (https ://support.10xgenomics.com/single-cell geneexpression/software/pipe lines/latest/what-is-cell-ranger, version 3.1.0). The scRNA-seq data were aligned to Ensembl genome GRCh38 reference genome. The Cell Ranger Seurat (version 3.1.1) was used for dimensional reduction, clustering, and analysis of scRNA-seq data. Approximately 200,000 cells passed the quality control threshold: genes that expressed in less than one cell were removed, number of genes expressed per cell >500 as low cut-off, UMI counts <500, and the percentage of mitochondrial-DNA derived gene-expression < 25 %.
In detail, the TCR Barcoded, full-length TCR V(D)J segments were enriched from the amplified 5 ′ cDNA libraries using the Chromium Single-Cell V(D)J Enrichment kit according to the manufacturer's protocol (10× Genomics). The TCR transcripts were enriched in separate reactions from the same amplified cDNA material if both T cells were expected to be present in the partitioned cell population.
Cell Ranger pipeline (https://support.10xgenomics.com/single-ce ll-vdj) that processes Chromium single cell 5 ′ RNA-seq output for V(D) J was employed to assemble, quantify, and annotate paired V(D)J transcript sequences. The workflow of Cell Ranger started by demultiplexing the Illumina sequencer's base call files (BCLs) for each flowcell directory into FASTQ files. Cell calling was performed independently of V(D)J read filtering and assembly. Clonal expansion cells shared the same paired receptor sequences (presumably derived from the same progenitor cell). Cell barcodes were grouped into clonotypes if they shared the same set of productive CDR3 nucleotide sequences. The repertoire information of all samples was exhibited based on the following features: CDR3 abundance, CDR3 base length, variable (V) and joining (J) segment usage (the frequency of associated reads for each V/J gene presented in the samples), the length of V/J gene in CDR3 region, and V-J gene paired frequency in CDR3 junctions. We employed the hierarchical cluster (https://stat.ethz.ch/R-manual/Rdevel/library/ stats/html/hclust.html) and Kruskal's Non-metric Multidimensional Scaling (isoMDS, https://cran.r-project.org/web/packages/MASS /MASS.pdf) to estimate the Immune Repertoires Similarity. To avoid batch effects, we followed the standard operational procedures (SOP) recommended by the manufacturer (10× Genomics). In particular, all the datasets in this study were generated by the same operators at the same laboratories using the same experimental protocols and sequencing platforms.

Statistical analysis
All the statistical analyses were performed using SPSS (Statistical Package for the Social Sciences) (version 18.0 software, SPSS Inc.) unless otherwise stated. The clonal diversity of each sample was calculated using Shannon's entropy, which took into consideration of both abundance and the frequency distribution of each CDR3 sequence [25]. Entropy was calculated by summing the product of the frequency of each clone with the log2 of that frequency via CalcDiversityStats function of VDJTools, which was calculated by several re-sampling steps (downsampling to the size of the smallest dataset) [26]. The normalized entropy of clonotype frequency distribution was calculated by dividing by log [number of clonotypes]. The higher Shannon's entropy, the more diverse the distribution of CDR3 clones. The categorical variables were described as frequency rates and percentages, whereas continuous variables were described using mean or median with interquartile range (IQR) values. For analysis of two groups, independent group t-test was used when the data were normally distributed; otherwise, the Mann-Whitney test was used. For analysis of three and more groups, Oneway Anova test was performed when the data were normally distributed; otherwise, Kruskal-Wallis test was used. In addition, Bonferroni correction was used for the multiple comparisons. P < 0.05 was considered as statistically significant.

Data availability
The raw sequence data reported in this paper have been deposited in GEO, under accession code GSE165080 and GSE169441 are publicly accessible at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi? acc=GSE165080 [27] and https://www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?acc=GSE169441. Other supporting raw data are available from the corresponding author upon reasonable request. Source data are provided with this paper.

Identification of eleven T cell subtypes by single-cell RNA sequencing (scRNA-seq)
We recruited eleven healthy controls (HC), five asymptomatic subjects (AS) and 33 symptomatic patients (SM) in this stud and summarized the demographic and clinical features of the study subject (Fig. 1,  Tables 1, S1). Of note, the asymptomatic subjects were tested positive for the nucleic acids of SASR-CoV-2 by RT-PCR but were lack of clinical symptoms and abnormal chest CT images. On the other hand, all 33 symptomatic patients were positive for the viral nucleic acids and presented clinical symptoms and abnormal chest CT images. In addition, we revealed that the CD3 + , CD4 + , CD8 + , and CD19 + cell counts were significantly higher in AS than SM patients (1529. 25  To profile the immune cell landscape of COVID-19 patients with different clinical presentations, we performed deep single-cell RNA sequencing on PBMCs from the healthy subjects and COVID-19 patients [27]. We detected eleven T cell subclusters, including activated effector T, double negative T (DNT), effector memory T, γδT, mucosal associated invariant T (MAIT), memory T, naive T, NKT, pro-NKT, T helper cell 17 (Th17), and regulatory T (Treg) cells (Fig. 1A). We also showed the representative marker genes for the corresponding T cell subsets, such as activated effector T (GZMH), DNT (DDX17), effector memory T (GZMK), γδT (TRDV2), MAIT (TRAV1-2), memory T cells (DDIT4), naive T (CCR7), NKT (KLRD1), pro-NKT (MK167), Th17 (CRIP2), and Treg (FOXP3) (Fig. 1B, Table S2). The ridge plot shows the distribution identity of these marker genes in the eleven T cell clusters (Fig. S1B). We further investigated TCR repertoire on T cells from the PBMCs among AS, SM, and HC by single-cell TCR sequencing (scTCR-seq). Of note, except some of the CD4 − CD8 − T cells, >90 % of cells in all subsets exhibited TCR (Fig. S1A).

Distinct T cell clonal expansions associated with AS and SM patients
In the cell-mediated immune responses, activated T cells specifically recognized the amino acid fragments of antigen presented by antigen presenting cells (APC) through TCRs, by which the immune clearance was achieved and distinctive TCR repertoire was greatly enriched. Therefore, it is important to understand the difference in TCR repertoires between asymptomatic and symptomatic patients of COVID-19 for the prognosis and optimizing the treatments. To illustrate the difference between clonal expansion and TCR bias when T cells were activated by SARS-CoV-2, we performed scTCR-seq with complete TCR α and β chains. We characterized the T cell clonal expansion status in the HC, AS and SM, respectively. We found that more clonal expansion occurred in MAIT cell cluster of AS and SM compared to HC (Figs. 2A, S1C). As expected, most T cell clones were unique (Fig. 2B, Table S4). We then calculated the relative proportion of the clonal T cells in HC, AS and SM groups, respectively. Our data suggested that there was no significant difference in the relative proportion of the clonal T cells among the three groups (Fig. 2C). However, the number of the clonotypes in SM (ranged from 3 to 191) and the clonal diversity were significantly different from the HC group. Clonal diversity in each sample was calculated using Shannon's entropy and considered both the CDR3 abundance and the frequency distribution of CDR3 sequence (Fig. 2D, E, Table S3). To clarify whether different TCR repertoire in each group was due to the different T cell subsets, we calculated the proportion of each T cell subsets in the mean level of HC, AS, SM, respectively (Fig. S1D). Overall, there was no significant difference of T cell subsets in all groups and the differences of TCR repertoire in each group was not due to the different T cell subsets. Finally, the Venn diagram showed that there were 42,198 TCRs in SM patients, 5248 in AS and 12,424 in HC, respectively. Interestingly, there were 1393 shared TCRs between SM and HC, but only two shared TCRs between SM and AS and three shared TCRs among SM, AS, and HC (Fig. 2F).

A number of significant variable and joining genes of the TCR αβ chain were correlated with disease symptoms
CDR3 is the focus of the immune repertoire because it is encoded by three genes: variable (V), diversity (D), and joining (J) genes, which leads to the highest degree of diversity in CDR3 [28]. Each immune protein in the immune repertoire has a wide variety of subtypes and plays a vital role in health. Therefore, the more subtypes of immune proteins, the more effective they are against pathogens. Here, we investigated the bias of V and J gene segments in AS and SM, and discovered that the frequency of several significant V and J genes in the TCR αβ chain were increased in AS, including TRAJ16 and TRBJ2-1, whereas the frequency of TRAV2, TRAJ8, TRAJ40, TRBV3-1 and TRBV5-1 were the highest in SM, and TRBJ2-5 was highest in HC (Fig. 3A, Table S5). Moreover, TRAV2, TRAJ8, TRAJ40, TRBV3-1 and TRBV5-1 were found to be correlated with the disease symptoms (Fig. 3B). Taken together, these results revealed the TCR V J bias in AS and SM, suggesting that there was a possible association between TCR V and J gene segments with the symptoms of COVID-19.

Discussion
In this study, we first reviewed and analyzed the clinical laboratory data and found that the CD3 + , CD4 +, and CD8 + T cell counts were significantly lower in SM compared with AS and normal reference values, suggesting that T cell exhaustion might be associated with the disease symptoms [29]. We then conducted the scTCR-seq and revealed >90 % of cells in all subsets with TCR information, demonstrating that direct sequence-based immune profiling is a powerful approach for understanding the dynamics of the immune responses [30].  We subsequently examined the T cell clonal expansion status in HC, AS and SM, respectively. We found that more clonal expansion occurred in the MAIT cluster from AS and SM compared to HC. Previous studies have shown that MAIT cells were associated with the mucosa immunity to a variety of microbial infections [31][32][33], suggesting that these clonal MAIT cells may play an important role in the immune responses to SARS-CoV-2 infection in AS and SM patients. The findings that clone patterns of TCRs were reused and the clonotype diversity was significantly different in AS and SM may help us to better understand the TCR repertoires in AS subjects and SM patients. Next, we also investigated the bias of the V and J gene segments in AS and SM and found the frequency of TRAJ16 increased in AS compared with HC, which was  consistent with the previous report in the graft-versus-host disease (GVHD) [34]. Interestingly, GVHD and COVID-19 shared similar pathological consequences [35]. In addition, TRBJ2-1 was highly expressed in AS compared with SM. Since SARS-CoV-2 infection was shown to affect the cardiovascular system, and previous studies have demonstrated the importance of TRBJ2-1 in HLA-I restricted tumor-associated peptide and acute myocardial infarction [36,37], suggesting that TRBJ2-1 might be associated with the cardiovascular disorders in COVID-19 patients [38,39]. TRBJ2-5 was the highest in HC subjects whereas TRAV2, TRAJ8, TRAJ40, TRBV3-1 and TRBV5-1 were the highest in SM. These skewed V and J gene segments might be associated with the symptomatic COVID-19. Finally, we uncovered the common and unique the V-J pairing of α and β chain rearrangements in AS, SM and HC. In particular, we identified the most frequent αβ VJ pairs in AS subjects, including TRAV1-2-J33-TRBV6-4-J2-2, TRAV1-2-J33-TRBV6-4-J2-3, TRAV1-2-J33-TRBV6-1-J2-3, which may contribute to the antiviral immunity and symptom-free state in the AS group. Of note, McGregor et al. reported the TRAV1-2-TRAJ33 TCR-α rearrangement in peripheral T-cell lymphomas derived from innate-like T-cells [40]. The preference of αβ VJ pairs in AS may have important implications for development of novel therapeutic strategies and vaccine candidates against SARS-CoV-2 infection. Due to the frequent mutations of SARS-CoV-2, T cell-based vaccines will receive more and more attention. Indeed, several TCRbased vaccine candidates have been developed, highlighting the significance of the TCR-related studies [20,41].
Some of the limitations of this study include relatively small sample size, variations in the time-points for the sample collections and relevant statistical weakness. Another limitation was lack of experimental validations on some of the important skewed TCR genes and its TRA/B rearrangements in the peripheral blood samples to rule out the random differences in pairings due to diversity in genes in the population or influences from other individual differences such as viral infections, chronic illness, differences in HLA alleles.
The complexity and dynamics of the T-cell repertoires may render huge challenges for study of the immunopathogenesis of COVID-19. However, the scTCR-seq technology allowed us to conduct unbiased research to understand the contributions of T cell-mediated immune responses to.
SARS-CoV-2 infection. Taken together, this study revealed the basic features of TCR repertoires in response to SARS-CoV-2 infection and provided novel insights into the mechanisms of T-cell immunity, which may help to identify optimal TCRs for development of novel therapeutic strategies and vaccine candidates against SARS-CoV-2 infection.

Funding
This study is supported in part by the Science and Technology Department, Shaanxi Province (Grant No. 2020ZDXM2-SF-02) (CZ and BS) and the operational funds from The First Affiliated Hospital of Xi'an Jiaotong University (CZ and BS).

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Author Weikang Mao is employed by the LC-BIO TECHNOLOGIES (HANGZHOU) CO., LTD., China. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.