We used previously published (PBMC) TCR repertoire sequencing data for both TCRα (23) and TCRβ (24) chains. For TCRβ, a sample of 10 CMV donor repertoires selected ad hoc were chosen from the HIP cohort of the Emerson et al. (24) dataset for TCRβ analysis (sample IDs in dataset from 1 to 10). For TCRα we used all bulk TCRα PBMC data (10 samples) available in Heikkilä et al. (23). Only clonotypes that are supported by a single read (singletons) were used in subsequent analysis (9). The TCRβ sample consisted of 1,147,250 TCR sequences and the TCRα sample of 1,582,774 ( Supplementary Table 1 ). Note that generation biases were not controlled, as the datasets comprised sequences from multiple individuals, which masked the effects of individual generation biases. Additionally, generation biases are imprinted in the selection process [see reference (16)].

Repertoires of sorted CD4+ and CD8+ naive (post-selection) T-cell sequencing were taken from Qi et al. (25); this data is available for TCRβ chain sequencing only. TCR repertoire DP and CD8+ (SP) thymocytes (pre-selection) were taken from a recent Quiniou et al. (10) study. Naive repertoires of all nine donors in this dataset were combined together, resulting in 1,346,776 clones for CD8+ naive cells and 1,599,217 clones for CD4+ naive cells ( Supplementary Table 1 ). Note that for these datasets no read count information was used (i.e. all clonotypes were assumed to be singletons) in order to avoid potential amplification biases and to make it compatible with other datasets.

A dataset containing PBMC TCRβ repertoire sequencing for three pairs of monozygotic twins was taken from “Precise tracking of vaccine-responding T cell clones reveals convergent and personalized response in identical twins” by Pogorelyy et al. (26). Pre-vaccination and day 0 repertoires sampled prior to treatment were used as biological replicates for each twin, and only singletons were included in the analysis. HLA alleles for donors are listed in Supplementary Table 2 . An additional dataset containing PBMC TCRα and β repertoires from three pairs of monozygotic twins was taken from “Distinctive properties of identical twins’ TCR repertoires revealed by high-throughput sequencing” by Zvyagin et al. (27). Sorted naive CD4+ T-cells TCRβ repertoires from two pairs of monozygotic twins were also obtained from “Functionally specialized human CD4+ T-cell subsets express physicochemically distinct TCRs” by Kasatskaya et al. (19). Summary statistics for datasets are reported in Supplementary Tables 1 , 3 , and 4 .

Pre-selection TCRα and β repertoires were simulated based on a theoretical probabilistic model of the V(D)J rearrangement process using OLGA software (v1.2.4) as described previously (8). The software was executed with default runtime parameters and model probabilities, random seed was set to 100, and a sample of 10⁷ random rearrangements was generated for each TCR chain.

Single-cell datasets with assigned cell types totaling 178307 cells from PBMC samples of 88 healthy patients were taken from Lindeboom et al. (28). Cell type annotation was performed using CellTypist, as specified in the original study (29). Cluster abundance on a particular cell type was tested using the Fisher exact test.

Basic features of the TCR sequence such as V/J gene usage, single amino acid frequencies, k-mer (k=3) frequencies, and physicochemical properties of CDR3 regions were carried out using in-house scripts as described previously (30). The set of k-mers was not filtered based on their relative position in CDR3. Kidera factors (key amino acid features that describe most variance in polypeptide physicochemical properties) (31), charge, and hydrophobicity were calculated using the “peptides” python package (v0.3.2). The most informative Kidera factors, Kideras 2, 4, 6, and 8, correspond to Side-chain size, Hydrophobicity, Partial specific volume, and Occurrence in the alpha region respectively. Kidera factor values were compared using T-test; Cohen’s d was used as an effect size estimate.

Kidera factors were not z-score normalized because the datasets compared with OLGA-generated data were obtained using the same protocol as the data used to train the OLGA model ( Supplementary Table 1 ) and all samples in the thymocytes dataset were obtained in a single batch. Detection of TCR CDR3 sequence motifs was performed using the TCRNET algorithm implemented in VDJtools (v1.2.1) as described previously (11). This method defines TCR sequences of interest as those that are placed in the more dense regions of a CDR3 sequence similarity graph compared to a control (typically produced assuming V(D)J rearrangement model with no selection pressure) dataset: the number of 1-hamming distance neighbors is compared to the expected number of neighbors adjusted for sample and control sizes to produce an enrichment score and a P-value based on Binomial approximation.

Note that in order to infer TCR sequence clusters that were depleted by negative thymic selection, we simply swapped “background” (control) and “foreground” (our sample of interest), i.e. we searched for TCRs enriched in pre-selection data compared to post-selection.

