Single-cell mapping of γδ T cell activation following SARS-CoV-2 vaccination

To investigate the effect of SARS-CoV-2 vaccine exposure on γδ T cells, we performed a longitudinal scRNA-seq study on 6 SARS-CoV-2–naïve volunteers vaccinated with two doses of the Pfizer-BioNTech BNT162b2 mRNA vaccine (Fig. 1a and Supplementary Fig. 1a). PB samples were collected 1 day before (P0), 3- and 17-days after prime vaccination (P1 and P2, respectively), and 3-days and 3-months following vaccine-boost (P3 and P4, respectively). γδ T cells were identified by their expression of CD3 (CD3E) and TCR canonical δ constant (TRDC) and variable δ (TRDV) region-encoding segment (Supplementary Fig. 1b). Initially identified clusters of γδ T cells were poorly-resolved, necessitating further re-clustering to obtain purified and separable groups of γδ T cells (Method section). In total, we refined 12640 γδ T cells that Uniform Manifold Approximation and Projection (UMAP) analysis resolved into 10 different γδ Τ cell clusters (c0-c9) as visualized by the top 10 differentially expressed genes (DEGs) (Fig. 1b, c, and Supplementary Fig. 1c). Specifically, we identified clusters 3 and 4 (c3, c4) as Vδ1 T cells expressing TRDV1 segments; instead, clusters c0-c2 and c5-c9, were identified as Vδ2 T cells as those expressing TRDV2. The relative frequencies of each cluster in the 6 samples are represented in Fig. 1d. As expected31, Vγ9 emerged as the main Vγ chain (TRGV) associated with Vδ2-enriched clusters (Fig. 1e). On the other hand, previously defined as adaptive-like Vγ9negVδ2 T cells32, were enriched in c0 and c2. In fact, cells in c0 and c2 present high expressions of Vγ2 and Vγ4, respectively. These cells are rarely present in adults PB32, and were mainly presented in one of the analyzed individuals. Moreover, two main Vγ3 and Vγ2 chains were linked to Vδ1-enriched T cell clusters.

Fig. 1: Activation of γδ T cells following SARS-CoV-2 vaccination.

a Schematic overview of the experimental design. b UMAP clustering projection of the integrated PB γδ T cells from all subjects (s01-s06; a total of 12640 γδ T cells). c The dot plot (left panel) and feature plots (right panel) show the expression of canonical markers used for the annotation of γδ T cell subtypes. d The bar plot shows the frequency (%) of cell cluster (c0-c9) distribution across all subjects. Cell numbers were normalized to the total number of cells per subject loaded independently from the time point. e The violin plots show the expression of the main TCR γ chains (TRGV) associated with either δ1 (TRDV1) or δ2 (TRDV2) chains per cluster. f The bar plot shows the number of identified DEGs across different time points (P1-P4) after vaccination, compared to pre-vaccine P0. For all figures, DEGs were defined as follows: (i) absolute value of average |Log2- FC| ≥ 0.25 for up- and down-regulated genes; (ii) adjusted P values (adj.p) ≤ 0.05, and (iii) detected not less than 10% of cells (min.pct ≥ 10%).

The relative abundance of Vδ1 and Vδ2 T cell clusters in each individual, and at different time points did not display substantial changes (Supplementary Fig. 1d). However, significant variations were observed in their transcriptome profiles, in terms of kinetics and amplitude of DEGs, before and after immunization. Specifically, minor changes were detected in the transcriptome of Vδ1 and Vδ2 T cells early after the prime-vaccine at P1 and P2. Instead, their transcriptional profile substantially changed after the vaccine-boost, suggesting that the first vaccine dose primed γδ T cells into a “trained” status that the boost tuned to a higher responsiveness (Fig. 1f). We then investigated the features and dynamic changes in individual γδ T cell clusters in relation to DEGs and to the underlying transcriptional programs (Reactome). Heterogeneous responses were observed across different Vδ2 T cell clusters (Fig. 2a). Clusters c1 and c6 exhibited early induction of CD69 at P1, and various members of the AP-1 transcription factor (TF), particularly prominent in c1 (Supplementary Fig. 2a). Clusters c1 and c6 were also characterized by a higher response after the vaccine booster, sharing different activated profiles, including AP-1 TFs (i.e, JUN, JUNB, JUND, FOS, FOSB, FOSL2) and genes related to the NF-κB activities (i.e., NFKBIA, NFKB1, TNFAIP3, MAP3K8) (Fig. 2b). In fact, pathway enrichment analysis based on all DEGs calculated compared to P0 resulted in a notable transcriptional reprogramming marked in c1, and involving AP-1 and NF-κB TFs, interleukins (ILs), interferons (IFNs), and TNF pathways, and the TCR-complex signaling engagement (Fig. 2c). Cascades of Toll-like receptors (TLRs) and pathways regulating the response to SARS-CoV-2 infection show progressive activation in c1 and c6. Importantly, the activation pattern in c1 after the second vaccine dose indicated a proliferation signature, reprogramming of metabolism and the memory T cell preservation (TXNIP, GATA3, and IL7R) (Fig. 2c, d)33,34,35,36.

