CD4+ T cells exhibit distinct transcriptional phenotypes in the lymph nodes and blood following mRNA vaccination in humans

Abstract

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection and mRNA vaccination induce robust CD4⁺ T cell responses. Using single-cell transcriptomics, here, we evaluated CD4⁺ T cells specific for the SARS-CoV-2 spike protein in the blood and draining lymph nodes (dLNs) of individuals 3 months and 6 months after vaccination with the BNT162b2 mRNA vaccine. We analyzed 1,277 spike-specific CD4⁺ T cells, including 238 defined using Trex, a deep learning-based reverse epitope mapping method to predict antigen specificity. Human dLN spike-specific CD4⁺ follicular helper T (T_FH) cells exhibited heterogeneous phenotypes, including germinal center CD4⁺ T_FH cells and CD4⁺IL-10⁺ T_FH cells. Analysis of an independent cohort of SARS-CoV-2-infected individuals 3 months and 6 months after infection found spike-specific CD4⁺ T cell profiles in blood that were distinct from those detected in blood 3 months and 6 months after BNT162b2 vaccination. Our findings provide an atlas of human spike-specific CD4⁺ T cell transcriptional phenotypes in the dLNs and blood following SARS-CoV-2 vaccination or infection.

Extended Data Figure 1.. Sample preparation diagram and representative flow cytometry plots of cell purity following magnetic cell enrichment.

a. Blood and dLN samples from BNT162b2 mRNA vaccinated cohort. b. Peripheral blood mononuclear cell enrichment strategy for BNT162b2 mRNA vaccinated or SARS-CoV-2 infected donors.

Extended Data Figure 2.. Reference-based T cell annotations for UMAP in Figure 1b.

a. Density plots showing the relative distribution of ProjecTIL-based CD4⁺ T cell labels. b. Density plots showing the relative distribution of ProjecTIL-based CD8⁺ T cell labels. Individual gray dots indicate individual cells matching the label.

Extended Data Figure 3.. Performance metrics for Trex autoencoder models by approach and chain.

For the given hyperparameter, models were trained on 2e5 random sequences with 10 epochs for minimal Kullback-Leibler divergence value. a. Mean square error of models after training varying the latent dimensions (left panel) and batch size (right panel) with different learning rates. b. Kullback-Leibler divergence values of models after training varying the latent dimensions (left panel) and batch size (right panel) with different learning rates. c. Evaluations of fidelity of models to return unique values using novel sequences for TRA and TRB chains across all models in Trex. Novel sequences were randomly sampled and bootstrapped a total of 10 times. d. Distribution of computational time for model application across the models, chains, and bootstraps.

Extended Data Figure 4.. Comparison of Trex co-embedding approach with clonotype neighbor graph analysis (CoNGA).

a. Schematic representation of the CoNGA pipeline that generates nearest neighbors of clones using both edit-distance-based TCR networks and gene expression (GEX) networks. b. Resulting UMAPs for CoNGA-based dimensional reduction using gene expression or edit-distance-based TCR with denoted locations of previously identified spike-specific clones. c. Nearest-neighbor overlap using the Dice (left) and Jaccard (right) index of the 10 nearest neighbors defined by CoNGA and by the co-embedding with Trex. d. Breakdown and distribution of TCR-based clusters using CoNGA TCR output or Trex latent dimensions. Blue colored data indicate the relative proportion of clusters with spike-specific clones with a summary of the graphed values to the right of each bar chart. e. Trex-based latent dimensional clusters with proportion filled by the respective CoNGA TCR-based clusters. f. Distribution and relative size of the candidate TCRs and related sequences (edit distance ≤ 2) selected in Figure 3 for both the CoNGA-based TCR clusters (upper panel) and Trex-based clusters (lower panel).

Extended Data Figure 5.. Confirmation of TCR candidates’ specificity for SARS-CoV-2 spike.

Each TCR candidate’s variable gene regions were cloned with murine T cell receptor (mTCR) constant regions into a retroviral transduction vector and resultant retroviruses were used to transduce primary human CD4⁺ T cells. Positive results from intracellular cytokine stain mapping of the spike protein with overlapping peptides are shown. Gating was first performed on total live single cells, then on CD3⁺CD4⁺ T cells, and finally on mTCR beta chain (mTCRb) positive candidate TCR-transduced cells. Unstimulated background cytokine expression, positive control phorbol 12-myristate 13-acetate (PMA) and Ionomycin cytokine expression, and top cytokine expression to individual 17-mer peptides used for total spike proteome mapping are shown for each TCR candidate (a-e). Representative surface stain of unstimulated TCR2-transduced CD4⁺ T cells with the S_167-180 DPB1*04:01 HLA-class II tetramer is shown (right panel in b). Each experiment shown is representative of two independent TCR transduction and mapping experiments.

Extended Data Figure 6.. Confirmation of TCR candidate HLA restriction.

