Thymuses were harvested from C57BL/6 mice at 6–7 weeks of age. All experiments were approved by the St Vincent’s Hospital Animal Ethics Committee and performed under the Australian code for the care and use of animals for scientific purposes.

Total thymocytes were obtained by crushing the thymus through a metal sieve to generate a single cell suspension. The cells were washed with PBS and filtered through a 70 µm sieve to remove any clumps. CD4⁺- and CD8⁺-expressing thymocytes were depleted using anti-CD4 and CD8 magnetic-activated cell-sorting (MACS) beads (Miltenyi Biotec), according to the manufacturer’s instructions. The depleted thymocyte preparation was stained with surface antibodies for sorting on a FACSAria (BD Biosciences).

For the analysis of cell surface phenotype, cells were simply stained with antibodies. For the analysis of intracellular cytokine expression, cells were first restimulated in vitro with 50 ng/mL PMA + 2 µg/mL ionomycin (both Sigma-Aldrich) in the presence of Monensin (BD Biosciences) at 37°C for 2.5 h before staining with cell surface antibodies. The cells were then fixed with the Intracellular Fixation and Permeabilization buffer set (eBioscience) and stained with antibodies against cytokines. All antibodies were purchased from eBiosciences except for the anti-CD53 and TCR Vγ1.1 antibodies, which were purchased from BD Biosciences. The full list of antibodies can be found in Supplementary Table 1 . Flow cytometry data were acquired on LSR Fortessa III (BD Biosciences) and analyzed with the FlowJo v10.7.0 (Treestar) software. When only analyzing cell surface phenotype, dead cells were excluded using DAPI. For intercellular cytokine analyses, the cells were not stained with DAPI but were gated on live cells determined by size.

Sorted thymocytes were counted and checked for viability, then loaded onto the Chromium platform (10x Genomics) for scRNAseq library construction using the Single Cell V2 or V3.1 Reagent Kit according to the manufacturer’s instructions. Libraries were sequenced using 150-cycle/150-bp-reads NextSeq500 (Illumina) or 300-cycle/150-bp-reads Novaseq (Illumina). Sequencing files were demultiplexed and aligned to the Mus musculus transcriptome reference (mm10), and count matrix was extracted using the CellRanger Single Cell software v2.1.1 or v4.0.0 (10x Genomics) (22).

The Illumina sequencing output was pre-processed with Seurat (v2.3 or v.3.2.2) on R (v3.6.3 and 4.1.0). Cells with <500 genes, >5000 genes, or >7% mitochondria gene expression were filtered out as low-quality cells. Following normalization and removal of confounders, highly variable genes were identified and selected using the VST selection method (23). Unsupervised linear principal component analysis (PCA) was performed on these highly variable genes to group them into 20 principal components. Cell clustering was implemented using the number of components that retain >90% of variance of gene expression in the data.

DoubletFinder (v2.0.3) was applied to remove likely sequencing doublets before downstream analyses (24). The expected number of doublets was calculated as 0.75% of recovered cells. The remaining cells were re-clustered and visualized with t-distributed stochastic neighbor embedding (t-SNE) or uniform manifold approximation and projection (UMAP) dimensional reduction. Differential expression between subclusters was carried out using the FindAllMarkers function, with default parameters; differentially expressed genes with adjusted p-value <0.05 and fold change >0.5 or <-0.5 (log2FC) were considered unless otherwise stated.

Cell cycle genes specific to the G1, S, or G2/M stages were used to perform cell cycle scoring and assign cells to their respective stage of the cell cycle (25). Cell cycle genes were regressed out using Seurat’s built-in regression model.

To compare cell types and proportions across three independent sequencing runs, the datasets were integrated as described at https://satijalab.org/seurat/archive/v3.0/integration.html (26). The Seurat package (v.3.2.2) was used to assemble multiple distinct scRNAseq datasets into an integrated dataset and cell cycling scores were calculated. To remove technical variability, the datasets were pre-processed and normalized using SCTransform (27). To correct for experimental batch effect, integration anchors were identified between the experiments then merged using canonical correlation analysis. Linear dimensional reduction was applied and principal components that retain >90% of variance of gene expression in the data were included for downstream analysis. Unsupervised clustering was implemented on the integrated data.

