Introduction

The thymus is essential for T-cell development and TCR repertoire selection, but undergoes changes in output and function across the life of a mouse(Kondo et al., 2019, Montecino-Rodriguez and Dorshkind, 2023). In this study, we compared the αβTCR repertoire generated in the foetal thymus to that of the young adult, to investigate how these life-stages influence V(D)J gene usage and TCR repertoire diversity and distribution.

A diverse pool of αβTCR is essential for adaptive immunity and is generated during T-cell development in the thymus by the RAG-dependent joining of variable (V), joining (J), and diversity (D) gene segments at the β locus, and V and J gene segments at the α locus, with additional diversity coming from non-template nucleotides (junctional diversity)(Schatz and Swanson, 2011).

The thymus and the parathyroid organs, emerge from the third pharyngeal pouch and cleft during mid gestation of the foetal mouse(Gordon and Manley, 2011). Haematopoietic progenitors migrate from the foetal liver through the blood stream to colonise the thymus from embryonic day (E) 10-12. The first CD4-CD8-double negative (DN) thymocytes are already present at E11 and proliferate over the following days, producing a larger population by E14. CD4+CD8+ double positive (DP) and CD8+ immature single positive (ISP) thymocytes are first detectable at ∼E16, while mature CD4+CD8-CD3+ single positive (SP4) and CD4-CD8+CD3+ single positive (SP8) thymocytes can be detected at day E18(Solanki et al., 2018, Solanki et al., 2020). After birth, the thymus continues to grow rapidly and reaches its peak size at 4 to 6 weeks of age(Xiao et al., 2003).

During this developmental progression, the T cell receptor (TCR)β and TCRα loci are sequentially rearranged to produce T-cells with a diverse TCR repertoire. Successful rearrangement of TCRβ and pre-TCR expression are required for commitment to the αβT-cell lineage and differentiation from the CD25+DN (DN3) stage to DP cell, whereas cells that successfully rearrange the γδTCR commit to the γδT-cell lineage at the CD25+DN stage. TCRα rearrangement and signalling through the αβTCR are essential for MHC-restriction and differentiation from DP to SP cell (positive selection) and then also for deletion of T-cells bearing self-reactive TCR (negative selection). Pre-TCR signalling is required for survival, expansion and differentiation of cells that have completed TCRβ rearrangement, to produce a pool of cells which express a single functional TCRβ chain, in which TCRα rearrangement and αβTCR repertoire selection will occur(Dutta et al., 2021). In general, allelic exclusion of the TCRβ chain locus is believed to prevent developing thymocytes from simultaneously rearranging and expressing two TCRβ chains, whereas in DP cells the TCRα chain locus undergoes biallelic Vα-Jα rearrangement, with many possible rounds of rearrangement on each allele(Carico and Krangel, 2015, Carico et al., 2017, Genolet et al., 2012). Although the sequence of this developmental programme is the same between foetal and adult αβT-cell development, progression beyond the DN3 stage is less tightly linked to TCRβ expression in the foetal thymus(Hager-Theodorides et al., 2007).

There are several other fundamental differences between adult and foetal thymus that may affect the generation and selection of the TCR repertoire. In contrast to adult αβT-cells which are produced in a thymus structure that is mature, during foetal T-cell development, the thymus is still forming: the structure of the medulla and differentiation of the medullary thymic epithelial cell (mTEC) population, which express tissue restricted antigens for tolerance induction, are dependent on the haematopoietic compartment, with different requirements in foetal and adult thymus(Desanti et al., 2012). There are also differences in the progenitor cells that seed the thymus at different life-stages(Ramond et al., 2014), which arise from either the foetal liver or postnatal bone-marrow, and enter the thymus through different locations(Montecino-Rodriguez and Dorshkind, 2023). Some foetal haematopoietic stem cells (HSC) have been shown to have increased capacity to proliferate and to differentiate preferentially into innate-like lymphocytes compared with adult cells(Beaudin et al., 2016). Additionally, CD8 T-cells derived from HSCs from the foetal liver are functionally distinct from adult CD8 T-cells(Wang et al., 2016).

The foetal TCR repertoire is known to be less diverse than the adult repertoire. TdT is absent from the foetal thymus, and is first expressed a few days after birth, so foetal repertoires contain less non template nucleotide additions(Bogue et al., 1992). This has been confirmed by next generation TCR sequencing from mouse and human foetal and adult repertoires, which showed increased diversity after birth that was attributed to the postnatal increase in non-template nucleotide additions(Sethna et al., 2017, Britanova et al., 2016, Pogorelyy et al., 2017).

In this study we compared the TCR repertoire that is generated and selected in the foetal thymus with that from the young adult thymus. In addition to the reduction in non-template nucleotide additions seen in the embryonic repertoire, we identified biases in VJ usage and VxJ pairing of both TCRβ and TCRα repertoires. The foetal TCRβ repertoires were less diverse and less equally distributed with more clonal expansions than the young adult repertoires.

Results

Foetal thymocyte populations show clonotypic TCRβ expansion with a less diverse TCRβ repertoire than adult