NFAT-GFP reporter Jurkat T cells transduced with candidate TCR expressing retrovirus were sort purified and maintained as clonal cell lines. a. Reporter Jurkat lines or b. Transduced primary human CD4⁺ T cells were co-cultured with spike peptides identified in Extended Data Figure 5 presented in the context of various K562-based aAPC cell lines expressing single HLA class II alleles. Cells were gated on total live single cells, then on CD3⁺ cells. In b. the top panels show the frequency of retrovirally transduced (murine TCR beta constant region expressing, mTCR⁺) primary human CD4⁺ T cells that were gated on prior to evaluation of intracellular cytokine staining in the bottom panels. Red asterisks denote positive responses for each TCR line.

Extended Data Figure 7.. Circulating blood spike-specific CD4⁺ T cells induced early after primary SARS-CoV-2 infection were similar regardless of illness severity.

Comparison of circulating blood spike-specific CD4⁺ T cells during acute (day 18 to 36 post-onset of disease symptoms) infection between donors with moderate (350-041, 350-117 and 350-400, n=3) versus severe (350-065, 350-084 and 350-397, n=3) infection. Statistical significance was based on bootstrapping 1,000 times to form a null distribution. * adjusted two-tailed permutation test p-value < 0.05.

Figure 1.. BNT162b2 mRNA vaccination induces spike-specific CD4⁺ T cells with diverse transcriptional phenotypes in the blood and dLN.

a. Schematic showing sample collection time points for dLN FNA and peripheral blood collection from six donors 368-01a, 368-04, 368-13, 368-16, 368-20 and 368-22 at day 28, day 60, day 110 and day 201 post-first dose of the BNT162b2 vaccine and graph showing the number of cells isolated from each donor from dLN or blood. For each sample collection, a technical replicate was performed and sequenced. b. UMAP of 219,283 dLN and blood CD3⁺ T cells with paired TRA-TRB sequences that passed quality control filtering. Cluster annotation based on canonical subtype markers and automated annotation using SingleR and ProjecTIL. c. Gene-weighted density estimates overlaid on the UMAP coordinates for the T cell markers CD4, CD8A, CCR7, SELL, CXCR5, ICOS, FOXP3, IL2RA, CTLA4 and PDCD1. d,e. Relative cellular density at day 28, 60 110 and 201 (d) and in blood and dLN (e) in CD3⁺ T cells as in b. f. Localization of spike-specific TCRs identified in ref^, along the UMAP. S_167-180-specific TCRs are highlighted in blue; other spike-specific TCRs are in red. g. Alignment of TRA and TRB CDR3 motifs for S_167-180-specific and other spike-specific TCRs.

Figure 2.. Diverse CD4⁺ T_FH cell and CD4⁺ T_FM cell transcriptional phenotypes are detected in the dLN after BNT162b2 mRNA vaccination.

a. UMAP of the subset of CD4⁺ T_FH cells and CD4⁺ T_FM cells from Fig. 1b found in the dLN of 368-01a, 368-04, 368-13, 368-16, 368-20 and 368-22 on day 28, 60, 110 and 201 post-first dose of the BNT162b2 vaccine. b. Gene-weighted density estimates of the indicated transcripts overlaid on the UMAP as in a. c. Top 8 or fewer differentially expressed genes in clusters c0-c11 based on UMAP as in a. Dot size represents the percentage of cells expressing the gene, and color is assigned based on scaled expression value. d. Heatmap of median gene set enrichment for the significantly altered gene sets in clusters c0-c4 and c6-c11 as in a., K1-K4 represent k-means clustering of gene sets with general summaries of groupings listed to the right of each group. Significance defined as adjusted p-value < 0.05 via two-way ANOVA. e. Normalized gene set enrichment values in clusters c0-c4 and c6-c11 at day 28, 60, 110 and 201 post-first dose of the BNT162b2 vaccine. Colors indicate individual gene sets. f. Circos plot showing overlap of unique individual TCR clonotypes at day 28, 60, 110 and 201 post-first dose of the BNT162b2 vaccine, with ribbons between clusters representing overlapping clonotypes. g. Relative CD4⁺ T_FH and CD4⁺ T_FM cells density at day 28, 60, 110 and 201 post-first dose of the BNT162b2 vaccine. S_167-180-specific TCR represented by white dots. Shown is the percentage of S_167-180-specific CD4⁺ T_FH and CD4⁺ T_FM cells among total CD4⁺ T_FH and CD4⁺ T_FM cells at each time point.

Figure 3.. Coembedding of single-cell TCR and RNA values from dLN CD4⁺ T_FH and CD4⁺ T_FM using Trex identifies spike-specific responses.

