1. Computational and Systems Biology
  2. Immunology and Inflammation

Primary and secondary anti-viral response captured by the dynamics and phenotype of individual T cell clones

  1. Anastasia A Minervina  Is a corresponding author
  2. Mikhail V Pogorelyy  Is a corresponding author
  3. Ekaterina A Komech
  4. Vadim K Karnaukhov
  5. Petra Bacher
  6. Elisa Rosati
  7. Andre Franke
  8. Dmitriy M Chudakov
  9. Ilgar Z Mamedov
  10. Yuri B Lebedev
  11. Thierry Mora
  12. Aleksandra M Walczak
  1. Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Russian Federation
  2. Center for Precision Genome Editing and Genetic Technologies for Biomedicine, Pirogov Russian National Research Medical University, Russian Federation
  3. Center of Life Sciences, Skoltech, Russian Federation
  4. Institute of Immunology, Kiel University, Germany
  5. Institute of Clinical Molecular Biology, Kiel University, Germany
  6. Masaryk University, Central European Institute of Technology, Czech Republic
  7. V.I. Kulakov National Medical Research Center for Obstetrics, Gynecology and Perinatology, Russian Federation
  8. Moscow State University, Russian Federation
  9. Laboratoire de physique de l’École normale supérieure, ENS, PSL, Sorbonne Université, Université de Paris, and CNRS, France
https://doi.org/10.7554/eLife.53704
Share this article
Cite this article
  1. Anastasia A Minervina
  2. Mikhail V Pogorelyy
  3. Ekaterina A Komech
  4. Vadim K Karnaukhov
  5. Petra Bacher
  6. Elisa Rosati
  7. Andre Franke
  8. Dmitriy M Chudakov
  9. Ilgar Z Mamedov
  10. Yuri B Lebedev
  11. Thierry Mora
  12. Aleksandra M Walczak
(2020)
Primary and secondary anti-viral response captured by the dynamics and phenotype of individual T cell clones
eLife 9:e53704.
https://doi.org/10.7554/eLife.53704
  2. Download BibTeX
  3. Download .RIS

Abstract

The diverse repertoire of T-cell receptors (TCR) plays a key role in the adaptive immune response to infections. Using TCR alpha and beta repertoire sequencing for T-cell subsets, as well as single-cell RNAseq and TCRseq, we track the concentrations and phenotypes of individual T-cell clones in response to primary and secondary yellow fever immunization — the model for acute infection in humans — showing their large diversity. We confirm the secondary response is an order of magnitude weaker, albeit ∼10 days faster than the primary one. Estimating the fraction of the T-cell response directed against the single immunodominant epitope, we identify the sequence features of TCRs that define the high precursor frequency of the two major TCR motifs specific for this particular epitope. We also show the consistency of clonal expansion dynamics between bulk alpha and beta repertoires, using a new methodology to reconstruct alpha-beta pairings from clonal trajectories.

Introduction

T-cells play a crucial role in the immune response to pathogens by mediating antibody formation and clearance of infected cells, and by defining an overall response strategy. The specificity of T-cells is determined by the T-cell receptor (TCR), a heterodimer of alpha and beta protein chains. Genes for alpha and beta chains assemble in a random process of somatic V(D)J-recombination, which leads to a huge variety of possible TCRs (Murugan et al., 2012). The resulting diverse naïve repertoire contains T-cell clones that recognize epitopes of yet unseen pathogens, and can participate in the immune response to infection or vaccination. One of the best established models of acute viral infection in humans is yellow fever (YF) vaccination. Yellow fever vaccine is a live attenuated virus with a peak of viremia happening around day 7 after vaccine administration (Miller et al., 2008; Akondy et al., 2009). The dynamics of primary T-cell response was investigated by various techniques: cell activation marker staining (Miller et al., 2008; Blom et al., 2013; Kohler et al., 2012; Kongsgaard et al., 2017), MHC multimer staining (Akondy et al., 2009; Blom et al., 2013; James et al., 2013; Kongsgaard et al., 2017), high-throughput sequencing (DeWitt et al., 2015; Pogorelyy et al., 2018) and deuterium cell labelling (Akondy et al., 2017). Primary T-cell response sharply peaks around 2 weeks after YFV17D (vaccine strain of yellow fever virus) vaccination (Miller et al., 2008; Akondy et al., 2009; Kohler et al., 2012; Pogorelyy et al., 2018; James et al., 2013). The immune response is very diverse and targets multiple epitopes inside the YF virus (de Melo et al., 2013; Co et al., 2002; Akondy et al., 2009; James et al., 2013; Blom et al., 2013). An essential feature of effective vaccination is the formation of immune memory. Although most of the effector cells die shortly after viral clearance, YF-specific T-cells could be found in the blood of vaccinated individuals years (Akondy et al., 2009; Kongsgaard et al., 2017; James et al., 2013) and even decades after vaccination (Fuertes Marraco et al., 2015; Wieten et al., 2016). While the immune response to the primary vaccination has been much studied, there is only limited data on the response to the booster vaccination with YFV17D. Both T-cell activation marker staining and multimer staining show that the secondary response is much weaker than the primary one (Kongsgaard et al., 2017), but their precise dynamics, diversity, and clonal structure are still unknown.

In summary, previous studies provide insight into the macroscopic features of the T-cell response, such as total frequency of T-cells with an activated phenotype, or T-cells specific to a particular viral epitope on different timepoints after vaccination. However, with recently developed methods it is now possible to uncover the microscopic structure of the primary and secondary immune response, such as the dynamics and phenotypes of distinct T-cell clones, as well as the receptor features that determine the recognition of epitopes.

TCR repertoire sequencing allows for longitudinally tracking individual clones of responding T-cells irrespective of their epitope specificity. Single-cell RNAseq (scRNAseq) enables simultaneous quantification of thousands of transcripts per cell for thousands of cells, providing an unbiased characterization of immune cell phenotype. Single-cell TCR sequencing produces paired αβ repertoire data, and thus could help discover conserved sequence motifs in one or both TCR chains. These motifs encode TCR structural features essential to antigen recognition (Dash et al., 2017; Glanville et al., 2017). Information about complete TCR sequences allows homological modeling of TCR structure (Schritt et al., 2019), which can be used for binding prediction with protein-protein docking (Pierce and Weng, 2013). We combine longitudinal TCR alpha and beta repertoire sequencing, scRNAseq, scTCRseq, TCR structure modelling and TCR-pMHC docking simulations to get a comprehensive picture of primary and secondary T-cell response to the yellow fever vaccine – the in vivo model of acute viral infection in humans.

Results

Secondary T-cell response to the YFV17D vaccine is weaker but faster than the primary response

We sequenced TCR alpha and TCR beta repertoires of bulk peripheral blood mononuclear cells (PBMCs) and different T-cell subsets at multiple timepoints before and after primary and booster vaccination against yellow fever of donor M1 (Figure 1A). Clonotypes responding to the primary YF immunization were identified using the edgeR software as previously described (Pogorelyy et al., 2018). Briefly, the biological replicates of bulk PBMCs were used to estimate the noise in the TCR mRNA counts. Clonotypes were assumed YF-responding if they increased in concentration more than 32-fold (p<0.01, see Materials and methods) between any two timepoints before the peak of the primary response (days 0, 5, 10 and 15).

Figure 1 with 2 supplements see all
Download asset Open asset
Primary and secondary response to yellow fever vaccination.

(A) Experiment design. Blood was taken at multiple timepoints before and after primary and secondary immunization against yellow fever virus. Two biological replicates of PBMCs and different cell subpopulations (indicated below each day of blood draw) were isolated at all timepoints. cDNA TCR alpha and TCR beta libraries were sequenced on Illumina platform. (B) The number of significantly expanded TCR alpha and TCR beta clonotypes for both donors in comparison to day 0. For donor P30 the number of significantly expanded clones is lower, than observed in primary vaccinations (see Figure 1—figure supplement 2). (C) The fraction of YF-responding cells as a proportion of all T-cells, measured by cumulative frequency of YF-responding TCR alpha and beta clonotypes of donor M1 after first (light blue and dark blue) and second immunization (dashed light blue and dark blue), and donor P30 (orange and yellow), which had a second immunization 30 years after the first. (D) The fraction of CD4+ and CD8+ YF-responding cells, as a proportion of all T-cells of donor M1 during the primary and secondary response to YFV17D. No novel major expansions were observed after secondary immunization, see Figure 1—figure supplement 1.

Figure 1—source data 1

List of all sequencing libraries with summary statistics.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig1-data1-v2.xlsx
Figure 1—source data 2

Number of significantly expanded TCR alpha and TCR beta clonotypes between in comparison to day 0.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig1-data2-v2.csv
Figure 1—source data 3

YF-responding TCR alpha and TCR beta clonotypes of donors M1 and P30 identified by edgeR.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig1-data3-v2.xlsx

Overall we found 1580 TCR beta and 1566 TCR alpha clonotypes significantly expanded after the primary immunization, respectively occupying 6.7% and 7.8% of the sampled TCR repertoire of bulk PBMCs in cumulative frequency at the peak of the response (Figure 1B,C). As expected, both the numbers of responding clones and their cumulative frequencies were very similar for expanded clonotypes identified in bulk TCR alpha and beta repertoires. For simplicity in the following sections we focus on TCR beta repertoires, unless stated otherwise. In accordance with previous studies (Miller et al., 2008; Blom et al., 2013; Akondy et al., 2009; Kongsgaard et al., 2017; Pogorelyy et al., 2018), we show that during the primary response T-cells expanded intensely (with cumulative increase of about 950-fold) within 2–3 weeks after YF immunization. They subsequently contracted, but still exceeded baseline frequency 18 months afterwards.

We then tracked these YF-responding clonotypes identified during primary immunization before and after the second vaccination 18 months after the first one. The cumulative frequency of these clonotypes increased ≈2.5-fold at the peak of the response after the second immunization, reaching 0.5% of the TCR repertoire (Figure 1D, blue curve). The secondary response was weaker, but happened much faster than the primary one, with a peak frequency of responding clonotypes occurring on day 5 instead of day 15 after vaccination. To check if there was also recruitment of new clonotypes in the secondary response, we applied edgeR to timepoints from the second immunization only. Although we identified 73 additional responding clonotypes, their impact on the magnitude of the secondary response was negligible and we did not use them for further analyses (see Figure 1—figure supplement 1). Backtracking of these novel clonotypes showed that they also slightly expanded during the primary response but not enough enough to pass our significance and magnitude thresholds. In summary, we found no evidence of substantial recruitment of naive clones in the response to the booster vaccination.

