Characteristics of T-cell receptor repertoire of stem cell-like memory CD4+ T cells

Datasets

In this study, we conducted analyses on high-throughput TCR repertoire datasets of CD4⁺ T cell subsets. Gomez-Tourino et al. (2017) used a stringent strategy to sort Tn (CD3⁺CD4⁺CD45RO⁻CD27⁺CCR7⁺CD95⁻), Tscm (CD3⁺CD4⁺CD45RO⁻CD27⁺ CCR7⁺CD95⁺), and Tm (CD3⁺CD4⁺CD45RO⁺CD27⁺) from eight healthy subjects (HD) and eight T1D patients by fluorescence-activated cell sorting (FACS). Then RNA was extracted and sequenced in parallel. The sequence data is immuneACCESS format (https://clients.adaptivebiotech.com/pub/peakman-2017-naturecommunications). We examined the number of clonotypes in T1D and HD. As shown in Fig. S1, no significant difference was presented between T1D and HD in any subset.

To generating a model to identify public clonotypes which can occur in more than two individuals, we used datasets from Emerson et al. (2017) for training and testing a support vector machine (SVM) model. The dataset includes data of two cohorts, and can be found at https://clients.adaptivebiotech.com/pub/emerson-2017-natgen.

Statistical analysis and plots

Statistical analyses were performed with R. The paired Willcox-ranked test was used to examine the difference between two groups. The Kruskal-Wallis rank-sum test was used to examine the differences among multiple groups, and then Nemenyi test was used for multiple comparisons. The p values of multiple tests were corrected by false discovery rate (FDR) method. A test with a p value < 0.05 was considered as a significance. The Spearman correlation method was used to examine the correlation between samples of two groups. Graphics were generated with R package ggplot2. Principal component analysis (PCA) was conducted with R package forcats. R package readr, dplyr and tidyr were used for statistics.

Definition of a clonotype

A clonotype was defined as the amino acid sequence identity of the TRB CDR3 region.

Determination of diversity

Renyi entropy was used in our study to evaluate the diversity with alpha value from 0 to 20. When alpha increases, clonotypes with a higher frequency will have a greater influence on the entropy. When alpha equal 0, the Renyi entropy is the logarithm of the number of clonotypes; When alpha is 1, Renyi entropy tends to the Shannon entropy. When the alpha approaches infinity, the Renyi entropy is determined by the most frequent clonotype, where a lower frequency of the most frequent clonotype will generate a higher Renyi entropy index.

Renyi entropy formula is ${\displaystyle H=\frac{1}{1-\alpha }ln\left(\sum _{i=1}^n{f}_i^{\alpha }\right)}$

Shannon diversity index formula is ${\displaystyle H=-\sum _{i=1}^R{p}_{i}ln{f}_i}$ where H is the diversity index, α is the alpha value, n is the total number of clonotypes, and f_i is the frequency of the ith clonotype.

Hierarchical clustering

We performed ‘complete linkage’ clustering algorithm on the correlation matrix, and visualized dendrograms using pheatmap, a R package [50]. The Euclidean distance was used as a distance metric. Pearson correlation method was used to measure the correlation of TRB repertoire structures.

SVM model for identifying public and private clonotypes

SVM analysis was performed using kernel-based analysis of biological sequences with the R package KeBABS (Palme, Hochreiter & Bodenhofer, 2015). Amino acid sequence of clonotypes was split into features with length k = 3. A cost parameter C = 100 was used for the misclassification of a sequence. A total of 320,000 public and 320,000 private clonotypes were randomly sampled from the total set, and then were split into training (80%) and test (20%) sets. SVM training was performed on the training set, and class prediction was performed on the test set. Prediction accuracy of classification was qualified by calculating $BACC=\frac{1}{2}\times \left(spec+sens\right)$, where specificity was calculated as $spec=\frac{TN}{TN+FP}$, and sensitivity was defined as $sens=\frac{TP}{TP+FN}$, (where TN = true negative, FP = false negative, TP = true positive and FN = false negative). The area under the receiver operating characteristic curve (AUC) was calculated, where the AUC = 0.5 means a random classification (BACC = 50%), and AUC = 1 means a perfect classification (BACC = 100%).

The dataset (Emerson et al., 2017) includes two cohorts: cohort 1 includes 666 individuals, and cohort 2 includes 116 individuals. We termed clonotypes occurred in no less than two individuals as public clones and ones occurred in only one individual as private clonotypes. Data of cohort 1 was used as a training set. For this dataset of cohort 1, we randomly sampled 20,000, 40,000, 80,000, 160,000 and 320,000 public clonotypes and an equal number of private clonotypes, train a SVM model, and tested the model by cross validation (Fig. S2). When sampling more than 160,000 public clonotypes, the increase in the sample size had limited improvement in model accuracy. To save computing source, we trained a model on the 320,000 public clonotypes and 320,000 private clonotypes. The model was then validated with data of cohort2. In order to further increase prediction accuracy, we tested two thresholds for public clonotype definition: (1) definition of public clonotypes occurred in at least three subjects; (2) definition of public clonotypes occurred in at least two subjects. Comparing to the second threshold, the definition of public clonotypes occurred in at least three subjects elevated the prediction accuracy from 82% to 88%. We therefore defined public clonotypes occurred in at least three subjects.

Prediction of epitope specificity of clonotypes

GLIPH2 (Huang et al., 2020) is a robust tool to predict the cluster of clonotypes targeting the same epitope. Here, we used this method to cluster clonotypes recognizing the same antigens. The reference of CD4⁺ T clonotypes, the clonotypes gene usage and length distribution of CDR3 were included in ref_CD4.txt, ref_V_CD4.txt and ref_L_CD4.txt downloaded from the official website of GLIPH2 (http://50.255.35.37:8080/). A filter with a high stringency (Fisher_score < 0.0001, number of subjects >3, and number_unique_cdr3 ≥ 3) was used for identify the number of potential antigens in each sample.

Tscm and Tm had different TRB repertoire structures

To characterize immune repertoire clonal structure, we used Renyi entropy (Greiff et al., 2015). The α-values represent weights, which means as α increases, higher frequency clonotypes are weighted more. For a given alpha value, a larger Renyi entropy means a more considerable diversity of the sample. Since each alpha value focuses on a different stretch of the immune repertoire, this method enabled the reliable capture of TRB repertoire clonal frequency distributions. Our results showed that as alpha value increased, the Renyi entropy of all memory subsets decreased. At all alpha values, the Renyi entropy of Tn was higher than that of both Tscm and Tm. At the alpha value from 0 to 2, the Renyi entropy of Tscm was less than that of Tm, while after alpha value 3, the Renyi entropy of Tscm was greater than that of Tm (Fig. 1A), which reflects that Tscm had a greater diversity than Tm in the abundant clonotypes. Since the Renyi entropy profiling can recover a large amount of immunodiagnostic fingerprints from TRB repertoire data, we used the Renyi entropy to classify the cell subsets with hierarchical clustering approach based on Pearson correlation (Greiff et al., 2015). The results showed that Tscm of 13 subjects were clustered together; Tm of 15 subjects were gathered; Tscm of only 3 subjects mixed in the cluster of Tm (Fig. 1B). This hierarchical clustering result suggested that Tscm and Tm had different TRB repertoire structures.

Top1000 clonotypes of Tscm and Tm used different genes

The gene usage of the TRB repertoire of memory T cells is heavily affected by antigen experience. We analyzed the gene usage to unveil the effects of antigen-experience on Tscm and Tm, respectively. Since the frequent clonotypes of memory T cells are expanded by chronic antigen stimulations, we also analyzed the gene usage of the top1000 abundant clonotypes. The PCA on the gene usage of the entire repertoire separated Tn, Tscm, and Tm from each other (Fig. 2A). Specially, the PCA on the genes of the top1000 clonotypes could achieve a better performance of classification (Fig. 2B). A further analysis on top1000 abundant clonotypes showed that Tscm and Tm differentially used 12 V-genes and 3 J-genes: the frequency of TRBV12-03, TRBV07-09, TRBV18-01, TRBV23-01 and TRBJ02-07 were less in Tm than in Tscm; the frequency of TRBV03-01, TRBV02-01, TRBV11-02, TRBV09-01, TRBV06-05, TRBV25-01, TRBV24-01, TRBV05-05, TRBJ01-02 and TRBJ02-02 was greater in Tm than in Tscm (Table S1). In further, we used Spearman correlation to quantify the similarity of the gene usage between Tscm and Tm. Because Tscm and Tm differentiate from Tn, we therefore used the correlation between Tn and Tscm as a contrast. For the entire TRB repertoire, the correlation between Tm and Tscm was similar to that between Tscm and Tn, but less than that between Tm and Tn (Fig. 2C); for the top 1000 abundant clonotypes, the correlation between Tscm and Tm was greater than that between Tscm and Tn (Fig. 2D). It suggests that Tscm and Tm had a large difference in the gene usage of the TRB repertoire, especially in the range of top1000 abundant clonotypes.

Tscm was different from Tm in CDR3 length distribution

The antigen experience has a selection on clonotypes which may change the distribution of CDR3 length. For the CDR3 of the total clonotypes, Tscm was significantly longer than Tn, and Tm as well (Fig. 3A). For the CDR3 of the top1000 abundant clonotypes, Tscm was obviously shorter than Tm, but longer than Tn. (Fig. 3A). According to the Spearman correlation analysis, Tscm had a similar length distribution to Tm rather than to Tn in the top1000 abundant clonotypes (Fig. 3B). Tscm and Tm also showed an increased correlation in gene usage, then we examined whether gene usage was related to the distribution of CDR3 length in Tscm and Tm. However, we did not find a high correlation (cor = 0.27) between the gene usage and distribution of CDR3 length in Tscm and Tm (Fig. 3C), suggesting that, between Tscm and Tm, the increased correlation of gene usage and elevated correlation of length distribution might be independent.

Tscm had special CDR3 sequence compositions

We examined the CDR3 sequence composition by decomposing kernels containing three amino acids. As showed by the PCA on the sequence composition of all clonotypes, Tn, Tm, and Tsm samples were partially separated (Fig. 4A). Exhibited by the PCA based on the k-mer of the top1000 abundant clonotypes, the samples of Tn, Tscm, and Tm were completely separated (Fig. 4B). These results suggested that Tscm and Tm had a great difference in the sequence composition of the top1000 clonotypes. We used Spearman correlation to quantify differences of sequence composition among subsets. For the total TRB repertoire, the correlation between Tm and Tscm was significantly weaker than that between Tn and Tm, and slightly weaker than that between Tscm and Tn (Fig. 4C). For the top1000 abundant clonotypes, the correlation between Tm and Tscm was significantly weaker than that between Tn and Tscm, and between Tn and Tm as well (Fig. 4D). Furthermore, to identify whether the correlation between Tscm and Tm reduced in the top1000 clonotypes, we randomly sampled 1,000 clonotypes from the entire repertoire as a contrast. Our results showed that the correlation coefficient of the top1000 clonotypes between Tscm and Tm was significantly lower than that of random subsamples. In contrast, the correlation of top1000 clonotypes sequence composition between Tn and Tscm, and that between Tn and Tm were significantly stronger than that of random samples between corresponding cell subsets (Fig. 4D). It suggested that Tscm and Tm were different in the sequence composition of the entire TRB repertoire, especially the top1000 abundant clonotypes.

Frequent clonotypes of Tscm were more public

Public clonotypes which are shared among subjects can be induced by gene recombination and antigen stimulation (Venturi et al., 2008). To evaluate the public clonotypes’ distribution in each sample, we trained a SVM model to identify public clonotypes with a large dataset (Greiff et al., 2017). This dataset includes two cohorts. The SVM model had a high prediction accuracy (BACC = 88% and AUC = 95%) in cohort 1 and presented a strong robustness when classifying public clonotypes on the data of cohort 2 (BACC = 82%) (Fig. S3). Using the SVM model, we identified public clonotypes from samples of Tn, Tscm, Tcm, and Tem. Since the public clonotypes may be unevenly distributed in different frequency ranges of the TRB repertoire, we calculated the percentage of public clonotypes for each frequency range. We ranked clonotypes by their frequency. Our result showed that as the ranks increased, the percentage of public clonotypes decreased, which was inconsistent with findings in global T cells (Soto et al., 2020). Notably, percentage of identified public top1000 clonotypes in Tscm were more than the percentage in Tm. (Fig. 5A). We also used the data of top1000 clonotypes of effector memory T cells and central memory T cells presented by Jiang et al. (2020), and found a similar trend that Tscm contained more public top1000 clonotypes than effector memory T cells (Tem) as well as central memory T cells (Tcm). (Fig. 5B). We further trained a SVM model to classify public clonotypes of Tn and Tm, and showed a prediction accuracy of 66%, which suggested that public clonotypes of Tn and public clonotypes of Tscm contained predictive high-dimensional features.

To further identify the antigen specificity of clonotypes in memory subsets, we used GLIPH2 to cluster clonotypes potentially targeting same antigens. A stringent cutoff was used to avoid potential mistakes when performed GLIPH2. In GLIPH2, a clonotype can be grouped in over one cluster. For top1000 abundant clonotypes, 1095 clusters were exhibited in Tn, 176 exhibited in Tscm, and 71 exhibited in Tm. For public clonotypes within top1000 clonotypes, we defined 541 clusters in Tn, 117 in Tscm, 67 in Tm (Tables S2 and S3).

Self-reactivated clonotypes expanded in Tscm rather than Tm in T1D patients

Since Tscm potentially play roles in infections and autoimmune diseases, we evaluated the pathogen and autoreactive clonotypes in Tscm and Tm, respectively. The clonotypes were annotated by a database including 1,885 pathogen-related clonotypes from VDJdb (Bagaev et al., 2020), 1,409 autoreactive clonotypes (Gomez-Tourino et al., 2017). Our results showed that autoreactive, private clonotypes could be found in Tscm of 11 subjects, and in Tm of 3 subjects; autoreactive, public clonotypes were shown in Tscm of 16 subjects, and in Tm of 16 subjects. The number of autoreactive clonotypes were similar in Tscm and Tm, however, there were more autoreactive public clonotypes than autoreactive private clonotypes in each sample (Fig. 6A). Considering clonotypes’ frequency in each subject, we showed that GAD-related clonotypes within Tscm became more frequent in T1D than in HD. In contrast, the frequency of GAD-related clonotypes in Tm was similar in T1D and HD (Fig. 6B). Meanwhile, we showed more pathogen-related public clonotypes than pathogen-related private clonotypes for Tscm and Tm (Fig. 6C). However, we did not find that pathogen-related clonotypes of either Tscm or Tm were more frequent in T1D (Fig. 6D). In conclusion, self-reactivated clonotypes expanded in Tscm of T1D patients.