A cell atlas of human thymic development defines T cell repertoire formation

Abstract

The thymus provides a nurturing environment for the differentiation and selection of T cells, a process orchestrated by their interaction with multiple thymic cell types. We used single-cell RNA sequencing to create a cell census of the human thymus across the life span and to reconstruct T cell differentiation trajectories and T cell receptor (TCR) recombination kinetics. Using this approach, we identified and located in situ CD8αα⁺ T cell populations, thymic fibroblast subtypes, and activated dendritic cell states. In addition, we reveal a bias in TCR recombination and selection, which is attributed to genomic position and the kinetics of lineage commitment. Taken together, our data provide a comprehensive atlas of the human thymus across the life span with new insights into human T cell development.

Constructing the human thymus cell atlas.

We analyzed human thymic cells across development and postnatal life using scRNA-seq and spatial methods to delineate the diversity of thymic-derived T cells and the localization of cells constituting the thymus microenvironment. With T cell development trajectory reconstituted at singlecell resolution combined with TCR sequence, we investigated the bias in the VDJ recombination and selection of human TCR repertoires. Finally, we provide a systematic comparison between human and mouse thymic cell atlases.

Fig. 1. Cellular composition of the developing human thymus.

(A) Schematic of single-cell transcriptome profiling of the developing human thymus. (B) Summary of gestational stage/age of samples, organs (circles denote thymus; rectangles denote fetal liver or adult bone marrow, adult spleen, and lymph nodes), and 10x Genomics chemistry (colors). (C) UMAP visualization of the cellular composition of the human thymus colored by cell type (DN, double-negative T cells; DP, double-positive T cells; ETP, early thymic progenitor; aDC, activated dendritic cells; pDC, plasmacytoid dendritic cells; Mono, monocyte; Mac, macrophage; Mgk, megakaryocyte; Endo, endothelial cells; VSMC, vascular smooth muscle cells; Fb, fibroblasts; Ery, erythrocytes). (D) Same UMAP plot colored by age groups, indicated by post-conception weeks (PCW) or postnatal years (y). (E) Dot plot for expression of marker genes in thymic stromal cell types. Here and in later figures, color represents maximum-normalized mean expression of marker genes in each cell group, and size indicates the proportion of cells expressing marker genes. (F) RNA smFISH in human fetal thymus slides with probes targeting stromal cell populations. Top left: Fb2 population marker FBN1 and general fibroblast markers PDGFRA and CDH5. Top right: Fb1 marker GDF10, FBN1, and CDH5. Middle left: Fb1 marker COLEC11 and FBN1. Middle right: Fb1 marker ALDH1A2, VSMC marker ACTA2, and FBN1. Bottom left: TEC(myo) marker MYOD1. Bottom right: Epithelial cell marker EPCAM and TEC(neuro) marker NEUROG1. Data are representative of two experiments. (G) Relative proportion of cell types throughout different age groups. Dot size is proportional to absolute cell numbers detected in the dataset. Statistical testing for population dynamics was performed by t tests using proportions between stage groups. The x axis shows age of samples, which are colored in the same scheme as (D).

Fig. 2. Thymic seeding of early thymic progenitors (ETPs) and T cell differentiation trajectory.

(A) UMAP visualization of ETP and fetal liver hematopoietic stem cells (HSCs) and early progenitors. NMP, neutrophil-myeloid progenitor; MEMP, megakaryocyte/erythrocyte/mast cell progenitor. (B) The same UMAP colored by organ (liver in blue, thymus in yellow/red). (C) UMAP visualization of developing thymocytes after batch correction. DN, doublenegative T cells; DP, double-positive T cells; SP, single-positive T cells; P, proliferating; Q, quiescent). The data contain cells from all sampled developmental stages. Cells from abundant clusters are downsampled for better visualization. The reproducibility of structure is confirmed across individual samples. Unconventional T cells are in gray. (D to F) The same UMAP plot showing CD4, CD8A, and CD8B gene expression (D), CDK1 cell cycle and RAG1 recombination gene expression (E), and TCRa, productive TCRb, and nonproductive TCRβ VDJ genes (F). (G) Heat map showing differentially expressed genes across T cell differentiation pseudotime. Top: The x axis represents pseudo-temporal ordering. Gene expression levels across the pseudotime axis are maximum-normalized and smoothed. Genes are grouped by their functional categories and expression patterns. Bottom: Cell type annotation of cells aligned along the pseudotime axis. Colors are as in (C). (H) Scatterplot showing the rate of productive chain detection within cells in specific cell types (x axis) and the ratio of nonproductive/productive TCR chains detected in specific cell types (y axis). Left: TCRβ; right, TCRa. (I) Graph showing correlation-based network of transcription factors expressed by thymocytes. Nodes represent transcription factors; edge widths are proportional to the correlation coefficient between two transcription factors. Transcription factors with significant association to specific cell types are depicted in color. Node size is proportional to the significance of association to specific cell types.