Using sequenced CD4+ and CD8+ T-cell subsets, we attributed a CD4 or CD8 phenotype to each responding clone (see Materials and methods) and thus could track these two subsets separately. After booster immunization in donor M1, YF-responding CD4+ cells peaked earlier (day 5 vs day 10) and expanded much more (≈8 times vs. ≈1.5 times) than CD8+ T-cells (Figure 1D, green and pink curves). During primary immunization, the difference in response dynamics between CD4+ and CD8+ subsets was less prominent, as they both peaked on day 15. However, by day 21 CD4+ responding clones contracted much more (to 43.6% of peak frequency) than CD8+ clonotypes (87% of peak frequency). These observations confirm previous reports that the CD4 response precedes the CD8 response (Blom et al., 2013).

Secondary response to booster vaccination after 18 months and after 30 years have similar features

To see how long-lived T-cell memory response to YF can be, we recruited an additional donor (P30), who received the first YF-vaccine 30 years earlier and has not been in YF endemic areas for at least 28 years. From this donor, we collected bulk PBMCs and several T-cell subsets before and after booster immunization. Both the numbers of responding clonotypes (204 for TCR beta and 201 for TCR alpha) and the maximum frequency at the peak of the response (0.69%) were much lower than for any primary vaccinee both from this and other studies (Figure 1—figure supplement 2). Most of these clonotypes were low frequency or undetected before the second immunization, although a few were sampled in the memory repertoire prior to vaccination.

The response to the booster vaccination was characterized by a large expansion between days 0 and 5, and a peak on day 10, for both CD4+ and CD8+ T-cells. Overall the dynamics and the magnitude of this response was very similar to the response to the booster vaccination after 18 months we observed in donor M1 (Figure 1C), suggesting that protection against the virus was maintained even after 30 years.

Diversity of clonal time traces in primary and secondary responses

Our approach allows us to estimate the contribution of individual clones to the total response. We already showed that the overall response strength to secondary immunization was an order of magnitude lower compared to the primary response. However, several clones showed remarkable expansion rates and peak frequencies, comparable to the ones observed in primary immunization. Such clones were observed in both donors upon secondary immunization after 18 months and 30 years (Figure 2A and B, Figure 2—figure supplement 1). We traced each single clone during primary and secondary response in donor M1. The concentration of clonotypes prior to the booster immunization correlated well (Pearson r = 0.46 p<0.0001) with their concentration on day 45 after primary immunization (Figure 2—figure supplement 2) suggesting a uniform contraction rate for all clones resulting in a half-life of 158 ± 12.7 days for the YF-specific T-cell subpopulation. Previously, Akondy et al. using deuterium labeling of cells specific to the immunodominant epitope NS4B214-222 (as determined by a A02-NS4B214-222-multimer binding assay) showed a very similar half-life of 123 days (Akondy et al., 2017).

Figure 2 with 4 supplements see all
Download asset Open asset
Diversity of individual clonal trajectories in primary and secondary responses.

(A, B) Frequency of each YF-responding clonotype in bulk TCR repertoire as a function of time. Individual clones show remarkable expansion after the primary response (A, left panel) and secondary response both 18 months (A, right panel) and 30 years (B) after the primary vaccination. The ten most abundant (by peak frequency) CD4+ and CD8+ YF-responding clonotypes are shown for each vaccination. Clonal traces for all YF-responding clonotypes are shown in Figure 2—figure supplement 1. Color indicates the time of the response peak for each clonotype: blue for a peak at day 5, pink at day 10, green at day 15 and purple at day 21. Despite overall heterogeneity in clonal traces, more clones peak at early timepoints during the secondary response. Heterogeneity in clonal traces allows for expanded clones identification and computational alpha-beta TCR pairing (Figure 2—figure supplement 4).

Figure 2—source data 1

Concentrations of YF-responding clonotypes for donor M1 on all timepoints.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig2-data1-v2.xlsx
Figure 2—source data 2

Concentrations of YF-responding clonotypes for donor P30 on all timepoints.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig2-data2-v2.csv

It was previously reported that only 5–6% of YF-responding clones are preserved as immune memory, with the preferential recruitment of large clones (DeWitt et al., 2015). By contrast, in our sample we could re-identify 96% of CD4+ and 88% of CD8+ clones that responded to the primary immunization in at least one sample after the booster immunization. This suggests that practically all the diversity of the responding repertoire is maintained in memory. The larger fraction of re-identified YF-responding clones in comparison to previous work may be explained by the sampling depth. Sequencing more T-cells will lead to the re-identification of even more YF-responding clonotypes.

We then wanted to characterize how these persistent clonotypes responded to the booster vaccination. Interestingly, we found that the largest YF-specific CD8+ clones did not expand in response to the booster vaccine. Instead, the most expanded clonotypes were rare prior to the booster immunization (Figure 2—figure supplement 3A). The situation was different for CD4+ cells: both high and low-frequency CD4+ clones expanded in response to the booster immunization (Figure 2—figure supplement 3B).

The specific features of clonal trajectories shared by YF-responding clones make it possible to distinguish them from non-expanding clones, using unsupervised clustering (see Figure 2—figure supplement 4AB and Materials and methods). This method shows good concordance with edgeR and works also without biological replicates. In addition, we demonstrated that the heterogeneity of clonal trajectories could be leveraged to computationally pair alpha and beta chains from from bulk alpha and beta sequencing data, by exploiting the similarity of trajectories of alpha and beta clonotypes belonging to the same clone (see Figure 2—figure supplement 4C and Materials and methods).

TCR sequencing shows the transition of clonotypes between memory subpopulations

Several studies have reported subsets of long-lived memory YF-specific T-cells, whose concentration remained stable for years (Fuertes Marraco et al., 2015; Akondy et al., 2017). It was shown that these long-lived memory cells are the progenies of effector cells, which divide vigorously during the peak of the response to the vaccine (Akondy et al., 2017). TCR sequences can be used as ‘barcodes’ to measure transitions between different memory subsets after YF immunization, defined by their surface markers revealed by flow cytometry.

We isolated with FACS (see Figure 3—figure supplement 1 for the gating strategy) and sequenced TCR repertoires of 3 conventional T-cell memory subpopulations (Fuertes Marraco et al., 2015; Appay et al., 2008): effector memory (EM, CCR7-CD45RA-), effector memory re-expressing CD45RA (EMRA, CCR7-CD45RA+), and central memory (CM, CCR7+CD45RA-) on days 0, 15, 45, and 18 months after the primary vaccination of donor M1 and on days 0, 15, and 45 after the second vaccination of donor P30. On day 45 we also isolated and sequenced the repertoire of the recently described Tscm (T-cell stem cell-like memory) subset (CCR7+CD45RA+CD95+).

On day 0, the concentration of almost all YF-responding clonotypes was too low to be detected in any of these subpopulations. However, we were able to calculate the distribution of YF-responding clonotypes between these phenotypes after immunization. In agreement with previous studies the memory status of T-cell clones was tightly correlated with their CD4/CD8 status (Sathaliyawala et al., 2013; Thome et al., 2014). The vast majority of CD4+ T-cells were distributed between EM and CM, with <1% in EMRA, while CD8+ T-cell clones were predominantly found in EM and EMRA with ∼2% in CM. This difference also held for YF-responding clones (Figure 3A). While for most CD8+ clonotypes in the total repertoire EM/EMRA phenotypes were stable between day 15 and day 45 (Figure 3B, and Figure 3—figure supplement 2A,C), the distribution of CD8+ YF-responding clones between memory subsets was significantly shifted towards the EMRA phenotype (Figure 3C). This shift results from two processes: the rapid decay of EM cells (Figure 3—figure supplement 2B) and the phenotype switch from EM to EMRA (Figure 3—figure supplement 2D). Almost all YF-responding CD8+ clones detected 18 months after the first immunization corresponded to the EMRA phenotype (among 71 clones found in more than three copies in bulk repertoire at day 0 before second vaccination, 41 were found only in the EMRA subset, four only in EM, and six in both). For CD4+ T-cells, we did not observe any trend in phenotype switching between days 15 and 45 after the vaccination. We hypothesize that switching from EM to CM phenotype was masked due to homing of CM cells to lymphoid organs, defined by the expression of the CCR7 chemokine receptor.

Figure 3 with 2 supplements see all
Download asset Open asset
Distribution of clonotypes in memory subsets.

(A) Each color bar shows the estimated distribution of T-cell clones between memory subpopulations for a set of CD4+ (left panel) and CD8+ (right panel) clonotypes for donors M1 (top) and P30 (bottom) on day 15 and day 45. Each panel shows the 10 most abundant YF-responding clones in each donor on day 45, which are present in at least one memory subpopulation on both day 15 and day 45. (B) Estimated concentration of CD8+ clones with a given phenotype at different timepoints in the bulk PBMC repertoire, for non-YF-responding clonotypes and (C) YF-responding CD8+ clonotypes (Mann Whitney U-test, EM: p-value = 2.1 · 10-12, EMRA: p-value = 1.2 · 10-6). Only clones with 30 or more Unique Molecular Identifiers (see Materials and methods) in bulk repertoires on day 45 were used for the analysis.

Figure 3—source data 1

Distribution of 10 most abundant CD4+ and CD8+ YF-responding clonotypes from donors M1 and P30 between memory subsets.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig3-data1-v2.xlsx
Figure 3—source data 2

Concentrations of non-YF-responding CD8+ clones in EM and EMRA subsets on day 15 and day 45.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig3-data2-v2.csv
Figure 3—source data 3

Concentrations of YF-responding CD8+ clones in EM and EMRA subsets on day 15 and day 45.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig3-data3-v2.csv

The response to a single immunodominant epitope can contribute to up to 60% of the total response

It was previously shown that in HLA-A02 donors the NS4B214-222 LLWNGPMAV immunodominant epitope elicits the strongest CD8+ T-cell response (Akondy et al., 2009; Wieten et al., 2016; Kongsgaard et al., 2017; Blom et al., 2013). Using an A02-pMHC-dextramer, we isolated NS4B-specific CD8+ T-cells from both donors (Figure 4—figure supplement 1A,B) and applied TCR sequencing to get their unpaired TCR alpha and TCR beta repertoires. We obtained ≈2100 alpha and ≈2000 beta functional receptor chains, one of the largest datasets for TCRs with a single specificity. YF-responding clonotypes identified by edgeR as expanded between timepoints are not restricted to any particular YF epitope and represent the repertoire targeted towards many different peptides presented by different HLA alleles. This allows us to quantify the relative contribution of NS4B-specific T-cells to the total anti-YF response. At the peak of the response, approximately 24% of all YF-responding CD8+ T-cells were specific to NS4B in the donor vaccinated 30 years ago (P30), and up to 60% in the first time vaccinee (M1) (Figure 4A). However, NS4B-specific clonotypes could not be distinguished from other YF-responding clonotypes from their time traces alone, as they both responded with similar dynamics (Figure 4—figure supplement 2).

Figure 4 with 4 supplements see all
Download asset Open asset
Response to the immunodominant yellow fever epitope NS4B214-222.