Fig. 2: Identity of the specific γδ T cell cluster activation following SARS-CoV-2 vaccination.

a The bar plot shows the number of up- and down- modulated DEGs at different post-vaccine time points (P1-P4) calculated for each γδ T cell cluster, compared to the pre-vaccine P0. b The volcano plots highlight DEGs of the specific Vδ2 T cell clusters c1 and c6, resulting in the highest number of DEGs at P3. Green boxes indicate common genes expressed in clusters c1 and c6. Red dots indicate genes up-regulated in P3 vs P0, and blue dots indicate up-regulated genes in P0 vs P3. c The dot plot shows a selection of significantly enriched pathways with FDR values <0.05, identified among DEGs at each time point (P1-P4) for specific Vδ2 T cell clusters (c0-c2, c6) using the Reactome pathway browser. Dots are colored by FDR values and sized by the number of DEGs enriched in each pathway. d The ridge plots show the expression levels (x-axis, log-UMI) and the frequency of cells (y-axis) of three DEGs found in cluster c1 at P4 compared to pre-vaccine P0. Null gene expression cells were excluded from the analysis. The dotted line indicates the gene expression level corresponding to the highest cells’ frequency at the baseline (P0). e The volcano plots highlight DEGs of the Vδ2 T cell cluster c0 at P3 compared to P0. Red dots indicate genes up-regulated in P3 vs P0, and blue indicate genes up-regulated in P0 vs P3.

Vδ2 T cell activation in c0 and c2 were related to their cytotoxic potential, as evidenced by increased DEGs such as FCGR3A, PRF1, ACTG1, S100A4, and CX3CR1, along with the down-modulation of AP-1 TFs, distinct from c1 and c6 (Fig. 2c, e). The remaining clusters of Vδ1 and Vδ2 T cells demonstrated reduced biological pathway activation in response to vaccine stimuli (Supplementary Fig. 2b).

Based on the heterogeneous response found among different clusters of γδ T cells, we reasoned that γδ T cells might exhibit varying responsiveness to SARS-CoV-2 vaccine strengthened upon immunization. Moreover, we identify AP-1 TFs as potential drivers that contribute to the increased Vδ2 T cell response at P3 after the second vaccine dose.

Transcriptional changes in Vδ2 T cells following SARS-CoV-2 vaccination

We found that the main source of variation across all clusters can be imputed to the differentiation phenotypes of γδ T cells regardless of the vaccination (Fig. 3a). Indeed, the cells in c8 with TN profile express higher levels of CCR7, LEF1, FYB1, SELL, and COTL137,38. Cells in c1 show TN/TCM profile that share several markers including IL7R, TCF7, and CD27. While cells expressing cytotoxic-related genes (i.e., GZMA, GZMB, and PRF1) represent γδ T cells that exhibit TEM/TEMRA patterns. We observed that TCM Vδ2 T cell clusters, enclosed in c5-c7, show a gradual transition to the TEM profile. Compared to more heterogeneous Vδ2 T cells, Vδ1 T cells aggregate in two clusters with TEM/TEMRA cytotoxic outlines that mainly differ in Vγ chain usage, expression of KLRC2 (NKG2C), and TFs expression such as ZNF683 (Supplementary Fig. 1c), thus indicating their diverse functional commitment.

Fig. 3: SARS-CoV-2 vaccination shapes the effectorness of Vδ2 T cells.