Filtered 10x data was imported into Monocle 2 by generating a cell dataset from the raw counts slot of the Seurat object. Cells were ordered into a branch pseudotime trajectory according to the procedure recommended in the Monocle 2 documentation (28). The highly variable genes identified by Seurat were chosen as ordering genes to recover pseudospatial trajectories using the setOrderingFilter, reduceDimension, and orderCells functions in Monocle 2 with default parameters. Differential expression between pseudotime states was determined using the Seurat function FindAllMarkers.

Slingshot (29) was also employed to infer developmental trajectories using the umap, clusterLabels = seurat_clusters, and start.clus = X functions, with DN1 cells fixed as the starting point.

Thymocyte subpopulations of interest were purified by MACS depletion followed by cell sorting and then plated onto OP9-DL1 monolayers (30). The OP9-DL1 cells were inactivated with Mitomycin C (Stem Cell) immediately prior to use. 10³ sorted thymocytes were seeded per well in 96-well plates in αMEM (Thermo Fisher) supplemented with 20% FCS (Bovogen Biologicals), penicillin/streptomycin/gentamycin (Sigma), 2 ng/mL murine IL-7 (Peprotech), and 5 ng/mL human FLT3L (Peprotech). The media were refreshed every 2 days and freshly inactivated OP9-DL1 cells were added every 4 days.

Sorted total DN1 thymocytes were pre-cultured on OP9-DL1 for 24 h. The cells were then transduced with barcode lentivirus library (31) in StemSpan medium (Stem Cell Technologies) by centrifugation at 900×g for 1.5 h at room temperature. A viral titer pre-determined to give 5–10% transduction efficiency was used to ensure that the cells were not transduced with multiple barcodes; 2 ng/mL murine IL-7 and 5 ng/mL human FLT3L were then added to each well and the cells were returned to the incubator. The following day, fresh αMEM with supplements was added. αβ (CD90.2⁺ CD8α⁺ TCRγδ^-) and γδ (CD90.2⁺ CD8α^- TCRγδ⁺) lineage cells were sorted after 14 d and 20 d of the OP9-DL1 co-culture.

Barcode library construction was performed as described previously (31). The cells were lysed in 0.5 mg/ml Proteinase K (Invitrogen) in Direct PCR Lysis Buffer (Viagen) at 55°C for 2 h. The Proteinase K was then inactivated at 85°C for 30 min and 95°C for 5 min. The lysate was split into two wells for technical replicate PCRs. A first round of PCR was performed using 1× Standard-Taq magnesium-free reaction buffer pack (NEB) with 2 mM MgCl2 (NEB), 0.2 mM dNTPs (made in house), 0.5 μM TopLiB forward primer (TGCTGCCGTCAACTAGAACA), and 0.5 μM of BotLiB reverse primer (GATCTCGAATCAGGCGCTTA) for 32 cycles (1 cycle at 95°C for 5 min, 30 cycles at 95°C for 15 s, 57.2°C for 15 s, and 72°C for 15 s followed by 1 cycle at 72°C for 10 min). A second round of PCR was then performed to add different Illumina index to each sample by amplifying the first-round PCR product with a sample specific Illumina forward index primer and a common Illumina reverse index primer for 32 cycles (1 cycle at 95°C for 5 min, 30 cycles at 95°C for 15 s, 57.2°C for 15 s, and 72°C for 15 s followed by 1 cycle at 72°C for 10 min). An aliquot of the PCR product was run on 2% agarose gel to check for barcode amplification, then the samples were pooled, and the DNA was cleaned using NucleoMag SPRI beads (Machery-Nagel).

The 75-cycle sequencing runs were performed on a NextSeq instrument (Illumina). The data was demultiplexed and aligned to the reference barcode library using the processAmplicons function from the edgeR package (32). The barcode counts were then processed in the following steps: 1) barcodes with less than two read counts in all sample/cell types were excluded from the analysis; 2) barcodes that were detected in the water control were also removed from the analysis; 3) the read count between technical replicates of the same sample was averaged and the total read count in each sample was then normalized to 100%. The normalized barcode profiles were checked for any biases by t-SNE clustering. DBSCAN (version 1.1-8) was then used to classify barcodes based on their corresponding t-SNE coordinates.