To compare the TCRβ chain repertoire generated in the embryonic and young adult thymus we FACS-sorted DP (CD4+CD8+), SP4 (CD4+CD8-CD3+) and SP8 (CD4-CD8+CD3+) populations from E18.5 (foetal) and 4 week-old (young adult) thymus and sequenced their rearranged TCR transcripts. After RNA extraction from purified populations we used a TCR sequencing protocol and analysis pipeline which includes single-strand DNA ligation that tags each molecule of TCR α-and β-chain mRNA with a unique molecular identifier (UMI), allowing PCR bias to be corrected for in later analysis(Oakes et al., 2017, Uddin et al., 2019, Thomas et al., 2013). Firstly, to visualise the frequency distribution of the TCRβ repertoire in each population, we plotted the TCRβ abundance (number of copies of each sequence) against their proportion of the repertoire (Fig.1A). For each of the developmentally-defined thymocyte populations, the foetal thymus showed an increase in the proportion of abundant clones compared to adult mice in DP, SP4 and SP8 populations (Fig.1A). The frequency distributions of TCRβ abundances were fitted to a discrete power law, where the power law exponent corresponds to the gradient on the log-log plots (Fig.1A). The power law exponent was significantly lower in each foetal population compared to its adult counterpart (Fig.1B), demonstrating a difference in distribution of abundance of TCRβ clones between foetal and adult, consistent with an increase in clonotypic expansion of abundant clones in the foetal TCRβ repertoire. To confirm this, we calculated the proportion of the TCRβ repertoire represented by the top 1% most abundant TCRβ clones in each population. In foetal DP and SP4 populations the proportion of the TCRβ repertoire represented by the top 1% most abundant TCRβ clones was more than twice that of the adult populations (Fig1C), whereas the mean abundance of TCRβ sequences detected in the top 1% most abundant (expanded) clones was significantly higher in the foetal SP4 population only, suggesting a difference in equality of distribution in the top 1% most abundant (expanded) clones in the DP and SP8 populations compared to SP4 (Fig.1D).

Figure 1.

We therefore calculated standard indices of diversity, richness and distribution: Shannon entropy and the Gini index. Shannon entropy is derived from information theory and represents the complexity of the system, taking into account both the number of unique clonotypes and the frequencies of each clonotype, to generate a measure of richness and diversity(Shannon and Weaver, 1949). For each thymocyte population, the Shannon entropy was higher in the young adult TCRβ repertoire than in the foetal, demonstrating a richer more diverse repertoire in the adult thymocyte populations (Fig1E-F). The Gini index is a measure of the equality (evenness) of the distribution, and is frequently used in economics to demonstrate the distribution of a population’s wealth, where a lower index represents a more equal distribution(Ceriani and Verme, 2012). The Gini index of the distribution of TCRβ clonotypes was significantly lower in each adult population compared to its foetal counterpart, confirming a more equally distributed TCRβ repertoire in the young adult thymus (Fig1G-H).

In parallel we analysed the repertoire diversity, richness and distribution of the TCRα repertoires (SFig1). The fitted power law exponent of the frequency distribution of TCRα abundances was significantly lower in the foetal DP population compared to adult DP (SFig1A-B). However, we detected no significant differences between foetal and adult populations in the proportion of the TCRα repertoire represented by the top 1% most abundant TCRα clones (Sfig1C); the abundance of TCRα sequences detected in the top 1% most frequent (expanded) clones (SFig1D); the Shannon entropy (SFig1E-F); and Gini index (SFigG-H).

Foetal and young adult TCRβ and TCRα repertoires use distinct TRAV, TRAJ, TRBV and TRBJ gene segments

We next tested if the foetal and young adult repertoires showed differences in usage of TCR V and J gene segments. We found many differences between foetus and adult in proportional V and J segment usage within unique TCRβ and TCRα sequences for DP (SFig2), SP4 (SFig3) and SP8 (SFig4) populations. Preferential gene segment usage according to life-stage appeared to correlate with the chromosomal location of the gene segment, particularly for TRAV and TRBV usage within the DP population, where the adult repertoire favoured the 5’ gene segments, and the foetus the 3’ gene segments. To visualise these differences, we plotted heatmaps of mean proportional gene segment usage for unique TCR sequences, clustering each life-stage and thymocyte population, but showing each gene segment in chromosomal order (Fig2A-D). The samples clustered by life-stage first before cell-type (E18.5 populations clustered together) for both TCRα and TCRβ, suggesting that VJ gene usage is altered in the foetus compared to young adult. Interestingly, both TRAV and TRAJ showed clear preference for chromosomal location, with foetal populations favouring 3’ TRAV segments and 5’ TRAJ segments, whereas adult populations showed enrichment of 5’ TRAV and 3’ TRAJ (Fig2A-B). To test if this bias in gene usage was present in the unselected TCRα repertoire, we separated the DP population into cells that have not yet entered positive selection (CD3^-/loDP/CD69^-DP) and cells that have entered TCR repertoire selection (CD3^+/hiDP/CD69⁺DP). Both foetal DP populations clustered together and showed similar bias towards 3’ TRAV and 5’ TRAJ, suggesting that the bias was independent of positive selection. In the case of the adult DP populations, however, the more immature DP population clustered separately from the other young adult populations and showed weaker enrichment for 5’TRAV and 3’TRAJ than the populations that had undergone or were undergoing positive selection (CD3^+/hiDP, SP4 and SP8). This suggests that in the adult thymus, positive selection changes proportional TRAV and TRAJ gene usage, whereas in the foetal thymus an inherent bias in VJα gene usage persists after positive selection.

Heatmaps for TRBV and TRBJ also broadly showed chromosomal location-bias for gene usage at the different life-stages, with foetal populations preferentially using 3’TRBV and 5’TRBJ, although the pattern appeared less pronounced (Fig2. C-D). Therefore, to evaluate further the differences in TCRβ gene segment usage between foetal and adult repertoires, we carried out Principal Component Analysis (PCA) of the β-chain variable and joining gene counts of the unique TCRβ chains for each population (Fig. 2E-G). PCA separated the repertoires for DP, SP4 and SP8 by life-stage on PC1, confirming differences in TRBV and TRBJ usage between life-stage for each population.

Figure 2.

Figure 3.

Foetal and adult TCRβ and TCRα repertoires use distinct combinations of V and J segments