a. Graphical representation showing the computational embedding of TCR CDR3 amino acid sequences and single-cell RNA from CD4⁺ T_FH cells to generate a heat-diffusion-based manifold of CD4⁺ T_FH cells. Matrices are rescaled based on nearest neighbors and the corrected values are then used for dimensional reduction. b. PHATE projection of the tri-modal (RNA, TRA and TRB) embedding of dLN CD4⁺ T_FH cells by clonotype. Total number of each unique clonotype represented by dot size. c. Representative RNA expression of CD4⁺ T_FH cell marker genes (BCL6, CXCL13, CXCR5, ICOS, MAF, PDCD1, EMP3 and KLF2) overlaid onto the PHATE projection as in b. d. The location of spike-specific^, TRA (upper) and TRB (lower) CDR3 in the PHATE projection. e. Alignment of TRA and TRB motifs in spike-specific TCRs^, from PHATE-defined clusters Trex-C0, Trex-C1 and Trex-C3. f. The location within the PHATE projection (left) and sequence (right) of five candidate spike-specific TCR clonotypes derived from PHATE-defined clusters Trex-C0 and Trex-C1 that have one TCR chain appearing in > 1 donor and have not been previously described as specific for SARS-CoV-2 spike.

Figure 4.. Spike-specific CD4⁺ T_FH cell transcriptional phenotypes detected in the dLN change over time.

a. Median gene set enrichment heatmap showing immune-related gene sets in dLN CD4⁺ T_FH cells at day 28, 60 and 201 post-first dose of the BNT162b2 vaccine with 5 distinct clusters defined by k-means (K1-K5). b. Volcano plot of differential gene expression in dLN spike-specific CD4⁺ T_FH cells at day 28 (n=94) and day 201 (n=70) post-first dose of the BNT162b2 vaccine. Size is based on the change in the percentage of cells expressing the gene at day 28 compared to day 201. Statistical testing performed using two-sided MAST testing without adjustment for multiple comparisons. c. Clonal proportion of the spike-specific dLN CD4⁺ T_FH cell repertoire at day 28, 60 and 201 post-first dose of the BNT162b2 vaccine. Top number on the bar plot indicates the number of clones at the specific time point shared across more than one time point, whereas the number on the bottom indicates the number of spike-specific clones unique to the time point.

Figure 5.. TCR sequencing reveal limited overlap in clonal TCR repertoire between dLN and blood 3-6 months after mRNA vaccination.

a. Representation of total T cell clonal overlap between blood and dLN in donors 368-01a, 368-13 and 368-22 (blood samples from all three donors on days 110 and 201, dLN samples from 368-01a on days 28, 60, 110 and 201, dLN samples from 368-13 on days 60 and 110, and dLN samples from 368-22 on days 60 and 110 post-first dose of the BNT162b2 vaccine). Numbers indicate unique TCR clones and J is the calculated Jaccard stability index. b. Rarefication and extrapolation of all TCR clones included in a, for donors 368-01a, 368-13 and 368-22. The dotted line indicates the point of extrapolation, and the ribbon is the 95% confidence interval. c. Scatter plot showing the dLN or blood location and proportion of total TCR repertoire for each TCR clonotype found in donors 368-01a and 368-22. Overlapping clonotypes are indicated in yellow.

Figure 6.. Total spike-specific blood CD4⁺ T cells are transcriptionally distinct from the total spike-specific CD4⁺ T cell population in the dLN.

a. Volcano plot of differential gene expression between all spike-specific CD4⁺ T cells in dLN (n=533) and blood (n=938). Size of points is based on the difference in the percentage of cells expressing each gene in dLN compared to blood. Statistical testing was performed using two-sided MAST testing with adjustment for multiple comparisons. b. Z-scaled median gene set enrichment heatmap for immune-related gene sets found in dLN and blood spike-specific CD4⁺ T cells at day 110 and day 201 post-first dose of the BNT162b2 vaccine.

Figure 7.. Circulating blood spike-specific CD4⁺ T cells induced by infection are transcriptionally distinct from those induced by mRNA vaccination.

a. Schematic showing the timeline for blood sample collection 3-months and 6-months post-first dose of the BNT162b2 vaccine in donors (368-01a, 368-04, 368-13 and 368-22, n=4) or following SARS-CoV-2 infection in individuals who developed moderate (350-041, 350-117 and 350-400, n=3) or severe (350-065, 350-084 and 350-397, n=3) disease. b. UMAP projection of all circulating blood spike-specific CD4⁺ T cells in all donors as in a., including 693 from SARS-CoV-2-infected and 804 cells from BNT162b2-vaccinated donors (1,497 total cells). c. UMAP projection (left) and proportion breakdown (right) of cluster composition for circulating blood spike-specific CD4⁺ T cells in donors 368-01a, 368-04, 368-13, 368-22, 350-041, 350-117 and 350-400 at 3 months and 6 months in SARS-CoV-2-infected (infected, n=289) and BNT162b2-vaccinated (vaccinated, n= 804) donors. Statistical significance was based on bootstrapping 1,000 times to form a null distribution. * corrected p-value < 0.05; ** corrected p-value < 0.01. d. Top 8 or fewer cluster-defining differentially expressed genes in clusters 0-8 as in b. e. TCR cluster assignments based on normalized Levenshtein distance of the CDR3 sequence across donors. Only cluster assignments with more than two clonotypes were retained. f. UMAP and proportion breakdown of circulating blood spike-specific CD4⁺ T cells in infected donors 350-041, 350-117 and 350-400 at day 18 to 36 (Early) and at 3-6 months (Late).