Fetal thymic lobes were isolated from embryos at gestational Day 14.5 following timed pregnancies of C57BL/6 female mice. They were cultured for 5 days on 0.8-mm isopore membranes (Millipore) atop surgical gelfoam sponge (Ferrosan Medical Devices) soaked in RPMI-1640 (Sigma) supplemented with 10% FCS (Bovogen Biologicals), 20 mM HEPES (Sigma-Aldrich), 50 mM 2-mercaptoethanol (Sigma), and 1.35 mM 2′-deoxygyanosine (dGuo, Sigma) to deplete endogenous thymocytes. The depleted thymic lobes were then transferred onto new sponges soaked in fresh media with supplements but without dGuo for 2 days before repopulation with thymocyte progenitors. To repopulate thymic lobes, they were placed in 20 mL hanging drop cultures on Terasaki plates (Sigma-Aldrich) containing 5×10² to 2×10³ sorted thymocytes for 24 h before returning to fresh sponges. The media were refreshed every 3–4 days. After 14 days, single cell suspensions of the thymic lobes were generated by passing through a 70-mm sieve for analysis by flow cytometry.

Statistical testing was performed with one-way analyses of variance (ANOVA) using Prism v9 (GraphPad). p-Values are shown as * < 0.05, ** < 0.01, *** < 0.001, and **** < 0.0001 where each statistical significance was found, and all data are represented as means ± S.E.M.

All scRNAseq datasets have been deposited in the NCBI Gene Expression Omnibus repository under accession number GSE188913.

To better characterize the heterogeneity of the early stages of T cell development, we analyzed the transcriptional landscape of DN and γδ thymocytes at single-cell resolution. Cells were sorted from the thymus of C57BL/6 mice for analysis by 10x scRNAseq over three independent runs ( Figures 1A, B ; Table 1 ; Supplementary Figure 1A ). The first run consisted of total DN (defined as CD4^-CD8^-B220^-CD11b^-CD11c^-NK1.1^-TCRβ^-) and TCRγδ⁺ thymocytes, the second consisted only of DN1 and DN2 (CD4^-CD8^-CD44⁺B220^-CD11b^-CD11c^-NK1.1^-TCRβ^-TCRγδ^-) cells, and the third involved sorting DN1+DN2, DN3 and TCRγδ⁺ cells separately and mixing back together post-sort at a ratio of 55% to 30% to 15%, respectively. The latter two runs were performed to ensure that a sufficient number of DN1 and DN2 cells were captured for high-resolution analysis that is not possible by analyzing total DN cells.

A total of 22,094 high-quality cells passed quality control checks across the three datasets. The datasets were first integrated by anchoring common cells (26) in order to assemble a global view of early T cell development. To recover biological distinction from the different replicates and minimize batch-associated variability, the pooled data was normalized using SCTransform. Following dimensional reduction, unsupervised clustering was performed using the first 12 PCs. Next, we performed clustering of the cells at different resolutions starting from 0.5 up to 3.0. We selected a minimum value that was sufficient to separate pre- and post-T lineage commitment cells (DN2a vs. DN2b) and pre- and post-β-selection cells (DN3a vs. DN3b) into different clusters; 2.0 was determined to be the minimum resolution. This identified 30 distinct clusters ( Figure 1C ; Supplementary Figure 1B ), which were assigned to a canonical DN stage or γδ thymocytes based on the expression of key marker genes ( Figure 1D ; Supplementary Figure 2 ). High Cd44 and Il7r expression but low T lineage gene expression, including Il2ra, Tcf7, Cd24a, Notch1, and Bcl11b, identified DN1 cells. DN2 cells were identified by upregulation of T lineage genes and downregulation of Il7r. DN3 cells were identified by Ptcra and further upregulation of T lineage genes. Low Cd8b1 and loss of Il2ra distinguished DN4 from DN3 cells. High Trdc, Id3, Sox13, and Rorc identified γδ thymocytes. This indicates that multiple subpopulations correspond to each of the canonical DN stages.