a The dot plot shows the expression of selected genes for each γδ T cell cluster at all the time points (P0-P4), together. Dots are colored by the average expression of each gene scaled across all clusters and sized by the percentage of cells within a cluster (min.pct ≥ 10%). Clusters are ordered according to the hierarchical clustering shown by the dendrogram. b Pseudotime trajectory of Vδ2 T cells, where each cell is colored by its cluster identity (c0-c2, c5-c9) (left panel) or its pseudotime value (right panel). c Effectorness gradient mapped on the pseudotime trajectory of Vδ2 T cells, with each cell colored by its effectorness value. d Heatmap of selected genes variable along the pseudotime trajectory (from Monocle) calculated before vaccine at P0 (left panel) or after all (P1 + P2 + P3 + P4) vaccination time points (right panel). The x-axis represents cells ordered by pseudotime (from left to right), and different colors correspond to the scaled (Z-scored) expression of each gene in each cell. e The box plot shows effectorness values at different time points (P0-P4), with the median represented by a line and the mean indicated by a “+” across the boxes. The vertical bars show the min-max range distribution. For the statistical analysis, the unpaired nonparametric Dunn’s test was used to perform multiple comparisons vs the control P0, statistically significant comparisons are represented as P values (*): ****P < 0.0001.

To further investigate relationships between Vδ2 T cell clusters, we applied pseudotime analysis on Vδ2 T cells including all the time points (Fig. 3b and Supplementary Fig. 3). The trajectory of Vδ2 T cells started from highly enriched TN/TCM markers (e.g., COTL1, GZMK, LEF1, CD27, and IL7R) and finished with the markers of highly differentiated cells expressing B3GAT1 (CD57), NCAM1 (CD56), FCGR3A (CD16), and CX3CR1, cytotoxic molecules (GZMB, GZMH, PRF1, GNLY), cytokines and chemokines (IFNG, CCL3-4, CCL4L2). Based on the pseudotime analysis, we observed that Vδ2 T cell clusters do not comprise discrete subpopulations but instead represent interrelated multiple differentiation stages with a continuous advancement from TN and TCM profiles, toward cells with TEM phenotype (Fig. 3c). We reasoned that the observed Vδ2 T cell continuum reflects the potential to initiate a rapid effector response upon stimulation. We term this property “effectorness”, based on what has been observed for CD4+ T cells39. To assess if the vaccine influences impacts the effectorness of Vδ2 T cells, we generated pre- and post-vaccine pseudotime trajectories (Fig. 3d). To define the effectorness score, we scaled the pseudotime values for each trajectory to the range [0; 1] (Method section). In both trajectories, we found equivalent routes defined by initial genes corresponding to the TN profile [0 = TN] and terminal genes defining TEM cells [1 = TEM], indicating that the effectorness gradient is detectable before and after vaccination. We observed that the vaccine increases the effectorness values of cells with TN/TCM phenotype. Indeed, cells with a TN/TCM profile increased their effectorness by shifting to the right within the effector range after vaccination. Notably, AP-1 family genes, including FOSL2, JUN, JUNB, FOS, and FOSB are the main TFs associated with the increased effectorness. To support the notion that Vδ2 T cell effectorness increases upon vaccination, we calculated the cell distribution at different time points based on cell effectorness (Fig. 3e). We found an increased number of cells with an enhanced effectorness from P1 to P3, thus suggesting that the effectorness of Vδ2 T cells is a consequent of their vaccination history.

As for Vδ1 T cells, their reduced pathway activation, compared to Vδ2 T cells, hinders a similar effectorness annotation. These findings imply that the distinct TCR repertoire-phenotype states of Vδ1 and Vδ2 T cells before vaccination are associated with their differing responses to the vaccine.

SARS-CoV-2 vaccinations lead to increased effector response and memory phenotype of Vδ2 T cells