Given this difference in TRBV and TRBJ gene usage between foetal and adult, we next tested if TCRβ foetal and young adult thymocyte repertoires use different combinations of V and J gene segments. We carried out PCA using the counts of each VxJ combination from unique TCRβ sequences in foetal and adult DP, SP4 and SP8 populations (Fig3A-B). PCA separated by life-stage on PC1, with all adult repertoires clustering tightly together (Fig3A). When we compared by thymocyte population, each adult population clustered together, with adult SP4 separating from adult SP8 on PC2, and DP cells scoring in between, suggesting that PC2 might correspond to MHC-restriction of the adult populations (Fig3B). In contrast, the foetal populations were more dispersed, and did not segregate fully by cell type on PC2 (Fig3B). When we examined which VxJ combinations contributed most to PC1, we found strong correlation to the proportional TRBV and TRBJ usage detected: TRBV13-1 and TRBJ2-7 both showed adult bias in the analysis of individual gene usage, and in combination also contributed strongly to the positive score on PC1; whereas TRBV20 and TRBJ1-3 both showed foetal bias in the analysis of individual gene usage and contributed strongly to the negative (foetal) score on PC1, reflecting chromosomal gene segment location (SFig2-4 and Fig2-D).

We next compared usage of all possible VxJ combinations between the adult and foetal thymocyte populations for both TCRβ and TCRα repertoires, identifying all combinations that showed a significant change in proportional usage in foetus compared to young adult (Fig4A-C). The DP foetal population showed many increases (shown in red) and decreases (shown in blue) in proportional usage of VxJ combinations for TCRβ repertoire compared to young adult, which after selection were reduced in the SP populations. In the case of the TCRα chain, the foetal DP repertoire favoured VxJ combinations from the 3’ location of TRAV and 5’ location of TRAJ, but this bias was less pronounced in the SP populations (Fig. 4B-C). Concomitantly, the DP foetal population showed clear proportional decrease in usage of combinations of 5’ TRAV and 3’ TRAJ in comparison to young adult.

Figure 4.

Shorter CDR3 length and reduced non-template insertion length in foetal TCRβ repertoires

We next tested if the percentage of non-productive TCRβ rearrangements was different in foetal and young adult populations. The proportion of non-productive rearrangements was higher in the foetal SP8 population than adult (Fig5A). Furthermore, we examined non-template additions/deletions to V and J gene segments, which in the TCRβ chain include TRBD nucleotides. TRBD segments are difficult to identify as they are highly homologous and short. As previously reported(Sethna et al., 2017), the mean TCRβ non-template nucleotide insertion length was lower in the foetal DP and CD4SP populations compared to adult. To confirm this we then investigated the length of the complementary determining region 3 (CDR3) of the TCRβ chain as it is the most hypervariable region and primary site of antigen recognition. The mean predicted CDR3 length was shorter in all three populations in foetal TCRβ repertoires compared to young adult, confirming previous studies(Sethna et al., 2017) (Fig5B-C).

In contrast, when we compared predicted amino acid CDR3 sequences for the TCRα repertoires, we found no significant differences in the percentage of non-productive TCRα rearrangements between each adult and foetal population, or in the mean length of non-template insertions (SFig5A-B). The mean α-chain CDR3 length, however, was greater in the adult DP population than in foetal DP population (SFig5C).

Sharing of CDR3 sequences and CDR3 sequential amino acid triplets (k-mers) differs between adult and foetal TCRβ and TCRα repertoires

In order to investigate sharing between the predicted CDR3 repertoires within the same cell population from different embryos/mice of the same life-stage, we used the Jaccard Index of Similarity. This index measures the size of the intersection divided by the size of the union of finite sample sets(Jaccard, 1912). We calculated the Jaccard Index of Similarity for each CDR3 repertoire from each thymocyte population, comparing different mice/embryos to determine how similar each foetal repertoire is to the repertoire of the same thymocyte population from another foetus, and likewise for young adults. For each population, the Jaccard index was significantly higher in the foetal samples compared to adult, indicating that the CDR3 repertoires in each population from different embryos shared more sequences than the CDR3 repertoires from different individual adults (Fig5D-E). This increased sharing of CDR3 sequences in the foetal repertoires is consistent with shorter CDR3 length and shorter non-template insertions, leading to a repertoire that contains more genomic sequence than adult. We then compared the similarity in β-chain sequential CDR3 amino acid triplets (k-mers) for each population in the foetus and the adult. Conserved sequential amino acid triplets play an important role in determining antigen specificity and can aid in clustering together CDR3s that have similar specificities (Thomas et al., 2014, Glanville et al., 2017). Although the foetal repertoires from different individuals shared more CDR3 sequences with one another in all three populations than the adults, the adult populations shared more sequential motifs within their CDR3s with one another in DP and SP8 than the foetal populations (Fig5G-H). In support of this, when we carried out PCA, using the frequency distributions (counts) of each k-mer within CDR3 as input, the adult populations clustered together tightly scoring positive on PC1 and negative on PC2, whereas the foetal samples were more dispersed (Fig5I-J). Given that the foetal CDR3 repertoires contained fewer non-template nucleotides and showed more CDR3 sharing than adult, the increased adult CDR3 amino acid motif sharing is likely to be the result of selection of recurrent motifs by contact with pMHC during positive selection. In the case of the α-chain repertoires, we did not detect differences in life-stage sharing of TCRα predicted CDR3 sequences between individuals (SFig5D-E), but we found the adult DP, SP4 and SP8 repertoires shared more CDR3 amino acid triplets (k-mers) between one another than the foetal populations (SFig5F-G). However, in PCA, using the frequency distributions (counts) of each k-mer within the TCRα CDR3 as input, although most adult samples clustered together, the samples failed to segregate on either PC1 or PC2 by life-stage or by population (SFig5H-I).

Figure 5.

Different CDR1 and CDR2 usage between foetal and adult TCRβ and TCRα repertoires: less MHC-restriction in foetal repertoires