Because there is a massive expansion of cell number between the ETP and DP stages, we also tested the effect of cell cycle gene expression on the clustering. Cell cycling scores were calculated from the integrated data and regressed using Seurat’s built-in regression model ( Supplementary Figure 3 ). There was a slight variation in the number of output clusters, but we did not observe any biological variability that could be explained by cell cycle status. The genetic profiles were nearly identical between the outputs with cell cycle genes left in and regressed out. Thus, the transcriptional heterogeneity of DN thymocytes observed is not simply a result of being in a different phase of the cell cycle. This also meant that exclusion of cell cycle genes was not necessary for downstream analyses.

We next employed Monocle 2 (28) to infer a potential developmental pathway from the scRNAseq data of DN and γδ thymocytes by ordering the cells based on tracking gene expression in pseudotime analysis ( Figure 2A ). This assembled the cells along an asymmetric trajectory that divided into six states ( Figure 2B ). Each state was then analyzed for the expression of signature genes ( Figure 2C ) to assign to a stage in development. State 2, comprising DN1 and some DN2a cells, was identified as pre-T lineage committed cells, and therefore the starting point. State 3 contained T lineage committed cells and cells that had passed β-selection, and therefore corresponded to the main αβ pathway, with DN4 as the endpoint. States 5 and 6 contained T lineage committed cells but terminated with DN3 cells. State 1 corresponded to the γδ branch. Thus, based simply on the transcriptomic profile of the cells, Monocle 2 appears to infer that γδ cells develop directly from DN1.

We also performed trajectory analyses with another algorithm, Slingshot (29), on the integrated dataset ( Figure 2D ). This also suggested that γδ cells develop directly from DN1. Moreover, the main αβ pathway, ending with DN4, and the alternate branch that terminated with DN3 cells were also observed. However, this algorithm also predicted a DN1 branch.

We suspected that the alternate branch that terminated with DN3 cells represented cells that had not passed β-selection. We therefore performed differential gene expression analysis between the State 4/5/6 cells with the State 3 cells. First, Cd27 expression was upregulated in State 3 but not State 4/5/6 cells, which was previously shown to delineate pre- and post-β-selected DN3 cells (2). Furthermore, there was an enrichment of differentially expressed genes associated with “Apoptosis”, “Senescence”, “DNA replication” and “Cell Cycle” KEGG terms ( Figure 2E ), among others. Together, these strongly suggest that this branch represents cells that had failed β-selection.

Our de novo assembly of a model of early T cell development based on scRNAseq data thus suggests that the current model where γδ development branches from the DN2b>DN3a stage is likely to be incomplete.

If the αβ and γδ lineages can develop independently, we might expect to see each lineage being derived from distinct DN1 thymocytes. To investigate this, we sorted total DN1 (CD44⁺CD25^-CD4^-CD8^-B220^-CD11b^-CD11c^-NK1.1^-TCRβ^-TCRγδ^-) thymocytes, which includes both CD117⁺ ETPs and the CD117^-/lo subpopulations. The cells were tagged with a lentiviral library of unique DNA barcodes (31). Once tagged, that barcode is inherited by the progeny of that cell and thus we can estimate how frequently αβ cells and γδ cells originate from the same starting DN1 cell ( Figure 3A ). If an αβ cell and a γδ cell inherit the same barcode sequence, it means that they were derived from the same progenitor.

First, we performed a time-course of αβ vs. γδ differentiation from DN1 thymocytes to determine the optimal times for this analysis ( Figures 3B, C ). Cell surface CD8 was used as the marker of αβ lineage cells because the fixation required to detect intracellular TCRβ interferes with downstream analyses. Late DN4 to DP-staged cells and the few CD8 (TCRβ⁺) single positive cells that are generated were captured by this strategy. We also checked CD4 expression (not shown), which correlated entirely with the DP stage. CD4 (TCRβ⁺) single positive cells do not develop in these cultures due to lack of class II MHC presentation (33). TCRγδ⁺ cells first appeared in these DN1 OP9-DL1 cocultures on Day 9 and reached a maximum at Day 12. CD8⁺ cells started appearing on Day 14 and reached a maximum at Day 20. Beyond this time point, mature cells started dying off. We thus chose Days 14 and 20 to sort the two lineages for barcode analysis.