Fig. 3. Identification of GNG4⁺ CD8αα T cells in the thymic medulla.

(A) UMAP visualization of mature T cell populations in the thymus. Axes and coordinates are as in Fig. 2C. (The cell annotation color scheme used here is maintained throughout this figure.) (B) Dot plot showing marker gene expression for the mature T cell types. Genes are stratified according to associated cell types or functional relationship. (C) Scatterplot showing the ratio of nonproductive/productive TCR chains detected in specific cell types in TCRa chain (x axis) and TCRβ chain (y axis). The gray arrow indicates a trendline for decreasing nonproductive TCR chain ratio in unconventional versus conventional T cells. (D) Scatterplot showing the relative abundance of each cell type between fetal liver and thymus (x axis) and before and after thymic maturation (delimited at 10 PCW) (y axis). Gray arrow indicates trendline for increasing thymic dependency. (E to H) Scatterplots comparing the characteristics of unconventional T cells based on CD8A versus CD8B expression levels (E), KLRB1 versus ZBTB16 expression levels (F), TCRa productive chain versus TRDC detection ratio (G), and TRDV1 versus TRDV2 expression levels (H). Gray arrows or lines are used to set boundaries between groups [(E), (G), (H)] or to indicate the trend of innate marker gene expression (F). (I) RNA smFISH showing GNG4, TNFRSF9, and CD8A in a 15 PCW thymus. Lower right panel shows detected spots from the image on top of the tissue structure based on 4’,6-diamidino-2-phenylindole (DAPI) signal. Color scheme for spots is the same as in the image. (J) FACS gating strategy to isolate CD8αα(I) cells (live/CD3⁺/CD4^/CD137⁺) and Smart-seq2 validation of FACS-isolated cells projected to the UMAP presentation of total mature T cells from the discovery dataset (lower left). GNG4 expression pattern is overlaid onto the same UMAP plot (lower right).

Fig. 4. Recruitment and activation of dendritic cells for thymocyte selection.

(A and B) UMAP visualization of thymic DC populations (A) and dot plot of their marker genes (B). (C) Heat map of chemokine interactions among T cells, DCs, and TECs, where the chemokine is expressed by the outside cell type and the cognate receptor by the inside cell type. (D) Schematic model summarizing the interactions of TECs, DCs, and T cells. The ligand is secreted by the cell at the beginning of an arrow, and the receptor is expressed by the cell at the end of that arrow. (E) Left: RNA smFISH detection of GNG4, XCR1, and FOXP3 in 15 PCW thymus. Right: Computationally detected spots are shown as solid circles over the tissue structure based on DAPI signal. Color schemes for circles are the same as in the image. (F to H) Sequential slide sections from the same sample are stained for the detection of LAMP3, AIRE, and XCR1 (F), LAMP3, ITGAX, and CD80 (G), and LAMP3 and FOXP3 (H). Spot detection and representation are as in (E). Data are representative of two experiments.

Fig. 5. Intrinsic bias in human TCR repertoire formation and selection.

(A) Heat map showing the proportion of each TCRβ V, D, and J gene segment present at progressive stages of T cell development. Gene segments are positioned according to genomic location. (B) Same scheme as in (A) applied to TCRα V and J gene segments. Although there is a usage bias of segments at the beginning of development, segments are evenly used by the late developmental stages, indicating progressive recombination leading to even usage of segments. (C and D) Schematics illustrating a hypothetical chromatin loop that may explain genomic location bias in recombination of TCRβ locus (C) and the mechanism of progressive recombination of TCRα locus leading to even usage of segments (D). (E) Principal components analysis plots showing TRBV or TRAV and TRAJ gene usage pattern in different T cell types. Arrows depict T cell developmental order. For TRBV, there is a strong effect from beta selection, after which point the CD4⁺ and CD8⁺ repertoires diverge. The development for TRAV+ TRAJ is more progressive, with stepwise divergence into the CD4⁺ and CD8⁺ repertoires. (F) Relative usage of TCRα V and J gene segments according to cell type. The z-score for each segment is calculated from the distribution of normalized proportions stratified by the cell type and sample. P value is calculated by comparing z-scores in CD4⁺ T and CD8⁺ T cells using t test, and false discovery rate (FDR) is calculated using Benjamini-Hochberg correction: *P < 0.05, **FDR < 10%. Gene names and asterisks are colored by significant enrichment in CD4⁺ T cells (blue) or CD8⁺ T cells (orange).