The CDR1 and CDR2 loops in the TCR are encoded in the V gene segments and typically contact the MHC’s conserved α-helices. They can determine MHC restriction during repertoire selection(Wong et al., 2019), so we examined CDR1xCDR2 combinations in the foetal and adult thymocyte populations. We carried out PCA, using the frequency distributions (counts) of each β-chain CDR1xCDR2 combination as input (Fig6A-C). The PCA clustered the populations by life-stage and cell-type, with the DP samples positioning between either SP population. To further visualise this, we generated a heatmap of the proportional usage of β-chain CDR1xCDR2 combinations for each foetal or adult thymocyte population (Fig. 6D). The foetal populations clustered together, showing preference for distinct CDR1xCDR2 combinations, which reflect the bias observed in their TRBV gene usage (Fig 2C). For both foetal and adult samples, DP populations were positioned between SP8 and SP4, suggesting a divergence in CDR1xCDR2 usage as a result of MHC restriction following positive selection. Indeed, referring back to the PCA (Fig6C), CDR1xCDR2 combinations that contributed strongly to the positive side on PC1, such as MSHET-SYDVDS (TRBV29) and SGHSN-HYEKVE (TRBV12-1), showed increased expression in SP8 and decreased expression in SP4 (consistent with MHCI-restriction), while combinations that contributed strongly to the negative side of PC1, such as NSQYPW-LRSPGD (TRBV31) and GKSSPN-SITVG (TRBV1), showed increased expression in SP4 and decreased expression in SP8 (consistent with MHCII-restriction).

Figure 6.

TCRα CDR1 and CDR2 may also be determined by MHC-restriction(Sim et al., 1996), so we additionally investigated TCRα CDR1xCDR2 combinations in the TCRα repertoires (Fig6E). Again, samples from each life-stage grouped together, and consistent with the strong life-stage dependent influence on TRAV use, several CDR1xCDR2 combinations showed an adult-bias in usage. PCA, using the frequency distributions (counts) of each predicted α-chain CDR1xCDR2 combination as input failed to segregate samples by life-stage on either PC1 or PC2 (SFig6A-B). However, in the case of the adult samples only, PC2 separated samples by population, with SP4 cells scoring negative on PC2 and SP8 cells scoring positive, suggesting that TCRα CDR1/CDR2 usage was influenced by MHC restriction in adult repertoire selection, consistent with other reports in mouse and human(Camaglia et al., 2023, Suo et al., 2023).

Young adult DP VJα gene segment usage becomes more foetal-like after hydrocortisone treatment

We observed distinct patterns of V and J gene usage between foetus and young adult, with foetal DP repertoires showing bias towards use of 3’ V and 5’ J rearrangements (Fig2-4).

During TCRβ gene rearrangement, these segments are closest to one another in the looping structure that enables V to DJ rearrangements (Skok et al., 2007), suggesting this might account for the foetal bias in proportional TRBVxTRBJ segment usage. In contrast to TCRβ locus rearrangement, the TCRα locus can undergo multiple rounds of rearrangement on each chromosome, sequentially moving 3’ to 5’ along the TRAV array, while moving 5’ to 3’ along the TRAJ array, with proximal pairs (3’ TRAV segments with 5’ TRAJ segments) initiating this sequence of rearrangements(Carico et al., 2017). We observed that the foetal DP repertoire showed TRAVxTRAJ bias by their chromosomal position. One possible explanation for this bias is that in the foetus progressive rounds of TCRα rearrangement are less common than in young adult, perhaps as a result of differences in the kinetics of T-cell development. Differentiation of foetal thymocytes occurs in a largely synchronized wave, whereas adult thymocyte production is unsynchronised, continuous and slower. Once the adult thymus has reached steady-state, homeostasis between thymocyte populations regulates the rate of differentiation to DP cell (Hager-Theodorides et al., 2007, Outram et al., 2009). To examine the influence of the rate of differentiation on VJ gene usage, we synchronized the differentiation of adult DP thymocytes by treating young adult mice with hydrocortisone (HC) to deplete the adult thymus of all but the most mature cells. At 2 days after treatment, the thymus was depleted of DP cells. Four days later (6 days after treatment), we FACS-sorted and TCR sequenced the replenished CD3^-/loDP (CD3^-/loCD4+CD8+), and CD3^+/hiDP (CD3^+/hiCD4+CD8+) populations. We plotted heatmaps of mean proportional V and J gene segment usage for unique TCR sequences in foetal, HC-treated and control young adult repertoires, clustering each life-stage and thymocyte population, but showing each gene segment in chromosomal order (Fig7A-D). Foetal and adult HC-treated TCRα repertoires clustered together, with untreated adult repertoires clustering separately (Fig7A-B). Interestingly, both foetal and HC-treated young adult populations favoured 3’ TRAV segments and 5’ TRAJ segments, in contrast to control young adult populations, that displayed bias for 5’ TRAV and 3’ TRAJ (Fig7A-B). Therefore, by synchronizing young adult thymus using HC, foetal and young adult gene segment usage became more alike, with 5’ TRAV and 3’ TRAJ bias. We next compared usage of all possible VxJ combinations between control and HC-treated young adult DP populations for all TCRα sequences, identifying combinations that showed a change in proportional usage in the HC-treated young adult compared to control young adult (Fig4A-C). The DP HC-treated young adult population favoured VxJ combinations from the 3’ location of TRAV and 5’ location of TRAJ along with clear proportional decrease in usage of combinations of 5’ TRAV and 3’ TRAJ in comparison to young adult control (Fig7C). In the case of the β-chain, analysis of proportional TRBV and TRBJ gene segment usage in unique β-chain rearrangements showed that samples clustered by life-stage first, with both untreated and HC-treated young adult repertoires clustering together (Fig7D-E). However comparison of proportional usage of all possible unique TCRβ VxJ combinations between HC-treated and control young adult DP populations showed a 5’ bias in TRBJ gene segment usage (Fig7F). These data indicate that when the adult young adult DP population arises quickly and synchronously its TRAV, TRAJ and TRBJ segment usage and TRAVxTRAJ combinations closely ressemble the foetal DP repertoire, suggesting that the kinetics of transition to DP cell are important in determining gene segment usage.

Figure 7.

Discussion