The DNA from the resulting αβ and γδ cells were amplified by PCR then Illumina-sequenced ( Figure 3D ). We found that at Day 14, only 45% of the detected barcode species were sequenced in both αβ and γδ populations, while at Day 20, only 33% of detected barcode species were sequenced in both. At both time points, the majority of unique barcode species were detected in only the αβ or γδ population, suggesting that these were derived from DN1 thymocytes that gave rise to only one lineage. Together with the trajectory analysis of the scRNAseq data, this suggests that there are different DN1 thymocyte subpopulations, some of which have bipotency for both lineages, while others appear to preferentially develop into either αβ cells or γδ cells. Thus, it is possible that at least some γδ cells can develop directly from DN1 thymocytes rather than following the canonical pathway and bifurcating from αβ cells at the DN2>DN3 stage.

To better delineate the earlier stages of the developmental model, we focused the second 10x run that was performed specifically on DN1 and DN2 thymocytes. This allowed for more precise delineation of these two stages than if analyzing all DN thymocytes together. Unsupervised clustering of the 8,851 cells that pass quality control checks identified 26 clusters ( Figure 4A ; Supplementary Figure 4A ). This was based on a resolution of 2.0, which was the minimum number that resulted in clusters containing only DN1, DN2a, or DN2b thymocytes. Specifically, we checked that no cluster contained a mixture of cells with DN2a (pre-T-specification) or DN2b (post-T-specification) identity, at least based on expression of conventional stage marker genes. DN1 cells expressed high levels of progenitor markers, including Hhex, Il7r and Cd44 but low levels of T-commitment markers, such as Il2ra, Tcf7, Notch1, Cd24a, and Myb ( Supplementary Figure 4B ). Late T-lineage commitment genes, including Bcl11b, Rag1, Rag2, Notch3, and Ptcra, were used to separate the DN2a and DN2b subpopulations ( Supplementary Figure 4B ). This resulted in eight DN1, five DN2a, and 10 DN2b subpopulations. Some clusters formed distinct subpopulations, particularly DN1 thymocytes, while others appeared to be divisions within a continuum, particularly DN2 thymocytes. Such heterogeneity within these early T cell developmental stages is consistent with that previously reported, at least for DN1 thymocytes (3). There were also three small clusters of non-T cells consisting of doublets and B cells that, for simplicity, were excluded from downstream analysis.

While various cell surface markers have been used to divide DN1 thymocytes into DN1a to DN1e, and DN2 thymocytes into DN2a and DN2b (3), these are insufficient to delineate the larger number of subpopulations that are suggested by the scRNAseq analysis. We therefore needed to identify useful cell surface markers that could be used for flow cytometry. Differential expression analysis was performed to select features (p < 0.05 and twofold difference) for cross-referencing with GO terms for cell surface proteins (not shown). We then tested commercial antibodies against these candidates. Protein does not always correlate with mRNA levels and therefore not all antibodies produced a staining pattern that matched the scRNAseq data. However, we were able to identify a minimal panel that could delineate the eight DN1 subpopulations and most of the DN2 subpopulations.