Our results indicate that the vaccine shapes the effectorness of Vδ2 T cells. To further characterize the role of SARS-CoV-2 vaccination in shaping the effectorness of Vδ2 T cells at different time points, we used a multiple linear regression (MLR) model to predict gene expression with effectorness, vaccination, and their reciprocal interaction as independent variables (Methods section)39. With this approach, we found 179 genes whose expressions were significantly modulated by effectorness and vaccination acting jointly due to their cross-reactive effect (Supplementary Table 1). Among them, we found genes with a high impact on the effector response such as IFNG, TNF, and XCL2, and several genes crucial for durable TCM phenotype preservation including IL7R, TCF7, GATA3, CXCR3, and EOMES34,35,40,41,42. Profiling the expression of these genes over different time points revealed their progressive correlation with the vaccine stimulations (Fig. 4). For example, genes regulating effector response (IFNG, TNF, and XCL2) increased (as indicated by the angular coefficient ‘m’) after the first dose of the vaccine; however, their magnitude peaked after the second dose. Similar changes have been observed for the memory-associated genes (IL7R, TCF7, GATA3, EOMES, CXCR3) that, moreover, resulted in a more durable memory phenotype detected in Vδ2 T cells over 3 months after vaccination. These data clearly show the critical role of the first vaccination that elicits enhanced immune response and transcriptional changes of Vδ2 T cells after the second dose.

Fig. 4: The interaction between effectorness and vaccination regulates Vδ2 T cell response.

Plots of gene expression (y-axis, log-normalized UMIs) as a function of effectorness (x-axis), stratified by time points (P0-P4). Examples of significant cross-reactive genes, selected for their activation (top panels), effector (middle panels), and memory-associated phenotype (bottom panels) of Vδ2 T cells. Each dot represents a single cell with a gene expression value >0. The “m” value above each graph represents the coefficients of linear regression. Significance was calculated using a t-test, and only significant P values (*) are shown (P < 0.05).

Overall, the transcriptomic status of Vδ2 T cells during SARS-CoV-2 vaccination is influenced by the mutual interaction between effectorness and the number of vaccinations, leading to an enhancement of their effector potential associated with memory profile.

Adaptation of γδ-TCR repertoire to SARS-CoV-2 vaccination

To gain a deeper understanding of the γδ T cell response to the vaccine, we linked our scRNA-seq data with the individual γδ-TCR clones generated as described in the Methods section. Among TRD or TRG sequences, detected in all barcoded cells, we found TRG and TRD pairs in about 40% of all single γδ T cell transcriptomes. For further analysis, we considered only the paired γδ TCRs. The diversity and distribution of paired DV and GV genes are shown in Fig. 5a, b, and Supplementary Fig. 4a. As expected, GV9-DV2 pairs represented the most abundant combination among all Vδ2 T cell clusters in all donors. On the other hand, DV1 mainly paired with GV2, GV9, and GV8. We then investigated the diversity of γδ-TCR clonotypes across all cell clusters, defined by UMAP, by calculating the Shannon index, which ranges from 1 (polyclonal) to 0 (monoclonal) (Fig. 5c). The two Vδ1 T cell-enriched clusters, c3-4, in concordance with their high effector phenotypes, show more oligoclonal features with a relatively low Shannon index. We observed that the gradual clonal expansion aligned with the progressive differentiation states of Vδ2 T cells, accompanied by increasing effectorness (Fig. 5c, d). Moreover, clonal Vδ2 T cell expansion correlates with the persistence of discrete subpopulations with heterogeneous phenotypes. Indeed, approximately 15% of the clonal composition of the most expanded (≥50 cells) DV2 clones overlaps with the clusters displaying TCM profile (Fig. 5e and Supplementary Fig. 4b).

Fig. 5: Single cell γδ-TCR repertoire analysis.

a Chain pairing of different GVs and DVs are displayed as chord diagrams, where ribbons connecting chains are proportional to the number of paired GV/DV chains. b The bar plot shows the frequency of the most representative paired GV/DV chains for each cluster (c0-c9). The total numbers of paired GV/DV chains per cluster are indicated at the top of each bar, and the most abundant pairs are color-coded. c γδ-TCR repertoire diversities were estimated by the normalized Shannon Index, ranging from 0 (completely monoclonal) to 1 (completely polyclonal). The bar plot shows the Shannon Index of each cluster (c0-c9). d Plot of the normalized Shannon Index (y-axis) as a function of the effectorness gradient (x-axis) of Vδ2 T cells, calculated using a sliding window approach (Methods). e The chord diagram (left panel) shows the overlap of highly expanded clones (≥50 cells) across different Vδ2 T cell clusters (c0-c2, c5-c9). Spatial visualization on the UMAP (right panel) of all identified GV/DV paired clones Vδ2 T cells, scaled from unique clonotypes to highly expanded clones (≥50).