(A) Fraction of all T-cells corresponding to CD8+ YF-responding TCRβ clonotypes (solid lines) and CD8+NS4B-specific clonotypes (dashed lines) as a function of time post-vaccination (x-axis). Sequence similarity networks for TCR alpha (B) and beta (C) of NS4B-positive cells. Each vertex is a TCR amino acid sequence, connected with an edge if they differ by fewer than two mismatches. The size of the vertex indicates its degree. Vertices of zero degree are not shown. Color and text boxes indicate V-segments that are significantly enriched (exact Fisher test, Benjamini Hochberg adjusted p<0.001) in usage in epitope-specific cells compared to the bulk repertoire. NS4B-specific TCR alpha (D) and TCR beta (E) chains (red histograms) have biases in CDR3 length in comparison to bulk TCR repertoire of CD8+ cells (overlayed blue histograms). (F) Network of single-cell paired TCR alpha (blue) and TCR beta (red) of NS4B-specific TCRs. Vertices of the same color are connected if there are less than two mismatches in TCR chain amino acid sequence. An edge between vertices of different color represents the pairing of alpha and beta. The biggest alpha cluster (blue in the center) corresponds to the TRAV12-2 cluster on B, and it pairs with many dissimilar beta chains. The biggest beta cluster (top left in red) corresponds to the TRBV9 cluster of C. (G) Pairing of V-segments of TCR alpha (left) to V-segments of TCR beta (right) in scTCRseq of NS4B-specific T-cells. The height of each box is proportional to the number of unique clones with this V-segment. The width of ribbons is proportional to the frequency of TRAV-TRBV combination. NS4B-specific TCRs have two main binding modes, defined by TRAV12 segment family paired to almost any TRBV (blue) and by TRAV27 segment paired preferentially with TRBV9 (pink).

Figure 4—source data 1

NS4B-specific TCR alpha and TCR beta clonotypes from donors M1 and P30.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig4-data1-v2.xlsx
Figure 4—source data 2

Paired NS4B-specific alpha/beta TCR clonotypes.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig4-data2-v2.csv

Sequence analysis and structural modeling of NS4B-specific TCRs reveals two motifs with distinct peptide binding modes

We next asked whether there are distinct features in the sequence of NS4B-specific TCRs, which might explain the immunodominance of this epitope. Figure 4B and C show sequence similarity networks for TCR alpha and TCR beta chains of NS4B-specific clonotypes. The TCR alpha repertoire shows biased V-usage and complementarity determining region 3 (CDR3) lengths (Figure 4D). TRAV12-2, TRAV12-1, TRAV27, and TRAV17 gene usage were significantly enriched in the NS4B-specific TCRs (exact Fisher test, Benjamini Hochberg adjusted p<0.001), with more than 45 percent of the clonotypes expressing TRAV12-2, in comparison to just 4.5% of TRAV12-2 in the total CD8+ TCR repertoire. Beta chains formed several distinct clusters of highly similar sequences, with significant but less marked V-usage biases towards TRBV9, TRBV15, and TRBV6-1/2, as well as some bias in the length distribution (Figure 4E). Almost 37% of NS4B-specific clonotypes used TRBJ2-7.

We next asked how these clusters of highly similar sequences in the alpha and beta NS4B-specific repertoires corresponded to each other. Prior to booster immunization, we isolated NS4B-specific T-cells from donor M1 (Figure 4—figure supplement 1C) and performed single-cell RNA sequencing (scRNAseq) and single-cell paired TCR sequencing (scTCRseq). We collected data from 3500 cells corresponding to 164 clonotypes (see Materials and methods). Figure 4F shows a joint similarity network for TCR alpha and TCR beta chains, with both intra-chain sequence similarity and inter-chain pairings. Alpha-beta pairing seemed to be mostly random, with some exceptions: for instance, specific TCRs using the most dominant TRAV12-2 alpha motif were paired with many different beta chains with a broad usage of V-segments (Figure 4G and Figure 4—figure supplement 3A), but with a restricted CDR3β length of 13–14 amino acids. TCRs using TRAV27 and TRBV9 segments were also preferentially paired with one another (Figure 4—figure supplement 3C). Clustering of paired sequences using the TCRdist measure (Figure 4—figure supplement 3B) resulted in two large clusters corresponding to these two major motifs with conserved V-usage.

The preferential usage of the TRAV12 family was reported before for TCRs responsive to the NS4B epitope (Bovay et al., 2018; Zhang et al., 2018). It was speculated (Bovay et al., 2018), that the CDR1α of this V-segment forms contacts with the peptide. To test this hypothesis, we modeled the 3D structures of clonotypes from scTCRseq using the Repertoire Builder server (Schritt et al., 2019) and then docked the resulting model structures using RosettaDock (Lyskov and Gray, 2008) to the HLA-A02 pMHC complex structure, recently solved using X-ray crystallography (Bovay et al., 2018), see Materials and methods for details. Models of TCR-pMHC complexes showed that the TRAV12-2 TCRs formed more contacts with the peptide using CDR1α loops, and fewer contacts with CDR3α loops, in comparison to TRAV27 TCRs (Figure 4—figure supplement 4A). Interestingly, CDR3α sequences of TRAV12-2 TCRs were very similar to the ones observed in the repertoire of the same donor prior to the immunization, suggesting absence of epitope-driven selection of the CDR3α of these TCRs (Figure 4—figure supplement 4B). Based on these results, we hypothesize that TCRs using TRAV12 and TRAV27 motifs represent two independent and distinct solutions to the binding of the NS4B epitope.

scRNAseq of NS4B-specific T-cells reveals two distinct cytotoxic phenotypes

Next we used the scRNAseq gene expression data to investigate the phenotype of specific T-cells in finer detail. While almost all NS4B-specific clonotypes 18 months after vaccination belonged to the conventional EMRA subset, scRNAseq revealed huge heterogeneity of gene expression inside this population. Unsupervised clustering by Seurat 3.0 software (Stuart et al., 2019; Butler et al., 2018) (see Materials and methods) revealed three sub-phenotypes of NS4B-specific cells (Figure 5A).

Figure 5 with 3 supplements see all
Download asset Open asset
Phenotypic diversity of NS4B-specific cells 18 months after yellow fever immunization.

(A) 2D t-SNE visualization of unsupervised clustering (Seurat analysis) of RNAseq data based on 2000 most variable genes shows three distinct clusters of NS4B-specific cells. (B) The heatmap of top 15 significantly enriched genes of single cells in clusters 1 and 2 defined by the MAST algorithm. The panel above the heatmap identifies the cluster identity of the cells. (C) Gaussian kernel density estimate for the relative fraction of cells belonging to cluster one for each clonotype. Blue distribution shows the theoretical prediction under the null hypothesis: clonotype labels were shuffled between cells (1000 permutations). The observed distribution is flatter than the theoretical one, indicating the presence of clonotypes with either a minority or a majority of cells belonging to cluster 1 (χ2 -test with MC-estimated p-value=0.0005).

Figure 5—source data 1

Differentially expressed genes between NS4B-specific cells 18 months after vaccination.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig5-data1-v2.csv
Figure 5—source data 2

Differentially expressed genes between NS4B-specific clonotypes 18 months after vaccination.

https://cdn.elifesciences.org/articles/53704/elife-53704-fig5-data2-v2.csv

Overall we found 166 genes that were differentially expressed according to the MAST algorithm (Finak et al., 2015) between these clusters (Figure 5B). Cells from cluster one showed high expression of cytotoxicity related genes GZMB, GNLY, GZMH, NKG7, PRF1, CX3CR1, SPON2, KLRD1, Hobit and T-bet transcription factors (Figure 5—figure supplement 1A). The combination of these genes also suggests that this cytotoxicity is mediated by the perforin pathway. The second cluster of cells is enriched in genes such as CCR7, TCF7, SELL, JUNB, LEF1, and especially IL7R which are essential for long-term survival and maintenance of memory T-cells (Figure 5—figure supplement 1B; Jeannet et al., 2010; Zhou et al., 2010; Kaech et al., 2003; Jung et al., 2016; Schluns et al., 2000). However, these cells also express unique markers related to cytotoxicity: GZMK, LTB as well as KLRG1, KLRB1, T-bet, and GZMH, albeit at lower levels than cells in cluster 1.

Very similar clusters of genes were found in single-cell RNAseq analysis of CD4-cytotoxic lymphocytes EMRA cells (Patil et al., 2018). The expression pattern of granzymes and killer-like receptors in our clusters suggests that cells in cluster two may be the precursors of cells in cluster 1. The expression of GZMK (enriched in cluster 2) was shown to be prevalent in early memory stages (Harari et al., 2009; Bratke et al., 2005), while high levels of GZMB, GZMH, KLRB1, KLRG1, and ADGRG1 (enriched in cluster 1) are associated with more terminally differentiated memory cells with higher cytotoxic potential (Truong et al., 2019; Takata and Takiguchi, 2006). Interestingly, cluster two has higher expression of genes encoding ribosomal proteins, which were recently reported to be a feature of memory precursor cells (Araki et al., 2017). The transition of cells between the two clusters is also supported by the existence of cluster 3, which shows intermediate gene expression of cluster 1 and 2 markers, and thus may represent cells gradually changing phenotype.

For each cell from the scRNAseq experiment, we obtained matched scTCRseq results. We wondered whether the TCR clonotype influenced cell gene expression profile. Interestingly, the distribution of clonotypes between clusters was not random (χ2 -test with MC-estimated p-value=0.0005): some clonotypes showed a clear preference for one of the phenotypes (Figure 5C). To match single-cell gene expression data with measurements of clonotype concentrations obtained with TCRseq, we averaged mRNA counts over the all cells of the same clonotype, and repeated the differential gene expression analysis (see Materials and methods). We obtained two clusters of clonotypes with the same enriched genes (Figure 5—figure supplement 2A) as observed for clusters of single cells (Figure 5B), confirming the association of phenotype and clonotype. Clonotypes from both clusters expanded following the second immunization, indicating that both phenotypes are capable of response. Clonotypes associated to cluster one had larger frequencies both on day 45 after the first vaccination (Figure 5—figure supplement 2B, left), and 18 months later before the booster shot (Figure 5—figure supplement 2B, right), than clonotypes associated to cluster 2. This result suggests that even for T-cells recognizing the same epitope, particular clones are linked to particular memory phenotype.

Discussion

In this study, we applied high-throughput sequencing of TCR alpha and TCR beta repertoires to track T-cell immune response to primary and secondary immunization with yellow fever vaccine. This approach does not require previous knowledge of TCR specificity and thus allows us to quantify and compare the response of individual T-cell clones recognizing different epitopes on the same scale.