CD24 and CD117 have previously been shown to divide DN1 thymocytes into five subpopulations (3). DN1a, DN1b, DN1c, and DN1d each corresponded to a single cluster in the scRNAseq analysis, but we identified four clusters (#1, 3, 5, and 22) that corresponded to DN1e thymocytes ( Figure 4A ). The inclusion of CD314 (NKG2D), CD317 (BST2), and Sca-1 delineated these four DN1e subpopulations ( Figures 4B, C ).

Differential CD90 and CD117 expression distinguished DN2a from DN2b cells, which were then subdivided further based on CD53, Ly-6d, and CD3e expression ( Figures 4D, E ). However, unlike DN1 cells that appeared to form distinct subpopulations, DN2 cells did not clearly partition into subpopulations. The pair of DN2a clusters 7/15 could not be separated by cell surface markers, nor could the pairs of DN2b clusters 4/23, 9/13, and 8/19 due to very similar gene expression profiles. Thus, these cluster pairs were combined. The resulting 11 DN2 clusters were annotated as DN2a-1 to 4 and DN2b-1 to 7 for identification purposes (not an indication of developmental order).

Using this panel of antibodies, we reliably identified the eight subpopulations of DN1 cells and 11 DN2 cells over experiments ( Supplementary Figure 5 ), and thus, it was employed to sort subpopulations of DN1 and DN2 thymocytes for functional analyses.

The transcriptional heterogeneity of the DN1 thymocytes may be an indication that only some subpopulations are true progenitors of T cells, as previously suggested (3). Alternately, these different subpopulations could represent different progenitors of αβ and γδ lineages, which would be consistent with the early branch point inferred by pseudotime trajectory analyses ( Figures 2A, B ). To test this, we sorted the DN1 and DN2 subpopulations using our antibody panels and accessed their αβ vs. γδ lineage potential in OP9-DL1 cultures after 14 and 20 days.

At Day 14, there were primarily αβ lineage or undifferentiated cells in DN1a, DN1b, and DN1c cultures, with very few TCRγδ⁺ cells produced ( Figures 5A–C ). DN1d and all DN1e cultures only generated TCRγδ⁺ cells. Notably, the number of TCRγδ⁺ cells produced in all DN1e cultures was substantially greater than the number of TCRγδ⁺ produced in DN1a and DN1b cultures ( Figure 5C ).

At 20 days, there was a substantial increase in the percentage and number of αβ cells produced in DN1a and DN1b cultures, while DN1c had given rise to both γδ cells and αβ cells. DN1d and DN1e subpopulations continued to exhibit a bias toward the γδ lineage, with DN1e-1 and DN1e-2 producing the highest percentage and number of TCRγδ⁺ cells. In terms of absolute numbers, DN1a and DN1b subpopulations produced αβ cells at a greater rate than the production of γδ cells from the DN1c, DN1d, or DN1e subpopulations. Starting from 10³ DN1a thymocytes, almost 2 × 10⁴ αβ cells were present after 20 days, whereas only 2–3×10³ γδ cells had been produced in all four DN1e cultures after this time ( Figure 5C ). This of course could be a reflection of γδ cells differentiating and dying more quickly, but is also implies that DN1a and DN1b subpopulations are proliferating more quickly (3) and preferentially producing αβ cells.

At Day 14, all DN2b subpopulations and DN2a-1 had produced a high percentage and number of αβ cells and few TCRγδ⁺ cells, while the rest of DN2a subpopulations produced mostly TCRγδ⁺ cells and some αβ cells ( Figures 5D, E ). By Day 20, all DN2 subpopulations had produced αβ lineage cells, while TCRγδ⁺ cells had been lost from the DN2a cultures. There was also a dramatic reduction in cell numbers in the DN2b cultures, which may be a result of the cells undergoing cell death after reaching the DP stage because these cultures poorly support the later stages of αβ maturation (33).

While OP9-DL1 cultures are a well-characterized system for analyzing T cell development, it is possible that the lineages biases observed for the different DN1 subpopulations may be exaggerated in this model. We therefore also assessed the differentiation of select DN1 subpopulations in fetal thymic organ cultures (FTOCs), by seeding the sorted subpopulations into dGuo-depleted fetal thymic lobes ( Figure 5F ). The lobes were analyzed at Day 14 after repopulation. Like in the OP9-DL1 co-cultures, DN1b cells preferentially produced αβ cells, DN1c cells primarily produced αβcells, and some γδ cells, while DN1d and DN1e-4 cells only produced TCRγδ⁺ cells ( Figure 5G ).

Thus, both OP9-DL1 and FTOC systems suggest that the DN1d and DN1e subpopulations are biased toward the γδ lineage. Moreover, although DN1a, DN1b, and DN1c can produce both αβ cells and γδ cells, they are heavily biased toward the αβ lineage and produce these cells in large numbers. The DN1d and DN1e subpopulations, together, make up more than three-quarters of DN1 thymocytes. Thus, on a per cell basis, it appears that more γδ cells are produced from CD117^- DN1 thymocytes than from what are traditionally considered the ETPs (DN1a and DN1b).

Consistent with two group of γδ progenitors identified by Sagar and colleagues (34), γδ thymocytes clearly segregate into two clusters (#11 and 19) in our scRNAseq dataset ( Figure 1C ). We therefore wanted to determine how these relate to the different DN1 subpopulations that produced γδ cells. Differential gene expression analysis revealed substantial differences between the two γδ clusters, with 93 genes expressed at significantly higher levels by cluster 11 cells, while 117 genes were expressed at significantly higher levels by cluster 19 cells ( Figure 6A ). Genes that were highly expressed by cluster 11 include Gzma, Blk, Maf, Sox13, Etv5, Gata3, Ccr9, Rorc, Sox4, Tcf12, Lgals9, Cmak4, and Bcl11b, which are all associated with the IL-17 effector phenotype (34). This suggests that cluster 11 cells are probably the γδ thymocytes that mature into the γδ17 subset in the periphery. Numerous interferon-related genes, such as Stat1, were highly expressed by cluster 19 cells, suggesting that these cells probably mature into the IFNγ-producing γδ subset. This is consistent with a previous study that suggested the eventual effector function of γδ T cells may already be acquired in the thymus (21).

Next, to investigate the relationship between these γδ thymocyte subpopulations and the DN1 subpopulations that differentiate into TCRγδ⁺ cells, the DN1c, DN1d, and DN1e subpopulations were analyzed for the expression of the 210 differential genes that distinguished the two γδ thymocyte populations. This revealed a similar transcriptional profile between DN1d and cluster 11 γδ thymocytes, while DN1c and all the DN1e subpopulations overlapped significantly with cluster 19 γδ thymocytes ( Figure 6B ). However, there were clearly also differences between the DN1e subpopulations.

Differential expression of transcription factors was notable. These are likely to be important because as regulators of gene expression they could potentially play key roles in hardwiring γδ effector outcomes in DN1 thymocytes. Sox13 was highly expressed by DN1d cells ( Figure 6B ), which was previously shown to be important for the differentiation of a subset of DN1 thymocytes into IL-17-producing γδ T cells (18). Interestingly, DN1e subpopulations also express transcription factors associated with IL-17-producing γδ T cells. Notably, DN1e-1 and DN1e-2 cells expressed high levels of Maf, while Gata3 was highly expressed by both DN1e-1 and DN1d cells ( Figure 6B ). We thus predict that the foundation of the γδ effector transcriptional network may already be in place in DN1 subpopulations, with DN1d going on to develop into IL-17-producing cells and DN1e potentially producing both effector subsets.

To determine if different DN1d and DN1e subpopulations might differentiate into distinct effector γδ T cell subsets, the TCRγδ⁺ cells that developed in OP9-DL1 co-cultures were analyzed for intracellular IL-17A and IFNγ expression and for expression of specific Vγ chains ( Figure 7A ). Cytokine production is highly associated with specific Vγ chain usage, with Vγ1.1⁺ cells enriched for IFNγ and Vγ2⁺ cells enriched for IL-17A production (19, 35).

DN1c thymocytes generated both Vγ1.1⁺ and Vγ2⁺ cells, with only a low percentage of Vγ1.1⁺ cells expressing IFNγ ( Figures 7B–D ). DN1d primarily produced Vγ2⁺ cells that express IL-17A. Similarly, DN1e-1 and DN1e-2 were biased towards a Vγ2⁺ IL-17A⁺ γδ effector subset. On the other hand, only DN1e-3 and DN1e-4 exhibited the plasticity to generate both IL-17A and IFNγ-expressing cells ( Figures 7B–D ). This suggests that the foundation of γδ effector programs is already in place in DN1 thymocytes. Indeed, the differential gene expression analysis clearly showed that each of γδ-producing DN1 subpopulations are transcriptionally distinct from each other, although DN1e-3 and DN1e-4 cells did appear somewhat similar to each other ( Figure 6B ).