Mapping the two distinct proliferative bursts early in T-cell development

Abstract

T-cell development occurs in the thymus and is tightly regulated to produce a diverse enough repertoire of mature T cells that can recognize any potential pathogen. The development of T cells is dependent on small numbers of uncommitted precursors that continually seed the thymus from the bone marrow. As they progress along the developmental pathway, there is a massive expansion in cell number to generate the necessary diversity in T-cell receptor chain usage. It is recognized that there are two proliferative bursts that occur early in T-cell development, one prior to β-selection and one after, and these are responsible for the expansion. While the proliferation following β-selection is well-characterized, the earlier proliferative burst has yet to be precisely defined. In this study, we employ single-cell RNA sequencing coupled to trajectory inference methods to pinpoint when in T-cell development thymocytes are induced into cell cycle. We show that the first proliferative burst is initiated in the double-negative (DN) 2a stage before T lineage commitment occurs, with cell cycling downregulated by the DN3a stage. A second burst is then initiated at the DN3b stage, immediately after β-selection. We subsequently employ fluorescence-activated cell sorting-based analysis for DNA content to confirm these two proliferative bursts.

Keywords: Cell cycle; T-cell development; proliferation; single-cell RNA sequencing.

Figure 1

Clustering and trajectory analysis of DN thymocyte scRNA‐seq data. (a) Three independent data sets of mouse DN thymocytes were analyzed: (1) total DN and TCRγδ⁺ cells, (2) only DN1 and DN2 cells and (3) sorted populations of DN1/DN2, DN3 and TCRγδ⁺ cells that were mixed back together after sorting at a ratio of 55% to 30% to 15%, respectively. See Supplementary figure 2 for the sorting strategy. Following processing of the 10× data in Cell Ranger, doublets were removed with DoubletFinder and cell cycle genes were regressed out. (b) The integration of the three data sets with Seurat using anchors following preprocessing and batch corrections implemented with SCTransform. Cell cycle genes were also regressed out. Clustering of each independent data set or the integrated data set was performed at a resolution of 2.0, which resulted in DN1, DN2, DN3 and DN4 thymocytes separating into distinct clusters, based on the expression of the marker genes indicated in Supplementary figure 3. Indicated are the clusters that corresponded to the broad DN1/2/3/4 stages in T‐cell development, TCRγδ⁺ cells or non‐thymocytes. (c) Contribution of each independent data set to the integrated data set. (d) Developmental trajectory inferred with Monocle 3. The individual cells are colored coded by position in pseudotime. (e) Color coding of the Monocle 3 trajectory by DN stage. DN, double negative; scRNA‐seq, single‐cell RNA sequencing; TCR, T‐cell receptor; Thy, thymocytes; t‐SNE, t‐distributed stochastic neighbor embedding.

Figure 2

Mapping the two proliferative bursts during early T‐cell development using scRNA‐seq data. (a) The integrated data set, (b) data set 2 (DN1 + DN2 thymocytes) and (c) data set 3 (DN1/DN2, DN3 and γδ thymocytes mixed after sorting) were first clustered. Data sets 2 and 3 were analyzed individually to achieve a greater resolution of the DN2 stage. The assignment of cells within each cluster to cell cycle phase was then performed in Seurat. Each line indicates the percentage of cells assigned to G1/0, S or G2/M phase within an individual cluster. The clusters within each analysis are ordered based on the trajectories inferred by Monocle 3 in Figure 1d and Slingshot in Supplementary figure 4b. Clusters that correspond to β‐selection failure are shown separately. DN, double negative; scRNA‐seq, single‐cell RNA sequencing.

Figure 3

FACS analysis for DNA content confirms that thymocytes cell cycling is upregulated at the DN2a stage and again at the DN3b stage. DN thymocyte subpopulations were identified with the gating strategy shown in Supplementary figure 5 and DAPI staining for DNA content was employed to quantify cells in G1/0, S or G2/M phase. (a) An example of DNA content analysis of DN subpopulations. (b) The mean ± s.e.m. of three experiments, with each dot indicating an experimental datapoint. The cells from up to three mice were pooled for each experiment. DAPI, 4′,6‐diamidino‐2‐phenylindole; DN, double negative; FACS, fluorescence‐activated cell sorting.

Figure 4

Upregulation in cell cycling at the DN2a and DN3b stages is associated with an increased expression of metabolic genes. The gene expression profiles of total DN1 versus DN2a, and DN3a versus DN3b were extracted from the scRNA‐seq data. Differential expression analysis was then performed for genes associated with (a) glycolysis (KEGG map 00010) or (b) TCA cycle (KEGG map 00020). Shown are only genes with significantly different expression (adjusted P < 0.05) between DN1 and DN2a or between DN3a and DN3b. The genes are ordered by largest to smallest average fold change. The genes found to be upregulated from the DN1 to DN2a stage were analyzed for transcription factor recognition motifs with RcisTarget. (c) An enrichment of c‐MYC recognition sites was identified in the upregulated glycolysis gene set, (d) while an enrichment of GATA recognition sites was identified in the upregulated TCA cycle gene. Shown are the recovery curves for each analysis, where the recovery of each specific transcription factor motif is indicated by the blue line. The average recovery of all motifs is indicated by the red line, while average recovery + standard deviation is indicated by the green line. Also shown are the NES for each analysis, which is calculated based on the AUC distribution for all motifs in the gene set. A score for > 3.0 is considered significant. AUC, area under the curve; DN, double negative; KEGG, Kyoto Encyclopedia of Genes and Genomes; NES, normalized enrichment score; scRNA‐seq, single‐cell RNA sequencing; TCA, tricarboxylic acid.