Here we employed a bulk population-based strategy to investigate the TCRβ and TCRα chain repertoires from developmentally defined thymocyte populations during T-cell development. Our study revealed many differences between the foetal and young adult thymus in TCR gene segment usage, repertoire diversity, and clonality by bulk sequencing separate TCRβ and TCRα repertoires. These data indicate that MHC restriction and positive selection have a weaker impact on the TCR repertoire in foetal thymus compared to adult. Overall, the differences between the foetal and adult thymus TCR repertoires are consistent with the foetal thymus producing αβT-cells with properties and functions that are distinct from adult T-cells, and more innate-like: their repertoire is less diverse, more closely encoded by genomic sequence, less governed by MHC-restriction, and with preference for particular gene segment usage. Several recent studies have demonstrated distinct and more innate-like features of foetal and neonatal T-cell populations in both humans and mice (Beaudin et al., 2016, Wang et al., 2016, Rackaityte and Halkias, 2020, Smith et al., 2018). The first wave of foetal αβT-cells that leave the thymus must provide early protection against infection in the neonatal animal but also need to be tolerant to both self and maternal MHC/antigens. It is possible that the distinct features and weaker MHC-restriction of the foetal TCR repertoire have evolved to provide these two important features of foetal adaptive immunity. It will be interesting to investigate if the foetal-bias in gene segment usage confers specific properties on foetal TCRs, with respect to their binding specificity, affinity and avidity.

In validation of this approach, we were able to quantify known differences in TCRβ repertoire diversity, CDR3 length and non-template insertion length between foetal and adult(Sethna et al., 2017). The foetal TCR repertoire contained fewer non-template nucleotides than adult, indicating that it was more closely aligned to genomic sequence. Consistent with this, the predicted foetal β-chain CDR3 repertoires shared more sequences with one another than their young adult counterparts. Our study also showed that the foetal TCRβ repertoire had a less equal distribution in DP, SP4 and SP8 populations compared to adult, and increased clonal expansion of the most abundant clones in DP and SP4 populations. Thus, the reduction in TCRβ repertoire diversity observed in the foetal thymus could not simply be attributed to the reduction in N-nucleotide insertions. These differences in TCRβ repertoire distribution and clonal expansion may be the result of differences in the kinetics of expansion and differentiation of thymocyte subsets in the foetus, leading to greater clonal expansion of TCRβ clones following β-selection. The transition from DP to SP cell in the foetus takes just 2-3 days at most, as DP cells first appear around E16, and the mature SP populations were FACS-sorted on E18.5, whereas in the adult thymus differentiation from DN3 to SP cell may be slower(Ross et al., 2014, Solanki et al., 2020).

We found many differences in proportional TRBV, TRBJ, TRAV, and TRAJ gene usage between foetal and adult thymocyte populations, with foetal populations showing proportional bias towards 3’ TRBV and 5’TRBJ usage. We observed a particularly striking bias towards use of 3’ TRAV and 5’ TRAJ gene segments in the foetus. In mammals, the gene segments encoding TCRδ and TCRα chains are found at a single genetic locus, with TCRδ rearrangement taking place at the DN stage of development (prior to TCRα rearrangement in the αβ-lineage committed DP population). The TCRα locus can undergo multiple rounds of rearrangement on each chromosome, with proximal pairs (3’ TRAV segments with 5’ TRAJ segments) initiating this sequence of rearrangements. However, as TCR8 gene segments are nested between the TRAV and TRAJ array and also amongst the TRAV array, TCR8 gene rearrangement splices out 3’ TRAV segments, facilitating 5’ TRAV segment pairing with TRAJ segments(Carico et al., 2017). It is therefore possible that the 3’ TRAV and 5’ TRAJ bias observed in the foetus may in part be the result of primary TCRα rearrangements taking place on chromosomes which have not already rearranged TCR8, so that 5’ TRAV segments are not favoured. In the future, it will be important to test experimentally if the combination of 3’ TRAV segments with 5’ TRAJ segments contribute to TCRs which facilitate functions associated with foetal T-cells, as this might account for the evolutionary conservation of the genomic organization of the TCRα and TCR8 loci across mammals.

Bias in TCRα gene segment usage by chromosomal location has been reported in human DN thymocyte populations by single cell RNAseq but was found to diminish as thymocytes mature and progressive recombination of the TCRα loci occurs(Park et al., 2020). In contrast, in the mouse foetal thymus, we found that the bias in proportional TRAV and TRAJ gene segment use by chromosomal position persisted in the positively selecting CD69⁺DP and mature SP4 populations, but was not apparent in the adult CD3^-/loDP population. Our data therefore suggest that progressive rounds of TCRα rearrangement are less common in the foetal DP population. Indeed, when we depleted the young adult thymus by HC-treatment, we found an increase in 3’TRAV and 5’TRAJ gene segment usage in the recently differentiated recovering DP population, making these repertoires more similar to the foetal DP TCRα repertoires. In support of this, the VJ gene usage of HC-treated young adult DP populations clustered together with foetal. Thus differences between foetal and control young adult TCRαV-J rearrangements in DP cells may be the result of slower differentiation in the adult thymus, allowing more time for multiple rounds of TCRαV-J rearrangements in the DP population. Consistent with this, the foetal DP population showed a different TCRα frequency distribution than the adult DP population, with a lower power law exponent.

Comparison of proportional use of TCRβVxJ and TCRαVxJ combinations between adult and foetal populations also highlighted the 3’ V gene segment preference in the foetal DP population, and the 5’ V gene segment preference in the adult DP population for both β-chain and α-chain. Interestingly, life-stage had a greater impact on TCRβ and TCRα gene segment usage than MHC-restriction in all populations.

PCA clustered β-chain VxJ combination counts by population in the adult samples, showing the impact of MHC restriction (positive selection) on the selected adult SP repertoires. PCA also clustered the adult counts of β-chain CDR1xCDR2 combinations by cell type, confirming impact of MHC restriction on Vβ gene usage. However, foetal populations failed to segregate fully by cell type in either PCA, suggesting that MHC restriction/positive selection has less influence on β-chain V and J gene segment use in the foetal SP populations. In support of this, young adult β-chain CDR3 repertoires shared more sequential amino acid motifs (k-mers) than foetal repertoires, suggesting a stronger influence of pMHC on adult repertoire selection.