We found that up to 60% of all responding CD8+ T-cells were specific to a single immunodominant peptide. Several studies reported high precursor frequency of T-cells reactive to this epitope (Zhang et al., 2018; Bovay et al., 2018). Bovay et al. (2018) recently suggested that recognition of antigenic peptide through the germline-encoded CDR1 loop of the TRAV12 segment is one of the main reasons for high precursor frequency. This hypothesis is supported by our TCR structural modeling and TCR-pMHC docking simulations, as well as by the analysis of the NS4B-specific T-cell repertoire. We also identified an additional motif defined by TRAV27+TRBV9+ TCRs. It will be interesting to investigate if these two motifs differ in binding affinity or are susceptible to potential escape mutations that can appear in the antigenic peptide. Another question is how diverse is the level of clonal response to the YF vaccine in HLA-A02 negative donors, and what fraction of the response is directed towards the most immunodominant epitopes in the context of other HLA types.

Most previous studies focused on TCR beta repertoires, partially because the diversity of TCR beta is higher, making it a better marker for clonal tracking (Chu et al., 2019). We found that TCR alpha may be used for clonal tracking as well, giving almost the same results as TCR beta in terms of the number of expanding clonotypes and their cumulative fractions on different timepoints. In the particular case of the response of HLA-A02 donors to the YF vaccine, the TCR alpha repertoire turned out to be even more informative, as T-cells responding to the immunodominant epitope predominantly use certain TRAV segments.

One of the major limitations of bulk TCR sequencing is that the resulting repertoires are unpaired, while TCR specificity in most cases is defined by the combination of alpha and beta chains. We show that the simultaneous sequencing of bulk alpha and beta repertoires performed on many timepoints allows us to make predictions on alpha-beta pairing. Even with the rise of single-cell sequencing, this method might still be of interest since most available single-cell platforms can only analyze limited numbers of 104-105 cells. In addition, these experiments are still expensive in comparison to the bulk TCR sequencing, which enables the profiling of millions of lymphocytes more cheaply.

We found that ≈ 90% clonotypes responding to primary immunization were present in peripheral blood 18 months after immunization. Recently, Akondy et al. (2017) showed using deuterium cell labeling that long-survived memory cells have a history of intense clonal expansion, and thus are likely to differentiate from effector cells after response.

Interestingly, we observed a very different response of CD4+ and CD8+ memory cells to the booster vaccination. It may be explained by differences in antigen presentation mechanisms: CD4+ T-cells may be activated well by antigen presenting cells phagocyting neutralized viral particles and presenting exogenous peptides on MHC-II complexes, while CD8+ memory cells can be more efficiently triggered by a productive viral infection resulting in the presentation of endogenously translated viral proteins on MHC-I. It was previously shown that the magnitude CD8 response depends on the viral load (Akondy et al., 2015).

It will be interesting to perform a similar study in donors vaccinated with YF backbone chimeric vaccines, where genes from other viruses substitute some of the YFV17D genes. It was shown that preexisting anti-YF immunity (Monath et al., 2002) does not affect the formation of neutralizing antibodies to the novel virus. This finding suggests that not only efficient reactivation of existing CD4 memory but also the formation of CD4 responses to novel epitopes is possible during the booster with slightly different antigen.

We found that, while the overall secondary response to the vaccine was much smaller both in terms of clonal diversity and cumulative frequency, a few clones still undergo strong clonal expansion. This may be indirect evidence for the programmed proliferation hypothesis (Moore et al., 2019) according to which a single encounter of a TCR with an antigen triggers several rounds of T-cell division. It was shown that the virus is undetectable in the peripheral blood after booster vaccination (Reinhardt et al., 1998), meaning that the amount of available antigen is much lower, and so is the number of encounters and responding clonotypes.

The transition between EM and EMRA phenotypes in CD8+ clones responding to yellow fever vaccine was previously measured using flow cytometry (Wieten et al., 2016; Fuertes Marraco et al., 2015). Here we confirm these reports with high-throughput sequencing, using TCR as a barcode to mark cells of the same clonal lineage. Furthermore, we identified two distinct cytotoxic phenotypes in NS4B-specific T-cells 18 months after primary immunization. It is not clear why the distribution of clonotypes between two these phenotypes was biased. Since we performed scRNAseq of clonotypes specific to the single antigen, these differences might be either the consequence of different TCR affinity or some phenotypic heterogeneity present in the precursor cells. Additional experiments at later timepoints would be required to estimate the longevity of these clonotypes.

To summarize, we show that vaccination with YFV17D leads to the recruitment of a diverse repertoire of T-cells, which is then available as immune memory for the secondary response years after the immunization. Even T-cells with the same epitope specificity show several distinct motifs in TCR and have different memory phenotypes. Such heterogeneity of cells might be crucial for individual immune response robustness, underlying cross-reactive responses to similar viruses, and the possibility to escape mutants, which could be tested directly in future studies. However, this diverse T-cell response is strongly focused on single HLA-A02 restricted epitope. An interesting question is how many distinct foci of response exist in the human population with a variety of HLA-types; and how this diversity of individual responses contribute to the defense from the infection at the population level. Systematic studies of donors with different genetic backgrounds and corresponding immunodominant epitope-specific repertoires will be able to address this question.

Materials and methods

Key resources table
Reagent typeDesignationSourceIdentifiersAdditional
information
AntibodyAnti-CD3-FITC (Mouse monoclonal)eBioscienceCAT# 11-0038-42FACS (5 ul per test)
Antibodyanti-CD45RA-eFluor450 (Mouse monoclonal)eBioscienceCAT# 48-0458-42FACS (5 ul per test)
Antibodyanti-CCR7-AlexaFluor647 (Rat monoclonal)BD PharmingenCAT# 560921FACS (5 ul per test)
Antibodyanti-CD95-PE (Mouse monoclonal)eBioscienceCAT# 12-0959-42FACS (5 ul per test)
Antibodyanti-CD3-eFluor450 (Mouse monoclonal)eBioscienceCAT# 48-0038-42FACS (5 ul per test)
OtherHLA-A*0201 (LLWNGPMAV) dextramerImmudexCAT# WB3584FACS (10 ul per test)
Commercial kitChromium Single Cell A Chip Kit10x GenomicsCAT# 1000009
Commercial kitChromium Next GEM 5’ Library and Gel Bead Kit10x GenomicsCAT# 1000014
Commercial kitChromium V(D)J Enrichment Kit, Human T Cell10x GenomicsCAT# 1000005
Commercial kitChromium Single Cell 5’ Library Construction Kit10x GenomicsCAT# 1000020
Commercial kitChromium i7 Multiplex Kit10x GenomicsCAT# 120262
Commercial kitDynabeads CD4 Positive Isolation KitInvitrogenCAT# 11331D
Commercial kitDynabeads CD8 Positive Isolation KitInvitrogenCAT# 11333D

Donors and blood samples

Request a detailed protocol

Blood samples were collected from two healthy donors (M1 male age 26, and P30 male age 39) on multiple timepoints before and after immunization with YFV17D vaccine. All donors gave written informed consent to participate in the study under the declaration of Helsinki. The blood was collected with informed consent in a certified diagnostics laboratory. The study was approved by the institutional review board (IRB) of Pirogov Russian National Research Medical University. HLA haplotypes of donors (Table 1) were determined by in-house RNA-based amplification and sequencing method.

Table 1
HLA-typing results for donors M1 and P30.
LocusM1P30
A02:01:01/24:02:0102:01:01/31:01:02
B15:01:01/39:01:0135:01:01/48:01:01
C03:04:01/12:03:0104:01:01/08:01:01
DQB102:01:01/03:02:0103:01:01/03:01:01
DRB103:01:01/04:01:0111:01:01/12:01:01
DRB302:02:0101:01:02/02:02:01
DRB401:03:01-

Isolation of T-cell subpopulations

View detailed protocol

We isolated PBMCs from the blood using standard Ficoll-Paque protocol. CD4 and CD8 fractions were isolated with CD4/CD8 Positive Selection Dynabeads Kits according to the manufacturer’s protocol. For isolation of memory subsets, we stained PBMCs with the mix of antibodies: anti-CD3-FITC (UCHT1, eBioscience), anti-CD45RA-eFluor450 (HI100, eBioscience), anti-CCR7-AlexaFluor647 (3D12, BD Pharmingen), anti-CD95-PE (DX2, eBioscience). Four subsets of cells were sorted into RLT buffer (Qiagen) on BD FACS Aria III: EM (CD3+CD45RA-CCR7-), EMRA (CD3+CD45RA+CCR7-), CM (CD3+CD45RA-CCR7+), Tscm (CD3+CD45RA+CCR7+CD95+). HLA-A02 dextramer loaded with the NS4B241-222 peptide (LLWNGPMAV) from YFV17D (Immudex) was used for epitope-specific T-cells isolation. Cells were stained with NS4B-dextramer-PE, anti-CD3-eFluor450 (UCHT1, eBioscience), and anti-CD8-FITC (SK1, eBioscience) according to the manufacturer’s protocol. RNA was isolated using standard TriZol protocol (for bulk PBMCs, CD4 and CD8, NS4B-specific and negative fractions) or RNAeasy Micro Kit (Qiagen) (for memory subsets). The amount of RNA was measured on Qubit 2.0 (Invitrogen). Information about all antibodies and commercial kits could be found in Key Resources Table.

Sample preparation for the single-cell gene expression and immune profiling

Request a detailed protocol

For 10x Genomics single-cell gene expression and immune profiling, we used PBMCs isolated from 60 ml of blood of donor M1 before the second immunization. PBMCs were stained with NS4B-dextramer-PE (Immudex) according to the manufacturer’s protocol. Additionally, cells were stained with anti-CD3-eFluor450 (eBioscience), and anti-CD8-FITC (eBioscience). Previous to FACS sorting procedure, we used propidium-iodide to mark dead cells. As the NS4B-specific cell frequency was very low (Figure 4—figure supplement 1C), we used anti-PE Ultra-pure MicroBeads (Miltenyi) for the enrichment. In brief, every milliard of PBMCs was incubated with 10 µl of magnetic beads for 15 min on ice. After a washing step with PBS 5% FCS, the cell suspension was applied on MS MACS Column (Miltenyi). Columns were washed three times with PBS 5% FCS and stained with propidium-iodide just before the FACS (FACS Aria II). This procedure resulted in a dramatic increase of NS4B-specific cell frequency in the sample (Figure 4—figure supplement 1C) and accordingly lead to reduced FACS procedure time. For single-cell immune profiling of bulk T-cell clonotypes from PBMCs, we stained the cells with anti-CD3-eFluor450 (Invitrogen) and propidium-iodide, thus selecting CD3 positive cells. Approximately 10,000 CD3+ cells were used for 10x Genomics VDJ T-cell receptor enrichment protocol.

High throughput T-cell repertoire sequencing

Request a detailed protocol