In order to produce a representative set of TCR clusters, we selected the top 10,000 neighbor-enriched CDR3 sequences based on enrichment P-value. Selected sequences were clustered by choosing connected components of the graph with edges connecting sequences that differ by a single amino acid substitution. Motifs for selected clusters were visualized using logomaker package (v0.8). The top five largest clusters were subsequently analyzed and numbered according to their cluster size rank. The number of clusters was selected ad-hoc.

The SoNNia (v0.2.3) model was additionally used to assess differences in amino acids’ occurrence probabilities, in particular their position in post- vs pre- selection repertoires (32). For calculating marginal probabilities of sequence features, we trained the SoNNia model using pairs of post-selection and pre-selection datasets (e.g. post-selection TCRβ and OLGA generated TCRβ data). Post-selection datasets were used as data for model inference and pre-selection datasets were used as data sampled from generative distribution. Each dataset was processed with methods from the “Processing” class and then passed to the SoNNia model with appropriate “pgen_model” parameters. Each model was trained for 50 epochs, with a batch size equal to 10⁴. Models were assessed with built-in plotting functions, as shown in the SoNNia tutorial (https://sonnia.readthedocs.io/en/latest/sonnia_tutorial.html).

In order to identify positively and negatively selected CDR3 clusters in the twins dataset, we subsampled each twin sample to 306,553 CDR3s (size of the smallest repertoire) and pooled all the samples together. We used ad hoc thresholds to select significantly enriched (log2 fold change > 2 and -log10 p > 12 for sample pool compared to simulated sequences as control) and depleted (log2 fold change > 1 and -log10 p > 12 for simulated sequences compared to sample pool as control) clusters after thymic selection clusters. Clusters containing more than 10 sequences were used for further analysis. Next, similarity between positively and negatively selected clusters was estimated by computing the Jensen-Shannon divergence between cluster frequencies defined as the number of clonotypes from a given cluster present in a given sample.

CDR3 loop structures for TCRs of interest were modeled using the TCRmodel web tool (33) and processed with Pymol (version 2.3.0). Our in-house “mir” software package was used to annotate the resulting PDB files (see (34)). As CDR3α is not known in most of our datasets, we used a generic CAGGSSNTGKLIF (TRAV27, TRAJ37) sequence that was the most commonly observed variant in the dataset from Heikkila et al. (23) as a dummy TCRα sequence. The TCR CDR3 backbone was visualized by applying PCA to Cα atom coordinates.

Experimental structures of TCR:pMHC complexes were obtained from VDJdb (35). We selected only records with available PDB IDs of human TCRs. The structures were required to harbor both TCRα and TCRβ along with a pMHC complex. Additionally, at least one CDR3β residue had to be within 5 Å of the peptide in a pMHC complex to ensure direct contact between the TCR β chain and the pMHC complex. A total of 154 structures in total were analyzed. Dihedral angles in these structures were calculated with Biopython (v. 1.85) Python package (36).

All code used in this study is available at https://github.com/LuppovDaniil/Thymic_selection_notebooks (Python version 3.11.5., R version 4.1.2).

AIRR-seq data for DP T-cells sorted from the thymus can be used to explore the initial space of rearranged TCRα and β sequences existing prior to positive and negative selection, similar to the work of Quiniou et al. (10). Recent studies demonstrate that a probabilistic model can accurately replicate the structure of V(D)J rearrangement space, generating TCR sequences with amino acid composition and frequencies resembling those produced in vivo (8). Here we use both Quiniou et al. and model datasets as a pre-selection repertoire.

There are several ways to acquire the TCR repertoires of post-selection T-cells that have not yet undergone strong antigen exposure and obtained memory phenotypes. One can either sort and sequence naive CD4+ and CD8+ T-cells as performed in Qi et al. (25) or use SP thymocytes as done in Quiniou et al. (10). Alternatively, one can select T-cell clonotypes detected only once (singletons) from unsorted PBMC AIRR-seq data, as they mostly represent naive T-cells (see Britanova et al. (9)). In the present study, we used the datasets from Qi et al. and Quiniou et al. and selected singletons from 10 samples chosen ad hoc from the Emerson et al. dataset as a TCR repertoire after selection but prior to any antigen exposure or subsequent expansion. The detailed description of these datasets and comparison analysis is given in the Method section.

Analysis of TCRβ CDR3 amino acid frequencies revealed a significant post-selection decrease in specific residues compared to those expected from a V(D)J rearrangement model ( Figure 1A ). Positively charged and physically large amino acids, such as arginine, histidine, and lysine, were likely reduced due to strong antigen binding or steric hindrance in antigen recognition (37), while Proline and Cysteine were negatively selected possibly due to their effects on TCR structure. These findings were generally confirmed in thymocyte data ( Supplementary Figure 1A ).

Analysis of 3-mer frequencies revealed that 3-mers with cysteine had less chance to survive selection ( Figure 1B ), consistent with observations at the single amino acid level ( Figure 1A ). 3-mers with the NX[S,T] motif, associated with N-glycosylation sites (38), are less likely to survive selection, with NPT and NPS being the least affected, consistent with the absence of glycosylation at these sites (39) ( Figure 1B ). We also analyzed sulfonation site motifs, which were described in Pospelova et al. work (40) for antibodies, such as DDD, DDY, YYY, and EEE. The tangible reduction effect was observed only for the DDD motif. Both glycosylation sites and Cysteine-linked effects and a lack of effect from putative sulfation sites were observed in the thymocytes data ( Supplementary Figure 1B ). Additionally, we assessed the positions of the glycosylation sites in the CDR3 sequences and discovered that, for CDR3β, most (85% and more for different datasets) glycosylation sites were located in the NDN segment, highlighting their somatic origin. For CDR3α, however, only 44-53% of these sites were located in the N segment.

In order to describe changes in physicochemical properties of TCRβ CDR3s, we harnessed Kidera factors (31). These factors represent the key physical properties of amino acids obtained by dimensionality reduction. We compared each VJ pair independently in order to lessen the bias caused by the choice of the V and J genes ( Figure 1C ). The Kidera factors negatively affected by selection were Kidera 2, which determines side chain size, Kidera 4, which is inversely associated with hydrophobicity, and Kidera 6, which determines partial specific volume. These results suggest that TCRβs with a physically small hydrophobic CDR3 have a greater chance of passing selection. The only factor that increased post-selection was Kidera 8 (occurrence in the α-helix structural region). Notably, the selection effect within a particular Kidera Factor was shared across all VJ pairs. Identical analysis of DP and SP thymocytes mirrored the above results, however, the effect for Kidera factor 6 was less pronounced ( Supplementary Figure 1C ).

The analysis of the physical properties of TCRβ CDR3 revealed a decrease in repertoire charge post-selection ( Figure 1D ), consistent with the observed reduction in frequencies of positively charged amino acids ( Figure 1A ), the winsorizing (resulting distribution having thinner tails) of lengths (too short and too long sequences both have less chance to pass selection), and the increase in hydrophobicity ( Figure 1F ). Except for length, which experienced shortening instead of winsorizing, the same effects were detected on thymocytes data ( Supplementary Figures 1D–F ).

Next, we identified functional clusters of TCRβ CDR3s enriched in naive cells compared to the modeled sample and vice versa (30). The largest enriched clusters in naive and model-generated samples are presented in Figure 1G . Model-enriched clusters can be interpreted as sequences which tend not to pass the selection and naive sample-enriched clusters are ones which are likely to survive thymic selection. For convenience, the terminology (enriched and depleted pre- and post-selection versus positively and negatively selected) used here and in the Supplementary Materials is explained in Figure 1H .

The majority of the post-selection enriched clusters contain poly-Glycine sequences in the middle of a TCRβ CDR3 ( Figure 1G ; Supplementary Figure 2 ), which is known to be one of the most flexible among the polypeptides chains (41). Additionally, k-mers analysis demonstrated that GGG 3-mer passed through the selection with unchanged frequency ( Figure 1B ; Supplementary Figure 1B ), confirming the above finding. CDR3β from pre-selection enriched clusters frequently contained Arginine and Proline (residues likely impacting CDR3 loop structure) and exhibited deviations in the CASS consensus sequence at the CDR3 start ( Figure 1G ; Supplementary Figure 2 ). Results of enriched cluster analysis for thymocytes data mostly resembled those obtained in the generated data ( Supplementary Figure 3 ).

In order to obtain pre- and post- selection probabilities of amino acid occurrence in a particular position in CDR3, we utilized the SoNNia software (32). Generally, we observed the nearly identical to enriched clusters analysis trends ( Supplementary Figures 4A, С ).

Gene usage analysis revealed a vague picture of little or no preference of selection to more frequent genes (or “rich get richer effect” (42)) on both generated and experimental data ( Supplementary Figures 5, 6 ).

Interestingly, while no “rich get richer effect” was observed for the gene frequencies, it was detected in terms of CDR3β generation probabilities (pgens) for DP and SP thymocytes calculated via the OLGA model (p < 10^-5, Mann Whitney test, differences of median log2 Pgen between SP and DP was 2.24, Supplementary Figure 1G ).

Analysis of TCRα CDR3s, conducted identically to the one for TCRβ, revealed both common selection effects and notable differences between the two chains ( Supplementary Figures 4, 7–12 ). For this analysis, a subsample of CDR3 singletons from PBMC (23) was utilized alongside a thymocytes CDR3α repertoire (10) and model-generated data.

Among the observed differences, the following points are noteworthy: (i) there were no clear trends in the selection preferences toward individual amino acids ( Supplementary Figures 7A, 8A ) (ii) there was no consistent direction of Kidera factor changes for TCRα VJ pairs ( Supplementary Figures 7C, 8C ); (iii) there was an unchanged post-selection repertoire hydrophobicity ( Supplementary Figures 7F, 8F ); and (iv) there was a more complex structure of clusters enriched and depleted after selection, with DS and NY 2-mers prevalent in thymocyte-enriched clusters but not in naive data ( Supplementary Figure 10 ) and SS 2-mers being more abundant in clusters enriched in both model and DP thymocytes ( Supplementary Figures 9, 10 ). Notably, the common selection effects observed for both α and β chains were less pronounced in the case of TCRα ( Supplementary Figures 7, 8 ).

Additionally, we stored TCRs and patterns for the largest clusters enriched and depleted after thymic selection, obtained above in a database available at https://github.com/LuppovDaniil/thymic_selection_motifs_database.

To trace the signs of TCR driven lineage commitment during thymic selection, we utilize single-cell sequencing PBMC data derived from 88 healthy donors, totaling 178307 T-cells of different lineages (28). We aimed to demonstrate that enriched and depleted clusters are linked with particular T-cell phenotypes. We intersected CDR3s from these clusters ( Supplementary Figures 2, 3, 9 , 10 ) with a single-cell TCR repertoire and tested the linkage of these clusters with each particular cell phenotype.

Generally, the expected results were obtained. TCRs from clusters depleted after thymic selection were 10- to 20-fold less abundant in single-cell data than those from enriched clusters ( Figures 2C, D, G, H ), supporting our ability to identify TCRs disfavored by selection. Clusters enriched post-selection were mostly linked with naive phenotypes ( Figures 2A, B, E, F ). SP thymocyte clusters were mostly represented by CD8+ naive cells, which is expected since SP thymocytes in our research are of CD8+ lineage ( Figure 2B ).

We also found that the TCRβ of CD8+ cells are more prone to neighborhood enrichment than TCRβs of CD4+ cells ( Figure 2A ). The cluster enriched in CD8+ Recent Thymic Emigrants (RTE) subtype stood out against the other clusters enriched in CD8+ Naive cells ( Figure 2A ). Such an observation could be the sign of peripheral selection effects on the repertoire. Cluster 4 of SP thymocytes is also of particular interest, since it demonstrated abundance in both CD8+ RTE and Tregs ( Figure 2B ).

Moreover, we observed the abundance of TCRα from enriched SP thymocyte clusters in MAIT cells subtype ( Figure 2F cluster 3), which is again expected since MAIT cells are characterized by their “semi-invariant” TCRα (43).

Additionally, we assessed our assumption considering mostly naive origins of singletons in PBMC data. We inferred the phenotype for singletons from TCRβ and TCRα datasets using single-cell data and found that TCRβ singletons were indeed enriched in the naive subset (p < 10^-4, Fisher exact test); however, we were unable to detect the same enrichment for TCRα. Despite this peculiar inconsistency, we still believe that our assumption is correct, since the share of memory cells in the single-cell data was 10 times lower than the share of naive cells and the majority of T-cells were once naive.

To consolidate our findings, we analyzed SP thymocyte clusters utilizing standard 10X Genomics datasets containing both single-cell gene expression and TCR sequencing data. Generally, similar results were obtained: (i) cells carrying TCRβ from CD8+ thymocyte-enriched clusters were mostly of CD8+ phenotype; and (ii) TCRα-enriched clusters were mostly associated with MAIT cells ( Supplementary Note 1 ).

These results demonstrate TCR driven lineage commitment during thymic selection both for a well-characterized case of MAIT cells as well as a more complex one of CD8+ cells.

Next, we compared CD4+ and CD8+ naive cell TCR repertoires taken from the Qi et al. study (25). The most notable result of this comparison is from clusters enriched in CD4+ or CD8+ CDR3βs. We employed the TCRmodel2 to analyze the structure of the most prevalent CDR3 in the largest clusters within each group (33) and discovered a difference in convexity between enriched CDR3β CD8+ and CDR3β CD4+ structures ( Figures 3A, B ). It appeared that CD4+ enriched clusters were structurally flat, whereas CD8+ clusters were more convex. The centers of mass for the contacting region of CD8+ CDR3β (excluding the first four and last five residues) were located further from the loop center compared to the center of mass of CD4+ CDR3β ( Figure 3B ). However, in our sample of 18 structures, this difference was not statistically significant.

This fact may be explained by the conformation of peptides in MHC class I and MHC class II grooves. MHC class I, which is recognized by CD8+ cells, tends to present peptides with a middle bulge, whereas MHC class II, recognized by CD4+ cells, tends to present flat peptides (44).

We further investigated the structural difference between CD4+ and CD8+ CDR3βs using data stored at VDJdb (35). We took Human TCRs with available TCR:pMHC complexes harboring contacts between CDR3β and peptide. In total 154 structures were analyzed. We visualized all CDR3β dihedral angles in these structures via a canonical Ramachandran plot and discovered a distinct region on the plot that was abundantly inhabited by CD8+ CDR3β residues ( Figure 3C , green rectangle). The structural conformation corresponding to this region of the Ramachandran plot is mostly available for Glycines (45). The majority of structures with residues in this region were represented in it only by one residue. CD8+ TCRs harboring such a conformation were 2.83 times more frequent than CD4+ TCRs (p = 0.017, Fisher exact test).

This region at the right bottom of the Ramachandran plot in the Hollingsworth and Karplus work (45) was referred to as P_II’. P_II’ is viewed as a mirror region of P_II (more commonly referred to as Polyproline-II Helix), which is known to maximize polypeptide chain entropy and expose all hydrogen bond capable backbone atoms to the water. Additionally, this structural region frequently forms protein-binding motifs (46), which is relevant for TCR recognition of antigens. Thus, we would expect this CD8+ specific conformation to be relatively convex as we described it above. Moreover, the majority of residues in the discussed region were Glycines, which were in direct contact with the peptide. The data for the Ramachandran plot with corresponding PDB ID is available in Supplementary Table 5 .

Next, we considered the HLA allele influence on the thymic selection. It is still a subject of debate whether HLA alleles affect the selection or not (16).

We used CDR3s bulk sequencing data of three pairs of twins aged 20 to 23 to address this question (26). This data contains TCRβ repertoires for three pairs of twins (assigned as S, P and Q) in two replicas: day of vaccination before the shot (day 0) and a day prior to vaccination day (pre-day).

We compared their naive CDR3 clusters to the modeled background repertoire. On the one hand, we expected their enriched naive CDR3 clusters to be closely related with each other since twins from the same pair have the same HLA alleles. On the other hand, we anticipated twins from different pairs to be different from each other in context of their naive TCR repertoire.

We extracted the same number of CDR3s from each twin sample and then clustered the enriched sequences from each sample together. Having done this, we obtained a number of clusters composed of sequences from different samples. Supposedly, each cluster contains CDR3s, which are functionally similar to each other (recognize similar peptides in the MHC-peptide complex). Thus, we would anticipate CDR3s derived from the relative pair of twins to fall into the same clusters.

To verify this hypothesis, we analyzed the proportion of CDR3s from a particular twin in the clusters obtained. In this representation, each sample can be viewed as a discrete probability distribution of appearing in a particular cluster. This allowed us to calculate the Jensen–Shannon divergence between samples. We discovered that samples within each twin pair are much closer than samples from different twin pairs ( Figures 4A, B ).

Notably, twin pairs with the homozygous HLA-A 02:01 allele (pairs P and Q, please see Supplementary Table 2 ) appeared closer to each other than to the S pair, which had only one copy of HLA-A 02:01 allele.

The analysis of the clusters that were depleted during the selection showed similar results ( Figures 4C, D ). Thus, HLA alleles affect not only which TCRs will be promoted by selection but also the TCRs which will be eliminated by it, so called “holes” in the repertoire. Notably, we detected the effect of homo/heterozygosity of HLA-A gene alleles on the T-cell repertoire again but this time it was stronger—twin pairs, which were homozygous by HLA-A (P and Q), were 2–3 times closer to each other than to HLA-A heterozygous twin pairs (S) ( Figure 4D ).

We repeated the analysis in the same manner for two other sets of twins’ TCRs, one of which represented CD4+ naive cells, and obtained similar results ( Supplementary Note 2 ).

These results suggest that one can observe the imprint of HLA-based selection on a naive repertoire. Moreover, this imprint was observed in both enriched and depleted clusters, indicating that particular HLA alleles can be favorable for one TCR and unfavorable for another.