Overall, these data indicate that Vδ2 T cell clonal expansion is accompanied by increased effectorness and the persistence of TCM phenotype.

Building upon this evidence, we tracked the fate of individual γδ-TCR clonotypes concerning the vaccination timeline. Among the total analyzed clonotypes, about 1% overlapped between more than one subject (Fig. 6a and Supplementary Fig. 5a). On the other hand, longitudinal tracking of the top 20 abundant clones in all analyzed time points revealed that the γδ-TCR repertoires in most subjects differed from those before vaccination (Fig. 6b).

Fig. 6: SARS-CoV-2 vaccine shapes γδ-TCR repertoire.

a The Venn diagram shows the counts (colored box) of overlapped paired γδ-TCR clonotypes among different subjects (s01-s06); the numbers represent subject-unique or shared γδ-TCR clonotypes among different individuals. b Longitudinal tracking of the 20 top-most expanded γδ-TCR clones for selected subjects (s01-s04), stratified for different time points (P0-P4). Each stratum represents a unique γδ-TCR clonotype highlighted by a different color. The colored bands between columns represent shared clones among time points.

A key question that has emerged was the timescale for the formation and persistence of the γδ-TCR repertoire following vaccination. To address this question, we identified all clonotypes numerically expanded during P1-P4 in relation to their basal level at P0, and which lasted 3 months after the recall of the vaccine (P4). In total, we identified 59 clones that met these criteria (Fig. 7a). The numerical cell-size distribution among expanded clones showed that the first vaccine dose prevalently induced the expansion of relatively small-sized clones (2–3 cells). In contrast, the booster prompted the expansion of larger 4–10 cell-size clones, which further increased over time (Fig. 7a and Supplementary Fig. 6a). The expansion of small cell-size clones was in line with a study on αβ T lymphocyte showing that mRNA-based SARS-CoV-2 vaccination led to the expansion of small cell-size clones (2–3 cells) compared to natural infection43. Importantly, this specific clonotype-tracking revealed that the first dose resulted in a limited clonal expansion that declined after 2 weeks. Instead, the second dose induced a higher clonal expansion which peaked 3 months later (Fig. 7b). While vaccine-preexisting clones could also expand, the prevalence (60%) of the expanded clones originated from unique clonotypes. We found the clonal expansion in 80% of analyzed individuals, including both females and males with ages ranging from 26 to 51 years old, and with both HCMV-seropositive and HCMV-seronegative status (Fig. 7c). Among the expanded clones, there was an enrichment (17% vs 1%) of public clonotypes shared by more individuals, with expansion occurring under subject-specific circumstances. Comparing the DV/GV gene usage in expanded clones, we found dominancy (76%) of GV9/DV2 paired clones that, as expected, show a high phylogenetic relationship compared to other DV1 and DV3 clones (Supplementary Fig. 6b). Importantly, the dissimilarity among individually expanded GV9/DV2 clones also resulted in significant genetic distances based on amino acid sequence homology, indicating the expansion of unique GV9/DV2 variants (Supplementary Fig. 6c). Additionally, vaccination was associated with an increase in the length of the complementary-determining region 3 (CDR3), particularly for CDR3δ (Supplementary Fig. 6d).

Fig. 7: Expansion of γδ TCR clones upon SARS-CoV-2 vaccination.

a UMAP visualization (left panel) of all identified and expanded γδ TCR clones with a cell size ≥2. Each color represents a different clonotype. The pie charts (right panel) show the percentage frequency distribution of the expanded γδ TCR clones grouped by cell sizes across all the time points (P0-P4). b The two alluvion plots show the longitudinal tracking of the expanded γδ TCR clones, grouped by the cell size, ≥3 - ≤10 (left panel) and >10 (right panel). Each stratum represents a unique γδ TCR clonotype highlighted by a different color. The time point distribution (P0-P4) and sequence of each clonotype are reported below. c The nested pie chart shows the frequencies (%) of subjects who responded to the vaccine with clonal expansion among all the analyzed individuals (inner gray chart). The distribution of HCMV seropositive status, age, and sex among responsive subjects is represented in the colored charts, progressing from the innermost to the outermost circles, which are colored in orange, green, and violet, respectively.