Libraries of TCR alpha and TCR beta chains were prepared as previously described (Pogorelyy et al., 2017). In brief, isolated RNA was used for cDNA synthesis with 5’RACE template switch technology to introduce universal primer binding site and Unique Molecular Identifiers (UMI) at the 5’ end of RNA molecules. Primers complementary to both TCR alpha and TCR beta constant segments were used for cDNA synthesis initiation. cDNA was amplified in two subsequent PCR steps. During the second PCR step, sample barcodes and sequence adapters were introduced to the libraries. Libraries for the fractions with low amount of cells (Figure 1—source data 1) were prepared using SMART-Seq v4 Ultra Low Input RNA kit (TakaraBio). Libraries were sequenced on Illumina platform HiSeq 2500 with 2×100 bp sequencing length or NovaSeq 2×150 bp sequencing length. Parallel single-cell alpha/beta TCR and 5’ gene expression sequencing was performed using 10x Genomics Kits (Chromium Single Cell A Chip Kit, Chromium Next GEM Single Cell 5’ Library and Gel Bead Kit, Chromium Single Cell V(D)J Enrichment Kit, Human T Cell, Chromium Single Cell 5’ Library Construction Kit, Chromium i7 Multiplex Kit) according to the manufacturer’s protocol. Libraries were sequenced on Illumina platform HiSeq 3000 with 2×150 bp sequencing length.

Repertoire data analysis

Raw data preprocessing

Request a detailed protocol

Raw repertoire sequencing data were preprocessed as described in Pogorelyy et al. (2017). Briefly, sequencing reads were demultiplexed and clustered by UMI with MIGEC software (Shugay et al., 2014). The alignment of genomic templates to the resulting consensus sequences was performed with MiXCR (Bolotin et al., 2015). Raw sequencing data obtained from RNAseq experiments were analyzed directly with MiXCR using default RNAseq analysis pipeline.

Identification of changed clonotypes by edgeR

Request a detailed protocol

To identify TCR alpha and TCR beta clonotypes that significantly expand after YF vaccination, we used the edgeR package (Robinson et al., 2010) as previously described (Pogorelyy et al., 2018). In brief, for each timepoint, we used two biological replicates of bulk PBMC. TMM-normalization and trended dispersion estimates were performed according to edgeR manual. We used an exact test based on the quantile-adjusted conditional likelihood (qCML) to identify clonotypes significantly expanded between pairs of timepoints. A clonotype with FDR adjusted p-value <0.01 (exact qCML-based test) was considered YF-responding if its log2-fold change estimate log2FC > 5 between any pairs of timepoints from 0 to the peak of the primary response (day 15). The usage of log2FC >5 threshold in addition to p-value threshold is important to filter statistically significant but small clonal expansions, which were previously shown to occur in healthy donors in the absense of vaccination on the timescale of one week, see Pogorelyy et al. (2018). The list of YF-responding clonotypes identified in alpha and beta TCR repertoires of donors M1 and P30 are in Figure 1—source data 3. CD4/CD8 in silico phenotyping was performed as suggested before (Pogorelyy et al., 2018): for each clone from bulk PBMC repertoire we assign CD4 phenotype if it is more abundant in the sequenced CD4 repertoire and vice versa. Over 98% of clonotypes were found exclusively in CD4 or CD8 compartment. However, a small group of clonotypes (1.4% for TCR alpha and 0.14% for TCR beta for day 15 timepoint of donor M1) was present in both compartments in comparable frequencies. These clonotypes have significantly higher TCR generative probabilities than others (p<0.001, Mann Whitney U-test) and thus are likely to arise from convergent recombination of the same TCR chain in both compartments.

To quantify the magnitude of the response on each timepoint we inferred the fraction of YF-responding cells as the proportion of all αβT-cells. To estimate this quantity from TCR repertoire data, for each susbset of interest (CD4+, CD8+, or NS4B-specific YF-responding clonotypes) we calculate the cumulative frequency of these clonotypes in TCR repertoire of bulk PBMCs in each timepoint.

Identification of YF-responding clonotypes by Principal Component Analysis (PCA)

Request a detailed protocol

We chose clonotypes that appeared in the top 1000 most abundant clonotypes at any timepoint after primary immunization. For these clonotypes, we made matrices of frequencies on all timepoints after primary immunization. Before applying PCA to these matrices, each value was normalized by dividing on maximum frequency for this clonotype. For cluster identification, we used hierarchical clustering with average linkage on euclidean distances between clonotypes. The number of clusters was set to 2. This analysis was performed for both alpha and beta chains of donor M1. For the twin donors (Pogorelyy et al., 2018), only replicate F1 was used for expanded clones identification.

Memory transition analysis

Request a detailed protocol

For this analysis, we used clonotypes that had at least 30 UMIs at day 45 after primary vaccination. The clonotype frequency in memory subset is multiplied by the number of cells obtained by FACS on this timepoint for this subset. Then adjusted frequencies are normalized across all subsets to get a partition of each TCR clonotypes across subsets. Obtained partitions were multiplied by the frequency of a clonotype in bulk at this timepoint to get the concentration of clonotypes with a particular memory phenotype in the bulk repertoire.

Computational decontamination of NS4B-specific repertoire

Request a detailed protocol

Since FACS sorting is not precise, TCR repertoires of the population of interest often contains abundant clonotypes from the bulk population. To obtain a list of NS4B-specific TCRs we took clonotypes that were enriched (at least 10 times) in the A02-NS4B-dextramer positive fraction compared to A02-NS4B-dextramer negative fraction. We also discarded TCR clonotypes that were more abundant in CD4 than CD8 subpopulation on day 0 (as only CD8 cells should bind to A02 which is a MHC I allele). Although ∼30% of resulting unique NS4B-specific clonotypes overlapped with the list of significantly expanded clonotypes identified with edgeR, they corresponded to ∼90% of NS4B-specific T-cells. See Figure 4—source data 1 for resulting list of NS4B-specific alpha and beta clonotypes for both donors.

Computational pairing of TCR alpha and TCR beta from bulk repertoires

Request a detailed protocol

For pair of clonal time traces we used a Euclidian distance between transformed frequencies:

D⁢(Cα,Cβ)=∑i(t⁢(Cα,i)-t⁢(Cβ,i))2,

where Cα,i and Cβ,i are the concentrations of an α and a β chain on the i-th timepoint. The transformation t of clone concentration C was chosen to address the overdispersion of frequencies at large concentrations (see Pogorelyy et al., 2018):

t⁢(Ci)=log10⁡(a+b⁢Ci+b⁢Ci),

where a=4.26×10-6 and b=3.09×10-3. To address possible systematic bias in expression between α and β chains in a clonotype, we introduce a log-fold shift λ in a trajectory with a quadratic penalty (µ=0.1):

Ds⁢(Cα,Cβ)=minλ⁡(D⁢(Cα,10λ⁢Cβ)+μ⁢λ2).

We calculated Ds distances between each pair of α and β clonotypes out of the 1000 most abundant ones in the bulk repertoires on day 15 post-vaccination. For each α clonotype, we picked the five closest β clonotypes as candidate pairings. As a benchmark, we used two single-cell TCR sequencing (scTCRseq) experiments using the 10x Genomics platform and obtained paired repertoires for samples of bulk T-cells (CD3+) and YF epitope-specific T-cells (CD8+NS4B-dextramer+). Note that these two samples are very different in their clonal time traces: NS4B-specific clones show very active response dynamics, expanding and contracting in the course of primary and booster immunization, while the CD3+ T-cell sample corresponds to the most abundant clones in the repertoire, which are largely stable between timepoints. A αβTCR clonotype from 10x Genomics experiment was considered correctly paired from bulk TCRseq data using the algorithm if the correct TCR beta was present among the five most probable TCR beta sequences predicted for the TCR alpha of this clonotype. Out of the 62 NS4B-specific clonotypes sampled in the 10x Genomics experiment, we were able to computationally identify 41 correct pairs from the bulk TCRseq data. Out of 26 CD3+ T-cell clonotypes, 20 were paired correctly.

Paired single-cell TCR sequencing

Request a detailed protocol

To investigate TCR chains pairing in YF-specific clonotypes, we performed single-cell immune profiling with 10x Genomics protocol. The analysis of the data with Cell Ranger 2.2.0 (10x Genomics) with default parameters resulted in 3244 cells corresponding to 986 clones. Many of these clones had multiple TRA/TRB chains and are likely to represent multimers of cells (Figure 5—figure supplement 3A). For further analysis, we chose only high-confident clones that had one TRA and one TRB chain and were present more than twice in the data. This procedure resulted in the list of ≈ 2000 cells corresponding to 164 TCR alpha/beta clones (Figure 4—source data 2).

TCR-pMHC complex modeling

Request a detailed protocol