In the case of the TCRα repertoire the DP population showed many increases in proportional use of 3’TRAV and 5’TRAJ combinations in the foetal thymus, whereas the adult DP population showed clear proportional bias to 5’TRAV to 3’TRAJ combinations. PCA clustered the adult α-chain CDR1xCDR2 combinations by cell type, confirming impact of MHC restriction on Vα gene usage in the SP populations. However, as also observed in the analysis of β-chain repertoires, the foetal samples did not cluster by cell-type, suggesting a weaker influence of MHC-restriction/positive selection on both the TCRβ and TCRα repertoires that are selected in the foetal thymus.

Our approach of bulk TCR sequencing from FACS-sorted thymocyte populations had the advantages of scale (allowing the sequencing of thousands of rearrangements from thousands of cells from each embryo or mouse), at relatively low cost. However, it had the disadvantage that we were unable to investigate which TCRβ rearrangements pair with which TCRα rearrangements in single cells. Clearly, in future it will be important and interesting to expand this analysis to include TCRβ/α pairing by single cell RNAseq.

Materials and methods

Mice

C57/BL6 mice were bred and maintained in specific pathogen-free conditions at University College London (UK) under UK Home Office regulations. For hydrocortisone (HC) treatment, mice were injected intraperitoneally with 0.6mg/gram of body weight pure HC sodium phosphate (Sigma Aldrich, cat. no: BP188) in sterile PBS, and analysed after 6 days.

Fluorescence activated cell sorting

Each foetal or young adult thymus was disaggregated into a single cell solution using the back of a syringe on a 70-M cell strainer. The cells were washed through the filter into falcon tubes using ice cold FACs buffer (1x AIM-V medium (research grade), AlbuMAX Supplement (Gibco, cat. no: 31035025)) to eliminate any large clumps. Cells were then counted using Accuri C6 flow cytometer and then pelleted by centrifugation for 5 minutes at 1400 rpm and supernatant was removed. The whole cell pellet from the entire organ was stained with a panel of directly conjugated antibodies supplied by Biolegend (San Diego, US) or eBioscience (San Diego, US). See Supplementary Table 1 and 2 for antibody panels. Firstly, an antibody mastermix of either 1:50 (foetal thymus) or 1:100 (adult thymus) in fresh ice cold FACS buffer was prepared. Subsequently, either 100 μL (foetal thymus), 400 μL (adult thymus) of antibody mastermix was then added to pelleted cells and the cells were then stained on ice for 30 minutes in the dark. The stained cells were then washed using 2 mL of FACS buffer and pelleted by centrifugation for 5 minutes at 1400 rpm. The cells were resuspended in either 100 μL (foetal thymus), or 500 μL to 1 mL (thymus) of fresh ice cold FACS buffer depending on the cell counts. To obtain CD3^-/loDP (CD4⁺CD8⁺CD3^-/lo), CD69^-DP (CD4⁺CD8⁺CD69^-), CD69⁺DP (CD4⁺CD8⁺CD69⁺), CD3^+/hi DP (CD4+CD8+CD3^+/hi), SP4 (CD4⁺CD8^-CD3⁺), SP8 (CD4^- CD8⁺CD3⁺) thymocyte cell suspensions were sorted using a BD FACS Aria III. CD69 expression was used to differentiate DP populations in the foetus, as anti-CD3 staining is dull in DP thymocyte populations at E18.5(Solanki et al., 2018). Only live cells were gated and collected using forward and side scatter.

The cells were collected in 100 μL (foetus) or 400 μL (young adult) of FACS buffer and were then pelleted by centrifugation for 25 minutes at 1400 rpm at 4 °C. The supernatant was removed and 100 μL of Extraction buffer from PicoPure RNA isolation kit (Applied Biosystems, cat. no. KIT0204) was added. Finally, the samples were incubated for 30 minutes at 42 °C on a shaking heatblock and then stored at -80 °C until the RNA was ready to be extracted.

RNA extraction

RNA was extracted by PicoPure RNA isolation kit (Applied Biosystems, cat. no. KIT0204) according to manufacturer’s instructions and eluted in 11 µL of elution buffer. RNA concentration was assessed using a spectrophotometer.

TCR amplification and sequencing

We used the protocol for TCR amplification and sequencing as described (Oakes et al., 2017, Uddin et al., 2019), with some minor modifications for mouse samples. See Supplementary Table 3 for the primer list with the sequences and Supplementary Table 4 for a full list of reagents.

RNA was first DNase treated. A maximum of 500 ng of RNA was mixed with 1 μL RQ1 10x Buffer (Promega, cat. no: M6101) and 1 μL RQ1 DNase (1 U/μL, Promega, cat. no: M6101) and then incubated at 37 °C for 30 minutes. 1 μL of RQ1 DNase Stop Solution (Promega, cat. no: M6101) was then added and the mix was incubated at 65 °C to stop the reaction. The RNA was then reverse transcribed by adding 1.5 μL of 10 μM TRAC2 primer, 1.5 μL of 10 μM TRBC3 primer, and 1.5 μL of 10 mM dNTPs (Promega, cat. no: U1515) and incubating at 65 °C for 5 minutes. The reaction was then rapid cooled down on ice and 6 μL of 5x FS Buffer (Invitrogen, cat. no: 18080044), 1.5 μL of 0.1 M DTT (Invitrogen, cat. no: 18080044), 1.5 μL of RNasin Ribonuclease Inhibitor (40 U/μL, Promega, cat. no: N2111) and 1.5 μL of SuperScript III RT (200 U/μL, Invitrogen, cat. no: 18080044) was added and it was incubated at 55 °C for 30 minutes, followed by 70 °C for 15 minutes. These reactions were purified by Minelute PCR purification (Qiagen, cat. no: 28004) following the manufacturer’s instructions and eluted in 10.5 μL nuclease-free water.