Hence, the γδ-TCR repertoire displayed a high clone-specific sensitivity to SARS-CoV-2 vaccination.

Multiple signals converge on memory Vδ2 T cell development following the SARS-CoV-2 vaccination

As we mapped the expanded clones onto the UMAP, we noted their overlap with the T cell clusters responsive to the vaccine. We then integrated signaling networks of activation and proliferation with the γδ-TCR repertoire changes induced by vaccination. Firstly we applied the known signature of cycling γδ T cells to measure proliferative changes in the expanded γδ-TCR clones (Fig. 8a)37. Importantly, we found the peak of proliferation occurring after the second vaccine dose. This was also supported by higher levels of G2M phase cell cycle genes at P3 (Supplementary Fig. 7a and Supplementary Table 2). Subsequently, we linked the transcriptional profile of γδ T cell activation to the TCM signature and vaccine stimulation in the identified clones. Following each vaccine dose, we observed an increased expression of genes like CD69 and AP-1 TFs (FOS and JUN) (Fig. 8b), confirming the activation signature of the most responsive clusters found in Fig. 2b. This clonal activation correlated with the acquisition of effector mediators (FNG, TNF, XCL2) and the previously identified (Fig. 4) TCM signature (IL7R, TCF7, GATA3, CXCR3, EOMES) that persisted at P4 (Fig. 8c). Importantly, by applying this memory signature to calculate its enrichment score in different γδ T cell clusters, we confirmed in vaccine-responsive clones the prevalence of TCM immune phenotypes (Fig. 8d).

Fig. 8: Establishment of memory-like response of Vδ2 T cells following repeated SARS-CoV-2 vaccination and peptide in vitro stimulation.

a The box plot shows the median values, represented by a line across the boxes, and vertical bars of the min-max range distribution of the cycling gene score calculated in the expanded γδ TCR clones at different time points (P0-P4). For the statistical analysis, the unpaired parametric ANOVA test was used, and only significant p values (*) are shown: *P < 0.05 and ****P < 0.0001. b, c The ridge plots show the expression levels (x-axis, UMI-log) and the frequency of cells (y-axis) of selected genes on the expanded γδ TCR clones at different time points (P0-P4), divided for activation, effector and memory signature. The dotted line indicates the gene expression level corresponding to the highest cells’ frequency at the baseline (P0). d The box plot shows the median values, represented by a line across the boxes and vertical bars of min-max range distribution, of TCM-associated gene score (IL7R, TCF7, GATA3, CXCR3, EOMES) identified in the expanded γδ TCR clones and calculated for all γδ T cell clusters (c0-c9). For the statistical analysis, the unpaired parametric ANOVA test was used, and only significant P values (*) are shown: ****P < 0.0001. e Functional analysis of Vδ2 T cell response upon SARS-CoV-2 Prot_S peptides stimulation in vitro. The box plots show the fold change in IFNγ expression (MFI) analyzed by flow cytometry, comparing the first and second stimulations with PepTivator® SARS-CoV-2 Prot_S peptides in αβ T cells (left panel) and Vδ2 T cells (right panel) (n = 6). The graphs show median values represented by a line across the boxes and vertical bars of the min-max range distribution. The fold changes were calculated as the ratio of the MFI value of the 1st stimulation over its unstimulated control and the 2nd stimulation over its un-restimulated control. For the statistical analysis, the paired nonparametric Student t-test was used. For all graphs, statistic values are represented as P values (*): *P ≤ 0.05 and ****P < 0.0001.

Overall, our data strongly support the development of TCM Vδ2 T cells as result of SARS-CoV-2 vaccination.

To provide additional evidence for the Vδ2 T cell memory response identified through scRNA-seq analyses, we performed in vitro experiments by measuring the response of circulating Vδ2 T cells isolated from healthy donors upon stimulations with SARS-CoV-2 Spike-peptides (Supplementary Fig. 7b, c). There was significant donor-to-donor variation in the Vδ2 T cell response, however, we observed that cells after the second exposure to Spike-peptides exhibited slightly increased production of IFNγ compared to the first stimulation. In concordance, we detected an improved reactivity of αβ T cells to Spike-peptides following repeated stimulation, resulting in higher IFNγ production (Fig. 8e and Supplementary Fig. 7d).