Models for each paired alpha-beta TCRs from 10x Genomics data were constructed using the RepBuilder server (https://sysimm.org/rep_builder/) (Schritt et al., 2019), and then docked to HLA-A02-LLWNGPMAV complex using rosettaDock2 (https://www.rosettacommons.org/software) routine (Lyskov and Gray, 2008). 152 TCRs passed the modeling step. For each TCR we obtained 1000 decoys in docking simulations. The thirty best decoys (by interface score) were used to calculate a contact map with the bio3d R package (Grant et al., 2006). It was previously shown (Pierce and Weng, 2013), that some docking decoys exhibit binding modes which are not found in natural TCRs. In the analysis, we only used decoys in which the root mean squared deviation between the centers of mass of the alpha and beta chains in the decoys, and the centers of mass of these chains in at least one published HLA-A02-TCR complex from ATLAS database (Borrman et al., 2017), were less than 4 Å. The number of contacts to the peptide was averaged over decoys that passed the threshold. Only clonotypes with ≥5 of resulting filtered decoys were used for the analysis (see Figure 4—figure supplement 4A).

Single cell gene expression analysis

Request a detailed protocol

For single-cell gene expression analysis, we pre-processed the data with Cell Ranger 2.2.0 (10x Genomics). We used GRCh38-1.2.0 reference genome for the gene alignment. The resulting gene count matrix was analyzed with Seurat 3.0 package (Stuart et al., 2019; Butler et al., 2018). Cells that had fewer than 200 features detected were filtered out. We also filtered out features that were present in fewer than 3 cells and genes of TCR receptors (e.g., TRAV, TRAJ, TRBV, TRBJ), as they are the source of unwanted variation in the data (Figure 5—figure supplement 3B). Then a standard data pre-processing was performed to remove low-quality cells and cells multiplets. We filtered out cells that had more than 2700 features or more than 8% of mitochondrial genes (Figure 5—figure supplement 3C). Feature expression measurements for each cell were normalized using default log-normalization in the Seurat package. Following the manual’s suggestion, the 2000 most variable features were selected for further analysis. Prior to dimensionality reduction, data were scaled so that the mean expression was 0 and the variance equals to 1. The first 10 dimensions of PCA were used for cluster identification with the resolution parameter set to 0.4. To identify differentially expressed genes between clusters we used the MAST algorithm (Finak et al., 2015) implemented in the Seurat package. We only tested genes that were present in more than 25% of cells in any group and that had at least a 0.25 log fold difference between the two groups of cells. The resulting list of differentially expressed genes is reported in Figure 5—source data 1.

We performed a similar analysis to identify differentially expressed genes between clonotypes (rather than individual cells). We created a matrix containing the mean gene expressions over cells within each clonotype, and treated it like normal single-cell results. In this case, we did not filter multiplet cells (with a high number of features and a high percentage of mitochondrial genes), as all our ‘cells’ were indeed computational multimers. The rest of the analysis was performed in the same way. The list of differentially expressed genes between clusters of clonotypes is reported in Figure 5—source data 2. To check the results we shuffled cell barcodes between the clonotypes and repeated the analysis. All cells ended up in a single cluster for this random control.

Data availability

Sequencing data have been deposited in SRA under accession code PRJNA577794.

The following data sets were generated
    1. Minervina AA
    2. Pogorelyy MV
    3. Komech EA
    4. Karnaukhov VK
    5. Bacher P
    6. Rosati E
    7. Franke A
    8. Chudakov DM
    9. Mamedov IZ
    10. Lebedev YB
    11. Mora T
    12. Walczak AMW
    (2019) NCBI BioProject
    ID PRJNA577794. Comprehensive analysis of antiviral adaptive immunity formation and reactivation down to single cell level.
    https://www.ncbi.nlm.nih.gov/bioproject/PRJNA577794
The following previously published data sets were used
    1. Pogorelyy MV
    2. Minervina AA
    3. Touzel MP
    4. Sycheva AL
    5. Komech EA
    6. Kovalenko EI
    7. Karganova GG
    8. Egorov ES
    9. Komkov AY
    10. Chudakov DM
    11. Mamedov IZ
    12. Mora T
    13. Walczak AM
    14. Lebedev YB
    (2018) NCBI BioProject
    ID PRJNA493983. Precise tracking of vaccine-responding T-cell clones reveals convergent and personalized response in identical twins.
    https://www.ncbi.nlm.nih.gov/bioproject/PRJNA493983

References

Decision letter

  1. Satyajit Rath
    Senior and Reviewing Editor; Indian Institute of Science Education and Research (IISER), India
  2. Benjamin Chain
    Reviewer; University College London, United Kingdom
  3. Philip Bradley
    Reviewer; Fred Hutchinson Cancer Research Center, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Thank you for submitting your article "Comprehensive analysis of antiviral adaptive immunity formation and reactivation down to single-cell level" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by Satyajit Rath as the Senior and Reviewing Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Benjamin Chain (Reviewer #1); Philip Bradley (Reviewer #2).

The reviewers have discussed the reviews with one another and with the Reviewing Editor. Clearly, the study is interesting and worthy of publication, despite the limitation that it is based on only two individuals. However, the manuscript is poorly constructed, and hence very hard for a non-specialist to understand the message. Substantial revision and focusing are essential. Both reviewers have made these points in specific ways. Therefore, instead of a synthesised decision, both reviews are being provided in full below for ease of specific addressing to help you prepare a revised submission. The major issues raised by both reviewers must be substantively addressed in order to increase enthusiasm for publication.

Reviewer #1:

This paper contains a descriptive study of TCR sequencing and some single cell RNAseq on one primary, and two secondary responses to yellow fever vaccination in humans. There are definitely some interesting and intriguing observations here, although the study essentially remains a case-report, with many results supported by data from one single individual. How strongly these findings will generalise remains to be determined. The impact of the paper is weakened by a very loose and anecdotal writing style, and an attempt to include a large number of disparate findings. The paper could be much improved by more effort to clearly identify a small number of key messages and support them with data.

The authors focus on large clonal expansions (>32 times frequency post-primary mmunisation) and show that the secondary response contains many fewer large clone expansions than the primary response. Interestingly, a large proportion of clones expanded in primary responses can be resampled in the secondary. Almost as an aside they claim they can pair α and β chains by timecourse similarity alone (but see comment below). I found the PCA analysis and the clone pairing extraneous to the message of the paper, and merely a distraction.

They next demonstrate that most YF-specific clones switch from EM to EMRA phenotype post immunisation. This is an intriguing observation, given that the bulk population of EM is much larger than the EMRA one, but this observation is not followed up at all. Strangely there are no CD8 CM cells. Next, the authors analyse TCRs specific to the immunodominant epitope, which they show accounts for 60% of the total number of CD8 clones responding in primary immunisation. Next the authors analyse the TCR sequences of the immunodominant epitope, which they demonstrate fall into at least two distinct clusters of sequences. They support their claim that the two sets of TCRs bind via different modes by using some homology modelling, although this is buried in the supplementaries, and it is not clear how much weight to give to these results. Finally the authors use sc RNAseq to look at these clones, and show that they are heterogeneous, and specific clones are biased towards specific phenotypes.

Specific points:

I didn't follow exactly how the authors defined responding TCRs. They mention both a p value, and a cut-off of 32 fold enrichment. How are these related ?

The Y-axis in many panels is given as% YF_responding cells. This is very ambiguous – do they mean the fraction of cells which respond as a proportion of all cells ? Do all panels mean the same thing. This needs clarification.

In the section of TCR pairing, the authors state that they can correctly identify a large proportion of scTCRs using this approach on bulk sequences. But they also say, they pair each α with five nearest betas. So what do they mean by "correct paring": one of the five they computationally pair is the correct pair as determined by scTCR ? This is all very unclear.

Reviewer #2:

This manuscript reports a detailed characterization of T cell responses to Yellow Fever vaccination in two individuals, one over the course of primary and secondary (after 18 months) vaccination, and one a secondary vaccination 30 years after initial vaccination. The main strengths of the study are (1) the detailed, longitudinal characterization of the TCR repertoires (unpaired α and β chain) in bulk PBMCs as well as sorted T cell subsets (>15 time points over the two individuals), (2) bulk α/β and paired TCRseq data for T cells that bind an immunodominant A*02:01-restricted YFV epitope (NS4B), (3) paired scTCRseq and scRNAseq data for several thousand T cells positive for NS4B from one of the donors prior to their secondary vaccination, and (4) new and interesting analytical approaches for TCR repertoire analysis (time-course clustering, computational α/β pairing, scRNAseq profiles of T cell clonotypes). In my view the main weaknesses are the small number of donors and the mostly confirmatory nature of the biological results in light of the existing literature on YFV vaccination. Nonetheless I am enthusiastic about the study: the repertoire and RNAseq data will benefit others working in this area, and the new methods look like they will have a wide range of applications. I have a few questions/clarifications/potential typos, as detailed below.

– Clonotypes were assigned CD4/CD8 according based on the subset in which they had the highest frequency: did any clones span both the CD4+CD8 compartments, robustly (in terms of UMI counts)? Would you attribute these to sorting noise?

– My read of the paper is that for donor M1, YFV-responding clonotypes were defined in terms of their counts in the primary response. These clonotypes were then used to assess the magnitude of the secondary response. This begs the question of whether any new clonotypes were recruited into the secondary response. Are any new YFV-responding clonotypes identified by edgeR in the secondary response? Or are the frequency changes just too low to call new responding clones?

– How many of the NS4B-specific clones were called as YFV-responding by edgeR? Were there "genuine" NS4B-binders that did not expand? And did they have any special features? If so, were TRAV12 and TRAV27 NS4B+ clonotypes equally likely to be called as responding by edgeR?

– "Overall the dynamics and the magnitude of this response was very similar to the response to the booster vaccination after 18 months we observed in donor M1", maybe this is splitting hairs but to me the P30 response looks a little slower. Is this true/significant?

– I'm a little confused by the memory subset sorting. Where would plain old T effector cells show up here? How would they be distinguished from the memory subsets? Is there also CD45RO data? Or is there an assumption about composition based on the timepoints analyzed?

– "Interestingly, we found that the largest YF-specific CD8+ clones did not expand in response to the booster vaccine. Instead, the most expanded clonotypes were rare prior to the booster immunization (Figure 2—figure supplement 3A).", I think this should be "Figure 2—figure supplement 3B" if you are really talking about CD8 here.

– "While for most CD8+ clonotypes in the total repertoire EM/EMRA phenotypes were stable between day 15 and day 45 (Figure 3B, and Figure 3—figure supplement 2A, C), the distribution of CD8+ YF-responding clones between memory subsets was significantly shifted towards the EMRA phenotype (Figure 3C)." I don't see "CD8" in Figure 3B/C or Figure 3—figure supplement 2 (figure or legends), are there labels missing? Or are Figure 3C/Figure 3—figure supplement 2 actually PBMC?

– Figure 3A, legend, typo: "EM (CD45RA-CCR7+)"

– "Prior to dimensionality reduction, data were scaled so that the mean expression was 1 and the variance equals to 1", I think maybe this should say "mean expression was 0"?

https://doi.org/10.7554/eLife.53704.sa1

Author response

Reviewer #1:

This paper contains a descriptive study of TCR sequencing and some single cell RNAseq on one primary, and two secondary responses to yellow fever vaccination in humans. There are definitely some interesting and intriguing observations here, although the study essentially remains a case-report, with many results supported by data from one single individual. How strongly these findings will generalise remains to be determined. The impact of the paper is weakened by a very loose and anecdotal writing style, and an attempt to include a large number of disparate findings. The paper could be much improved by more effort to clearly identify a small number of key messages and support them with data.

The authors focus on large clonal expansions (>32 times frequency post-primary mmunisation) and show that the secondary response contains many fewer large clone expansions than the primary response. Interestingly, a large proportion of clones expanded in primary responses can be resampled in the secondary. Almost as an aside they claim they can pair α and β chains by timecourse similarity alone (but see comment below). I found the PCA analysis and the clone pairing extraneous to the message of the paper, and merely a distraction.

We agree that these results are supplementary to the key messages of the manuscript and could distract the reader. At the same time we agree with the second reviewer that both PCA time-course clustering and computational α-β pairing are potentially useful approaches to the analysis of the TCR repertoire data. We restructure the main text and move these results to figure supplements and Materials and methods.

1)We move Figure 2 CDE from the main text to Figure 2—figure supplement 4.

2)We removed “Identification of expanded clones by clustering of individual clonal trajectories” and “Computational α-β TCR pairing using individual clonal trajectories” from the main text.

3) We expanded the Materials and methods sections for α-β pairing (subsection “Computational pairing of TCR alpha and TCR beta from bulk repertoires”).

4) We added sentences to the main text to reference these methods (final paragraph in subsection “Diversity of clonal time traces in primary and secondary responses”).

5) We renamed the section “Time traces show a strong response of several clonotypes to booster vaccination” to “Diversity of clonal time traces in primary and secondary responses”

They next demonstrate that most YF-specific clones switch from EM to EMRA phenotype post immunisation. This is an intriguing observation, given that the bulk population of EM is much larger than the EMRA one, but this observation is not followed up at all. Strangely there are no CD8 CM cells.