The cDNA was then ligated to the M13 ligation primer which is composed of the Illumina SP2 primer, an 8 base spacer and a 12 base Unique Molecular Identifier (UMI). This UMI consists of two random hexamers separated by the 8 base spacer. On the 5’ end the primer is phosphorylated, and on the 3’ end it is blocked with a Spacer C3 moiety to prevent primer concatemerization. The ligation mix consisted of: 10 μL cDNA, 3 μL of bovine serum albumin (BSA)/hexamine cobalt chloride (HCC) mixture (1 mg/mL BSA, 10 mM HCC), 3 μL of T4 RNA ligase buffer (NEB, cat. no: M0204S), 1 μL of 10 mM ATP (NEB, cat. no: M0204S), 1 μL of 10 μM M13 Ligation primer, 2 μL of T4 RNA Ligase (10 000 U/mL, NEB, M0204S), and 10 μL of PEG 8000 50% (NEB, cat. no: M0204S). The ligation mix was then incubated for at least 18 hours up to a maximum of 23 hours at 16 °C and then heat inactivated at 65 °C for 10 minutes. 70 μL of nuclease-free water was added to the ligation products. This was then purified by adding 50 μL of AMPure XP Beads (Beckman Coulter, cat. no: A63881) to the 100 μL ligation product, with some modifications to manufacturer’s instructions: 300 μL of 80% ethanol was used for the washing steps and the product was eluted in 31 μL nuclease-free water.

The purified ligation product was then amplified in a PCR reaction consisting of 31 μL purified product, 10 μL of 5x HF Buffer (NEB, cat. no: M0530L), 2.5 μL of 10 μM m-alpha-RC1 primer, 2.5 μL of 10 μM m-beta-RC1 primer, 2.5 μL of 10 μM SP2-M13 primer, 1 μL of 10 mM dNTPs (Promega, cat. no: U1515), and 0.5 μL of Phusion Polymerase (NEB, cat. no:M0530L). The m-alpha-RC1 and m-beta-RC1 primers hybridize to the constant region of the α-chain and β – chain genes in the mouse respectively. The SP2-M13 primer hybridizes to the sequence introduced by the M13 ligation primer earlier during the ligation step. PCR-1 was run using the following conditions: initial denaturation at 98°C for 3 minutes, followed by 4 cycles of denaturation at 98°C for 15 seconds, annealing at 69°C for 30 seconds and extension at 72°C for 40 seconds, followed by a final extension at 72°C for 5 minutes.

The 50 μL PCR product was then purified by adding 40 μL of AMPure XP Beads (Beckman Coulter, cat. no: A63881) with some modifications to manufacturer’s instructions: 200 μL of 80% ethanol was used for the washing steps. At this stage the samples were split into two: TCRα and TCRΠ, each eluted in 31 μL nuclease-free water.

The purified PCR product was then amplified in a PCR reaction to add indices and Illumina adaptors P5 and P7 and sequencing primer 1 (SP1) sequence. The PCR reaction mix consisted of: 31 μL purified product, 2.5 μL of 1 μM SP1 index (mSP1-6N-I-X-αRC1 or mSP1-6N-I-X-βRC1, for TCRα and TCRβ respectively), 2.5 μL of 10 μM SP2 index (P7-LX), 2.5 μL of 10 μM SP1-P5 primer, 10 μL of 5x HF Buffer (NEB, cat. no: M0530L), 1 μL of 10 mM dNTPs (Promega, cat. no: U1515), and 0.5 μL of Phusion Polymerase (NEB, cat. no: M0530L). PCR-2 was run using the following conditions: initial denaturation at 98°C for 3 minutes, followed by 6 cycles of denaturation at 98°C for 15 seconds, annealing at 69°C for 30 seconds and extension at 72°C for 40 seconds, followed by a final extension at 72°C for 5 minutes. The 50 μL PCR product was then purified by adding 40 μL of AMPure XP Beads (Beckman Coulter, cat. no: A63881) with some modifications to manufacturer’s instructions: 200 μL of 80% ethanol was used for the washing steps, and the product was eluted in 30 μL nuclease-free water.

The purified PCR product was then amplified in a qPCR reaction, so that the reaction can be monitored and stopped in real time when the reaction passes the threshold of 0.1 ΔRn to prevent overamplification and minimise the introduction of PCR artifacts. SYBR Green I Nucleic Acid Gel Stain, 10,000x concentrate in DMSO (ThermoFisher Scientific, cat. no: S7567) was diluted firstly 1:100 in DMSO and this stock solution was stored at -20 °C. Prior to the assembly of the PCR mix, this stock solution was diluted a further 1:50 in nuclease-free water. The PCR reaction consisted of: 28 μL purified product, 10 μL of 5x HF Buffer (NEB, cat. no: M0530L), 5 μL of SYBR Green (1:50), 1.25 μL of 10 mM dNTPs (Promega, cat. no: U1515), 1 μL of Rox (Invitrogen, cat. no: 12223012), 2.5 μL of 10 μM P5 primer, 2.5 μL of 10 μM P7 primer and 0.5 μL of Phusion Polymerase (NEB, cat. no: M0530L). The qPCR was run using the following conditions: initial denaturation at 98°C for 3 minutes, followed by a variable number of cycles of denaturation at 98°C for 15 seconds, annealing at 69°C for 30 seconds and extension at 72°C for 40 seconds, followed by a final extension at 72°C for 5 minutes.

The reaction was stopped when the samples passed the threshold of 0.1 ΔRn. The 50 μL PCR product was then purified by adding 40 μL of AMPure XP Beads (Beckman Coulter, cat. no: A63881) with some modifications to manufacturer’s instructions: 80% ethanol was used for the washing steps and the product was eluted in 30 μL nuclease-free water.