We modified the text to show that CD8 CM and CD4 EMRA cells were also found in our analysis but with very low frequency. This agrees with previous results for these populations in peripheral blood (see Figure 2B in Thome et al., 2014 https://www.cell.com/cell/pdfExtended/S0092-8674(14)01314-2 and Figure 3 in Sathaliyawala et al., 2013 https://www.cell.com/action/showPdf?pii=S1074-7613%2812%2900521-3), although there was some variability between donors (see table S2 in Sathaliyawala et al., 2013 https://www.cell.com/cms/10.1016/j.immuni.2012.09.020/attachment/6aa930ec-8c28-41b9-8f5bb891297d211b/ mmc1.pdf).

We modified the main text (paragraph three in subsection “TCR sequencing shows the transition of clonotypes between memory subpopulations”) to clarify this point.

Next, the authors analyse TCRs specific to the immunodominant epitope, which they show accounts for 60% of the total number of CD8 clones responding in primary immunisation. Next the authors analyse the TCR sequences of the immunodominant epitope, which they demonstrate fall into at least two distinct clusters of sequences. They support their claim that the two sets of TCRs bind via different modes by using some homology modelling, although this is buried in the supplementaries, and it is not clear how much weight to give to these results. Finally the authors use sc RNAseq to look at these clones, and show that they are heterogeneous, and specific clones are biased towards specific phenotypes.

Specific points:

I didn't follow exactly how the authors defined responding TCRs. They mention both a p value, and a cut-off of 32 fold enrichment. How are these related ?

We require clone to pass both thresholds on p-value and on log2FC estimate. This is done to filter out statistically significant yet small changes in the repertoire, which occur within a week in a healthy donor in the absence of vaccination, as we observed in our previous work on yellow fever vaccination (Pogorelyy et al., 2018).

We changed Materials and methods text to clarify our strategy (subsection “Identification of changed clonotypes by edgeR”).

The Y-axis in many panels is given as% YF_responding cells. This is very ambiguous – do they mean the fraction of cells which respond as a proportion of all cells ? Do all panels mean the same thing. This needs clarification.

On all figure panels Y-axis indeed mean “the fraction of cells which respond as a proportion of all cells”, which was computed as cumulative frequency of all YF-responding clonotypes (or a subset of YF-responding clonotypes with specific features such as CD4+/CD8+/Dextramer+) in total TCR repertoire of PBMCs on different timepoints.

To explain how these numbers were calculated we added a sentence in Materials and methods section describing the procedure (subsection “Identification of changed clonotypes by edgeR”).

We modified the main text (paragraph two of subsection “Secondary T-cell response to the YFV17D vaccine is weaker but faster than the primary response”) to clarify what fraction we report on figures.

We also modified captions of Figures 1C,1D, 2AB, 4A to clarify this point.

In the section of TCR pairing, the authors state that they can correctly identify a large proportion of scTCRs using this approach on bulk sequences. But they also say, they pair each α with five nearest betas. So what do they mean by "correct paring": one of the five they computationally pair is the correct pair as determined by scTCR ? This is all very unclear.

We consider a clonotype from scTCR sequencing correctly paired if its β is present among 5 predicted by the algorithm for its α. The reason for pairing α to multiple betas is less diversity in TCR α chain, leading to more frequent convergent recombination events in the TCRalpha repertoire, and thus merging of α chains from many distinct clones into the single clonotype in the bulk RepSeq dataset.

We expanded Materials and methods section to explain how we defined correct pairing (subsection “Computational pairing of TCR alpha and TCR beta from bulk repertoires”).

Reviewer #2:

This manuscript reports a detailed characterization of T cell responses to Yellow Fever vaccination in two individuals, one over the course of primary and secondary (after 18 months) vaccination, and one a secondary vaccination 30 years after initial vaccination. The main strengths of the study are (1) the detailed, longitudinal characterization of the TCR repertoires (unpaired α and β chain) in bulk PBMCs as well as sorted T cell subsets (>15 time points over the two individuals), (2) bulk α/β and paired TCRseq data for T cells that bind an immunodominant A*02:01-restricted YFV epitope (NS4B), (3) paired scTCRseq and scRNAseq data for several thousand T cells positive for NS4B from one of the donors prior to their secondary vaccination, and (4) new and interesting analytical approaches for TCR repertoire analysis (time-course clustering, computational α/β pairing, scRNAseq profiles of T cell clonotypes). In my view the main weaknesses are the small number of donors and the mostly confirmatory nature of the biological results in light of the existing literature on YFV vaccination. Nonetheless I am enthusiastic about the study: the repertoire and RNAseq data will benefit others working in this area, and the new methods look like they will have a wide range of applications. I have a few questions/clarifications/potential typos, as detailed below.

– Clonotypes were assigned CD4/CD8 according based on the subset in which they had the highest frequency: did any clones span both the CD4+CD8 compartments, robustly (in terms of UMI counts)? Would you attribute these to sorting noise?

The vast majority of clonotypes are present in only one CD4/CD8 subset. But there are some clonotypes that are present in both CD4 and CD8 compartments (see Author response image 1). This clonotypes can be further divided into a very small group of clonotypes that are significantly enriched in one of the subsets (more than 10-fold, shown in blue on the figure) and clonotypes that are present in both subsets at similar frequencies (<10-fold difference in frequency, pink). We believe that the first ones represent contamination during isolation of CD4 and CD8 T cells, and our procedure provides correct phenotyping for this group.

The second group of clonotypes is much more interesting, because it consists of clonotypes with high generation probability (Pgen, as calculated by OLGA, Sethna et al., 2019), which convergently recombined in both subsets (see Author response image 1, right panel). This effect is much more prominent in α chains, due to their lower diversity. In our analysis, these clonotypes will be attributed to CD4/CD8 compartment randomly, but as they represent only a tiny fraction of all clonotypes we deliberately ignored it.

Author response image 1
Download asset Open asset

We expanded Materials and methods section (subsection “Identification of changed clonotypes by edgeR”) to explain this point.

– My read of the paper is that for donor M1, YFV-responding clonotypes were defined in terms of their counts in the primary response. These clonotypes were then used to assess the magnitude of the secondary response. This begs the question of whether any new clonotypes were recruited into the secondary response. Are any new YFV-responding clonotypes identified by edgeR in the secondary response? Or are the frequency changes just too low to call new responding clones?

To address this question we applied edgeR to second vaccination timepoints of donor M1. We were able to identify 91 YF-responding TCR β clones using second immunization, 18 of which intersect with the ones identified during the primary response. We were able to identify 103 YF-responding TCR α clones using second immunization, 30 of which intersect with the ones identified during the primary response. However, the contribution of these 73 additional clonotypes to the magnitude of response was low (see new Figure 1—figure supplement 1), and some of them also clearly expanded to the first immunization, but were not called significant by edgeR due to their low peak abundance. To summarize, we have no evidence of recruitment of naive clones during secondary response and no strong novel clonal expansions.

We expanded the main text to include this new results (paragraph three of subsection “Secondary T-cell response to the YFV17D vaccine is weaker but faster than the primary response”) and added a new figure to visualize them (see Figure 1—figure supplement 1)

– How many of the NS4B-specific clones were called as YFV-responding by edgeR? Were there "genuine" NS4B-binders that did not expand? And did they have any special features? If so, were TRAV12 and TRAV27 NS4B+ clonotypes equally likely to be called as responding by edgeR?

About 30% of NS4B-specific clones were also called as YF-responding by edgeR (537 out of 1790 in TCRbeta; 525 out of 1983 in TCRalpha). However these 30% clones correspond to 90% of all NS4B-cells (and thus non-overlapping clones have small abundances).

Interestingly, before the computational decontamination procedure we observed a small group of 60 CD4+ clonotypes and this group extensively used TRAV12-2 segment (26% of clones vs 2.6% in bulk subpopulation). This suggests that TRAV12-2 indeed may drive non-specific dextramer staining, even in CD4+ cells, which should not recognise pMHC class I dextramers.

After the decontamination procedure we used to filter out sorting noise (including CD4+ genuine binders mentioned before, see subsection “Computational decontamination of NS4B-specific repertoire”), we see almost no binders which do not expand between days 0 and 15 (however some appear at day 15 at low frequency and thus are not called significant by edgeR).

We expanded Materials and methods text to include these results (subsection “Computational decontamination of NS4B-specific repertoire”).

– "Overall the dynamics and the magnitude of this response was very similar to the response to the booster vaccination after 18 months we observed in donor M1", maybe this is splitting hairs but to me the P30 response looks a little slower. Is this true/significant?

This is true, the secondary T-cell response of M1 seems to peak on day 5 (Figure 1C) vs day 10 for donor P30. However, on Figure 1D we see that this effect is likely driven by a more active and rapid CD4 response in M1, while CD8 response also peaks on day 10 and thus dynamics is more similar to what we see in P30. As a speculation, delayed secondary CD4+ response in P30 may be explained by lower initial level of YF-responding T cells, but more donors would be needed to support this claim. Another interesting speculation is that some responding clones in P30 are actually newly recruited naive clonotypes, which take more time to expand. Such naïve clones might be produced after the first vaccination of P30 in childhood, and thus did not have a chance to participate in the primary response. Experiments on additional donors with deeper naive and memory repertoire sequencing on every time point could support or disprove this hypothesis.

– I'm a little confused by the memory subset sorting. Where would plain old T effector cells show up here? How would they be distinguished from the memory subsets? Is there also CD45RO data? Or is there an assumption about composition based on the timepoints analyzed?

We define CM, EM and EMRA subpopulations by CCR7 and CD45RA expression (similarly to Fuertes Maracco et al., 2015), but it is known that CD45RO expression is the opposite of CD45RA expression (as an example see Sathaliyawala et al., 2013, Figure 2A https://www.cell.com/action/showPdf?pii=S1074-7613%2812%2900521-3), so we expect it CD45RA-negative EM and CM compartments.

Term “effector” doesn’t have a well defined set of markers in humans and could mean either (1) cells that were recently primed by antigen or (2) cells that have effector functions. This issue is discussed in the Review by Appay et al., see https://onlinelibrary.wiley.com/doi/full/10.1002/cyto.a.20643). We can say that “effector cells” (defined as cells with effector functions) could be found in both EM and EMRA compartments, and additional markers would be necessary to distinguish short-lived effector cells (our results suggest that the majority of them correspond to EM, see subsection “TCR sequencing shows the transition of clonotypes between memory subpopulations”).

We added references to the review on T cell memory population and the study we used as an example for memory subpopulation partition.

– "Interestingly, we found that the largest YF-specific CD8+ clones did not expand in response to the booster vaccine. Instead, the most expanded clonotypes were rare prior to the booster immunization (Figure 2—figure supplement 3A).", I think this should be "Figure 2—figure supplement 3B" if you are really talking about CD8 here.

We now fixed this mistake.

– "While for most CD8+ clonotypes in the total repertoire EM/EMRA phenotypes were stable between day 15 and day 45 (Figure 3B, and Figure 3—figure supplement 2A, C), the distribution of CD8+ YF-responding clones between memory subsets was significantly shifted towards the EMRA phenotype (Figure 3C)." I don't see "CD8" in Figure 3B/C or Figure 3—figure supplement 2 (figure or legends), are there labels missing? Or are Figure 3C/Figure 3—figure supplement 2 actually PBMC?

This typo is now fixed.

– Figure 3A, legend, typo: "EM (CD45RA-CCR7+)"

This typo is now fixed.

– "Prior to dimensionality reduction, data were scaled so that the mean expression was 1 and the variance equals to 1", I think maybe this should say "mean expression was 0"?

This typo is now fixed.

https://doi.org/10.7554/eLife.53704.sa2

Article and author information

Author details

  1. Anastasia A Minervina

    Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
    Contribution
    Conceptualization, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    aminervina@mail.ru
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9884-6351

  2. Mikhail V Pogorelyy

    1. Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
    2. Center for Precision Genome Editing and Genetic Technologies for Biomedicine, Pirogov Russian National Research Medical University, Moscow, Russian Federation
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    m.pogorely@gmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0773-1204

  3. Ekaterina A Komech

    1. Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
    2. Center for Precision Genome Editing and Genetic Technologies for Biomedicine, Pirogov Russian National Research Medical University, Moscow, Russian Federation
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared

  4. Vadim K Karnaukhov

    Center of Life Sciences, Skoltech, Moscow, Russian Federation
    Contribution
    Formal analysis, Investigation
    Competing interests
    No competing interests declared

  5. Petra Bacher

    Institute of Immunology, Kiel University, Kiel, Germany
    Contribution
    Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared

  6. Elisa Rosati

    Institute of Clinical Molecular Biology, Kiel University, Kiel, Germany
    Contribution
    Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2635-6422

  7. Andre Franke

    Institute of Clinical Molecular Biology, Kiel University, Kiel, Germany
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Project administration
    Competing interests
    No competing interests declared

  8. Dmitriy M Chudakov

    1. Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
    2. Center for Precision Genome Editing and Genetic Technologies for Biomedicine, Pirogov Russian National Research Medical University, Moscow, Russian Federation
    3. Center of Life Sciences, Skoltech, Moscow, Russian Federation
    4. Masaryk University, Central European Institute of Technology, Brno, Czech Republic
    Contribution
    Conceptualization, Resources, Funding acquisition, Project administration, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0430-790X

  9. Ilgar Z Mamedov

    1. Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
    2. Masaryk University, Central European Institute of Technology, Brno, Czech Republic
    3. V.I. Kulakov National Medical Research Center for Obstetrics, Gynecology and Perinatology, Moscow, Russian Federation
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Project administration
    Competing interests
    No competing interests declared

  10. Yuri B Lebedev

    1. Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
    2. Moscow State University, Moscow, Russian Federation
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Methodology, Project administration, Writing - review and editing
    Contributed equally with
    Thierry Mora and Aleksandra M Walczak
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4554-4733

  11. Thierry Mora

    Laboratoire de physique de l’École normale supérieure, ENS, PSL, Sorbonne Université, Université de Paris, and CNRS, Paris, France
    Contribution
    Conceptualization, Formal analysis, Supervision, Investigation, Methodology, Writing - original draft, Writing - review and editing
    Contributed equally with
    Yuri B Lebedev and Aleksandra M Walczak
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-5456-9361

  12. Aleksandra M Walczak

    Laboratoire de physique de l’École normale supérieure, ENS, PSL, Sorbonne Université, Université de Paris, and CNRS, Paris, France
    Contribution
    Conceptualization, Formal analysis, Supervision, Investigation, Methodology, Writing - original draft, Writing - review and editing
    Contributed equally with
    Yuri B Lebedev and Thierry Mora
    Competing interests
    Senior editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2686-5702

Funding

European Research Council (Consolidator Grant no 724208)

  • Thierry Mora
  • Aleksandra M Walczak

Russian Science Foundation (15-15-00178)

  • Anastasia A Minervina
  • Mikhail V Pogorelyy
  • Ekaterina A Komech
  • Vadim K Karnaukhov
  • Yuri B Lebedev

Russian Foundation for Basic Research (18-29-09132)

  • Ilgar Z Mamedov

Russian Foundation for Basic Research (19-54-12011)

  • Ilgar Z Mamedov

Deutsche Forschungsgemeinschaft (Exc2167)

  • Petra Bacher
  • Elisa Rosati
  • Andre Franke

Deutsche Forschungsgemeinschaft (4096610003)

  • Elisa Rosati

Ministry of Science and Higher Education of the Russian Federation (075-15-2019-1789)

  • Dmitriy M Chudakov

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank JC Crawford and PG Thomas for assistance with TCRdist software and for helpful discussions. This work was funded by the European Research Council Consolidator Grant No. 724208 and RSF 15-15-00178. IZM was supported by RFBR 18-29-09132 and 19-54-12011. PB, ER and AF were supported by the Deutsche Forschungsgemeinschaft (DFG) through the Cluster of excellence Precision Medicine in Chronic Inflammation (Exc2167). ER was partially supported by DFG 4096610003. DMC was supported by grant 075-15-2019-1789 from the Ministry of Science and Higher Education of the Russian Federation to the Center for Precision Genome Editing and Genetic Technologies for Biomedicine under Federal Research Programme for Genetic Technologies Development for 2019–2027.

Ethics

Human subjects: All donors gave written informed consent to participate in the study under the declaration of Helsinki. The blood was collected with informed consent in a certified diagnostics laboratory. The experimental protocol was approved by the Ethical Committee of the Pirogov Russian National Research Medical University, Russia (FLU0108, granted January 29, 2016).

Senior and Reviewing Editor

  1. Satyajit Rath, Indian Institute of Science Education and Research (IISER), India

Reviewers

  1. Benjamin Chain, University College London, United Kingdom
  2. Philip Bradley, Fred Hutchinson Cancer Research Center, United States

Publication history

  1. Received: November 18, 2019
  2. Accepted: February 21, 2020
  3. Accepted Manuscript published: February 21, 2020 (version 1)
  4. Version of Record published: March 6, 2020 (version 2)

Copyright

© 2020, Minervina et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 3,296
    Page views
  • 505
    Downloads
  • 27
    Citations

Article citation count generated by polling the highest count across the following sources: Scopus, Crossref, PubMed Central.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Anastasia A Minervina
  2. Mikhail V Pogorelyy
  3. Ekaterina A Komech
  4. Vadim K Karnaukhov
  5. Petra Bacher
  6. Elisa Rosati
  7. Andre Franke
  8. Dmitriy M Chudakov
  9. Ilgar Z Mamedov
  10. Yuri B Lebedev
  11. Thierry Mora
  12. Aleksandra M Walczak
(2020)
Primary and secondary anti-viral response captured by the dynamics and phenotype of individual T cell clones
eLife 9:e53704.
https://doi.org/10.7554/eLife.53704

Categories and tags

Research organism

Further reading

    1. Computational and Systems Biology

    Resource allocation accounts for the large variability of rate-yield phenotypes across bacterial strains

    Valentina Baldazzi, Delphine Ropers ... Hidde de Jong
    Research Article

    Different strains of a microorganism growing in the same environment display a wide variety of growth rates and growth yields. We developed a coarse-grained model to test the hypothesis that different resource allocation strategies, corresponding to different compositions of the proteome, can account for the observed rate-yield variability. The model predictions were verified by means of a database of hundreds of published rate-yield and uptake-secretion phenotypes of Escherichia coli strains grown in standard laboratory conditions. We found a very good quantitative agreement between the range of predicted and observed growth rates, growth yields, and glucose uptake and acetate secretion rates. These results support the hypothesis that resource allocation is a major explanatory factor of the observed variability of growth rates and growth yields across different bacterial strains. An interesting prediction of our model, supported by the experimental data, is that high growth rates are not necessarily accompanied by low growth yields. The resource allocation strategies enabling high-rate, high-yield growth of E. coli lead to a higher saturation of enzymes and ribosomes, and thus to a more efficient utilization of proteomic resources. Our model thus contributes to a fundamental understanding of the quantitative relationship between rate and yield in E. coli and other microorganisms. It may also be useful for the rapid screening of strains in metabolic engineering and synthetic biology.

    1. Computational and Systems Biology
    2. Neuroscience

    Contrasting action and posture coding with hierarchical deep neural network models of proprioception

    Kai J Sandbrink, Pranav Mamidanna ... Alexander Mathis
    Research Article

    Biological motor control is versatile, efficient, and depends on proprioceptive feedback. Muscles are flexible and undergo continuous changes, requiring distributed adaptive control mechanisms that continuously account for the body's state. The canonical role of proprioception is representing the body state. We hypothesize that the proprioceptive system could also be critical for high-level tasks such as action recognition. To test this theory, we pursued a task-driven modeling approach, which allowed us to isolate the study of proprioception. We generated a large synthetic dataset of human arm trajectories tracing characters of the Latin alphabet in 3D space, together with muscle activities obtained from a musculoskeletal model and model-based muscle spindle activity. Next, we compared two classes of tasks: trajectory decoding and action recognition, which allowed us to train hierarchical models to decode either the position and velocity of the end-effector of one's posture or the character (action) identity from the spindle firing patterns. We found that artificial neural networks could robustly solve both tasks, and the networks'units show tuning properties similar to neurons in the primate somatosensory cortex and the brainstem. Remarkably, we found uniformly distributed directional selective units only with the action-recognition-trained models and not the trajectory-decoding-trained models. This suggests that proprioceptive encoding is additionally associated with higher-level functions such as action recognition and therefore provides new, experimentally testable hypotheses of how proprioception aids in adaptive motor control.

    1. Computational and Systems Biology

    Applying causal discovery to single-cell analyses using CausalCell

    Yujian Wen, Jielong Huang ... Hao Zhu
    Tools and Resources Updated

    Correlation between objects is prone to occur coincidentally, and exploring correlation or association in most situations does not answer scientific questions rich in causality. Causal discovery (also called causal inference) infers causal interactions between objects from observational data. Reported causal discovery methods and single-cell datasets make applying causal discovery to single cells a promising direction. However, evaluating and choosing causal discovery methods and developing and performing proper workflow remain challenges. We report the workflow and platform CausalCell (http://www.gaemons.net/causalcell/causalDiscovery/) for performing single-cell causal discovery. The workflow/platform is developed upon benchmarking four kinds of causal discovery methods and is examined by analyzing multiple single-cell RNA-sequencing (scRNA-seq) datasets. Our results suggest that different situations need different methods and the constraint-based PC algorithm with kernel-based conditional independence tests work best in most situations. Related issues are discussed and tips for best practices are given. Inferred causal interactions in single cells provide valuable clues for investigating molecular interactions and gene regulations, identifying critical diagnostic and therapeutic targets, and designing experimental and clinical interventions.