The purified products were then quantified using Qubit dsDNA high sensitivity reagents (ThermoFisher Scientific, cat. no: Q32854) and the size was confirmed with the TapeStation System HSD1000 (Agilent, cat. no: 5067-5584, cat. no: 5067-5585). The expected size is around 650 bp for a successful library preparation.

Sequencing

The concentrations acquired by the Qubit and the size determined by the TapeStation, were used to calculate the molarity of the samples in nM. The samples were then diluted to 12 nM and pooled into one 12 nM library containing samples with different indices. 84 samples were prepared per sequencing run (42 α-chain and 42 β -chain libraries). A Vacuum Concentrator Centrifuge was used to pellet the library and remove excess volume, so that 30 μL of the pooled library could then be run on the Pippin Prep System using a Pippin Gel Cassette 1.5% w/v agarose dye free 250 bp-1.5 kb DNA size range (Sage Science, cat. no: CDF1510) to size select for sequences between 350 to 750 bp, following the manufacturer’s instructions.

40 μL of library eluted from the Pippin Prep was then quantified using a Qubit and TapeStation System HSD1000 as detailed above. The expected molarity was between 2 nM and 5 nM. The library was then denatured and diluted to 1.2 pM according to standard Illumina’s protocols (for 4 nM libraries) and was spiked with 22% of 20 pM PhiX control (Illumina, cat. no: FC-110-3001). The 1.2 pM library was sequenced on the NextSeq using the NextSeq 500/550 Mid Output Kit v2.5 (300 Cycles) (Illumina, cat. no: 20024905) at UCL genomics.

Error correction and outputting

The NextSeq outputs files in the format named binary based call (.bcl) which were converted into FASTQ files using bcl2fastq for downstream processing using a pipeline of scripts described previously(Oakes et al., 2017, Thomas et al., 2013): Decombinator_v3.1 (available at: https://github.com/innate2adaptive/Decombinator/) in Python 2.7. This pipeline identifies the errors introduced by PCR amplification using the UMI attached during the ligation step and thus produces an output of a corrected abundance of the TCR in the sample along with determining V, J and CDR3 regions for each TCR.

Downstream analysis

For further downstream analyses, R Studio was used. The tidyverse set of packages were used for data manipulation and visualisations, in particular ggplot2(Wickham, 2016, Wickham H, 2019). Dotplots show mean±c.i (package ggpubr(Kassambara, 2023a)) and statistical comparisons were carried out by unpaired Student’s t-test or Welch’s t-test as appropriate using the package rstatix(Kassambara, 2023b). Pheatmap and ComplexHeatmap were used to generate Pearson correlation heatmaps(Gu et al., 2016, Kolde, 2019).

TCR abundance distribution

TCR repertoire abundance distributions are typically L-shaped, with a high number of distinct TCRs or clonotypes present only once, and a low number of hyperexpanded clones. Therefore, the TCR abundance (number of copies of each sequence) was plotted against their proportion of the repertoire in a log-log plot. This transformed distribution follows approximately a linear distribution apart from the largest clones and can be fitted to a discrete power law (f(x) = kx^-α) as seen in previous studies(Oakes et al., 2017, Joshi et al., 2019). The TCR abundance frequency was fitted to a discrete power law using maximum likelihood estimation(Bauke, 2007, Clauset et al., 2009) using the PoweRLaw package(Gillespie, 2015). The power law exponents (α) were then plotted and compared.

Diversity indices

Before calculating diversity indices, α and β identified TCRs or CDR3s were subsampled (rarefied) to a number lower than the smallest repertoire using the package vegan(Oksanen et al., 2022) in order to correct for differences in diversity due to sample size. The mean diversity indices (Shannon Entropy, Gini Index and Jaccard Index of Similarity) were then calculated from 1000 repeats of this random sampling. The Shannon Entropy was computed using the package vegan(Oksanen et al., 2022). Gini index was computed using ineq package(Zeileis, 2014). The base R function dist was used to calculate the Jaccard Index of Similarity.

Principle Component Analysis

Before PCA was applied, raw counts were first log₁₀ transformed to account for differences in sample size, using a pseudocount of 0.01 for any counts that were not observed. Z-scores were then calculated from these values by subtracting the mean and dividing by the standard deviation. PCA was then computed using the base R function prcomp and factoextra package was used to investigate the results(Kassambara A, 2020).

VxJ combinations

Proportional VxJ gene usage of total TCRs was calculated for each sample and only VxJ combinations that were detected in all samples were compared. Statistical comparisons were carried out by unpaired Student’s t-test followed by FDR-adjustment (5%, Benjamini-Hochberg procedure) of p values.

Sequential amino acid triplet analysis

Conserved sequential amino acid triplets (k-mers) can be important in determining antigen specificity and can help cluster CDR3 sequences according to specificity (Thomas et al., 2014, Glanville et al., 2017) and were extracted from the predicted amino acid CDR3 sequences using the immunarch package(Team, 2019).

Authors’ Contributions

JR, BC and TC conceived the study; JR and TC designed experiments and wrote the manuscript. JR analysed sequencing data; JR, C-IL, DY and SR performed experiments. C-IL and BC critically reviewed the manuscript.

Acknowledgements

This work was funded by the MRC (MR/P000843/1); MR/S037764/1) and BBSRC (BB/T020970/1). JR was supported by a studentship from the BBSRC London Interdisciplinary Biosciences Consortium (1903458). Research at UCL Great Ormond Street Institute of Child Health is supported by the NIHR Biomedical Research Centre at Great Ormond Street Hospital and UCL. The authors declare no competing financial interests.

Data Availability

TCR sequencing data will be made available on UCL Research Data Repository.

Supplementary Figure Legends

Supplementary Figure 1:

Supplementary Figure 2:

Supplementary Figure 3:

Supplementary Figure 4:

Supplementary Figure 5:

Supplementary Figure 6:

Supplementary Tables

Supplementary Table 1

Supplementary Table 2

Supplementary Table 3

Supplementary Table 4