While many studies have explored the heterogeneity of IELs using single-cell sequencing techniques, it has been a challenge to clearly distinguish between TCRαβ⁺ and TCRγδ⁺ cells at the transcriptional levels, in part due to the shared usage of Vα and Vδ gene locus segments. To better characterize the populations of TCRαβ⁺ and TCRγδ⁺ IELs, we performed scRNAseq on FACS-sorted TCRαβ⁺ and TCRγδ⁺ IELs collected from the small intestine of naïve healthy 6-8-week-old C57BL/6 mice from the Jackson Laboratory. Uniform manifold approximation and projection (UMAP) analysis was performed separately for the TCRαβ⁺ and TCRγδ⁺ cells, revealing 16 and 17 clusters, respectively ( Figures 1a, b , Supplementary Figures 1a, b ). Based on the expression level of Cd4, Cd8a, and Cd8b1 coreceptor genes, we noticed that minor IEL populations defined in flow cytometry experiments (TCRαβ⁺ DN, TCRαβ⁺ DP, and TCRγδ⁺ DN; Supplementary Figure 2a ) were not represented as distinct clusters ( Figures 1c, d ), possibly due to the relatively small population sizes or stringent quality control filtering. Among the conventional TCRαβ⁺ cells, CD4⁺ cells comprised clusters 9, 14 and part of cluster 3, while CD8αβ⁺ cells consisted of clusters 2, 8, 12 and part of cluster 3. Surprisingly, our analysis of natural TCRαβ⁺ CD8αα⁺ cells revealed two distinct cluster groups: clusters 1, 4, 6, 15, and clusters 5, 7, 10 ( Figures 1a, c ). Similarly, the TCRγδ⁺ CD8αα⁺ clusters formed diverging groups on the UMAP plot ( Figures 1b, d ).

To define the phenotype of each cluster group, we summarized the expression level of marker genes ( Supplementary Figures 2b, c ). The conventional TCRαβ⁺ CD4⁺ clusters consist of T-follicular helper-like cells expressing Cxcr5 and Il21 (cluster 9), Foxp3 ⁺ Tregs (cluster 14), and naïve cells expressing Sell, Ccr7, and Dapl1 (cluster 3) (40–43). The CD8αβ⁺ clusters 2 and 8 were enriched for effector-like T cell markers, including Prdm1 (which encodes the transcription factor BLIMP1) and effector molecules Gzmb and Gzmk (38) ( Figure 1e ). In contrast, cluster 12 was enriched for memory-like T cell markers, such as Tcf7 (encoding the transcription factor TCF1) and Id3. CD8αβ⁺ cells were also detected in the naïve cell cluster (cluster 3). Interestingly, within the TCRαβ⁺ CD8αα⁺ clusters, we also observed two distinct groups: clusters 1, 4, 6 and 15 were enriched for effector-like markers, including Prdm1, Tyrobp, Gzma and Gzmb, whereas clusters 5, 7 and 10 were enriched for memory-like markers Tcf7 and Id3 (37, 38) ( Figure 1e ). A similar separation between effector-like and memory-like profiles was found among the TCRγδ⁺ cells (clusters 1, 3, 4, 6, 7, 8, 11 and clusters 2, 5, 10, 13, respectively; Figure 1f ).

Collectively, we annotated cluster groups according to the expression of enriched transcription factors: the effector-like cluster groups were labeled as “Prdm1 ⁺” clusters, and the memory-like cluster groups were labeled as “Tcf7 ⁺” clusters ( Figures 1g, h ). Additionally, a unique effector-like population within the TCRγδ⁺ IELs characterized by Zeb2, Gzmk, and Ccr5, was labeled as “Zeb2⁺ ” cells (former clusters 9 and 15). Finally, cluster groups enriched in proliferation and cell cycle genes (Neil3, Bub3, Esco2, Aspm, Mxd3 for TCRαβ⁺ dataset and Dtl, Ung, Chek1, Cdc6 and Mcm10 for TCRγδ⁺ dataset) were annotated as “Mki67 ⁺” clusters, corresponding to former TCRαβ⁺ cluster 11 and TCRγδ⁺ clusters 12, 14, 17. Clusters with poorly defined identities or nonspecific cluster localization (former TCRαβ⁺ clusters 13 and 16, and TCRγδ⁺ cluster 16) were excluded from further analysis. In our final annotation, we retained 4,845 TCRαβ⁺ IELs and 10,418 TCRγδ⁺ IELs.

TCRαβ⁺ and TCRγδ⁺ CD8αα⁺ IELs are considered “non-conventional” T cells, due to their unique developmental pathways, TCR-independent activity, and other innate-like features. These cells play versatile roles in gut immunity, including the secretion of antimicrobial peptides during infections, exhibiting cytotoxic potential, and supporting repair after inflammation, while also displaying higher activation thresholds for self-reactivity (1, 3, 44, 45). Our observation of the distinct CD8αα⁺ subsets may provide insight into their functional diversity.

To better understand the differences between the CD8αα⁺ sub-clusters, we filtered out cells that expressed Cd4 or Cd8b1 coreceptor genes and re-clustered the populations ( Supplementary Figures 3a, b ). When confirming the expression of Cd8a, we noticed the re-clustered cells included those with no Cd8a expression in each of the annotated populations ( Supplementary Figures 3c, d ), The coreceptor negative cells contributed to 4.7% and 11.9% of the re-clustered TCRαβ⁺ and TCRγδ⁺ cells, respectively, and were localized to the edges of the UMAP plots ( Supplementary Figures 3e, f ). We speculated that these cells represent the population of DN IELs identified by flow cytometry experiments ( Supplementary Figure 2a ), which are known as precursors to CD8αα⁺ IELs (13). Accordingly, we performed the subsequent analyses separately for the Cd8a ⁺ (CD8αα⁺) IELs and coreceptor negative DN IELs, with the cluster labels identified in Figure 1 .

We first compared the transcriptional profiles of the CD8αα⁺ sub-clusters ( Supplementary Figure 4 , Supplementary Table 1 ). In both the TCRαβ⁺ and TCRγδ⁺ datasets, Prdm1 ⁺ clusters were significantly enriched for effector T cell-like transcription factors including Prdm1, Runx1, Runx2, Plek and Irf8. In contrast, Tcf7 ⁺ clusters showed marked enrichment for memory-like transcription factors such as Tcf7, Id3, Aff3, Bach2, Nfkb1 and Satb1 ( Figures 2a–d ). Additionally, several transcription factors displayed lineage-specific expression pattern. For example, Foxo3, Jazf1, Zbtb16 and Lef1 were differentially expressed among clusters within the TCRαβ⁺ lineage, whereas Zmat4 and Zfnx3 and Elk4 were specific for TCRγδ⁺ cells. Interestingly, within the TCRγδ⁺ dataset, the Zeb2 ⁺ cluster exhibited a transcriptional profile intermediate between Prdm1 ⁺ and Tcf7 ⁺ clusters, while also expressing unique TFs such as Zeb2, Bcl11b and Rora ( Figure 2d ). Effector-like and memory-like transcriptional profiles of the Prdm1 ⁺ and Tcf7 ⁺ clusters, respectively, were corroborated by unbiased gene ontology (GO) analysis. GO terms associated with activation, effector function and cell motility were enriched in Prdm1 ⁺ clusters in comparison to Tcf7 ⁺ clusters in both TCRαβ⁺ and TCRγδ⁺ T cell populations. In contrast, GO terms related to ribosome biogenesis, which have been associated with memory formation (46, 47), were predominantly overrepresented in Tcf7 ⁺ clusters relative to Prdm1 ⁺ clusters ( Supplementary Figures 4a, b ). The Zeb2 ⁺ population appeared more similar to TCRγδ⁺ Prdm1 ⁺ cluster and was enriched in effector-related GO terms ( Supplementary Figures 4c, d ). The CD8αα⁺ Prdm1 ⁺ and Tcf7 ⁺ clusters differentially expressed genes encoding surface proteins, with Il2rb (CD122) being enriched in the Prdm1 ⁺ clusters, and Cd160 (CD160) predominantly expressed in Tcf7 ⁺ clusters ( Figures 2e, f ). Both CD122 and CD160 have been previously studied for their roles in IEL maintenance, cytotoxic activity, or protective ability, but their heterogeneous expression within total IEL populations has not been fully clarified (44, 48–51). To validate whether CD160 and CD122 can be used to distinguish memory-like and effector-like subsets of CD8αα⁺ IELs, we performed flow cytometry on small intestinal IELs stained for these surface markers. The analysis revealed two distinct populations, with approximately 20% of CD8αα⁺ IELs being CD122^intCD160⁺ and 75% being CD122^hiCD160^- ( Figure 2g ). These results suggest that the effector-like and memory-like phenotypes are distinctly distributed among total IELs, with surface markers CD160 and CD122 providing a means to distinguish between these functional subsets.

The maintenance and activity of IELs can be modulated by both diet and gut microbiota, which vary along the intestinal tract (17, 18, 52–57). These variations may shape the distribution of IEL subsets to ensure appropriate immune responses in different regions of the intestine. Therefore, we next examined the composition of CD8αα⁺ IELs along five different sections of the intestine. In both TCRαβ⁺ and TCRγδ⁺ CD8αα⁺ IELs, the proportion of CD160⁺ cells gradually increased from 10% to 30% along the small intestine (duodenum, jejunum, ileum), while the proportion of CD122^hi cells decreased reciprocally from 90% to 60% ( Figures 2h, i ). However, the distribution of these populations shifted drastically in the large intestine. Among TCRαβ⁺ IELs, CD122⁺ cells were nearly absent in cecum and colon, leaving one major population of CD122^intCD160^int cells ( Figure 2h ). In contrast, despite collecting very few of TCRγδ⁺ IELs in the large intestine, an equal proportion of CD122^hiCD160^- and CD122^intCD160^int populations were observed ( Figure 2i ).

To confirm whether the CD160⁺ population in all five sections of the intestine correspond to the Tcf7 ⁺ cluster IELs, we used Tcf7-GFP mice, where Tcf7-expressing cells are concurrently labeled by GFP (58). Indeed, the CD8αα⁺ CD160⁺ IELs expressed higher levels of GFP compared to CD8αα⁺ CD122⁺ cells in most tissue sections and contained a greater frequency of GFP+ cells, for both TCRαβ⁺ and TCRγδ⁺ subsets ( Supplementary Figures 5a, c ). These observations suggest that some yet unknown factors influence the abundance of CD8αα⁺ IELs along the intestine, promoting a high and similar prevalence of effector-like IELs in TCRαβ⁺ and TCRγδ⁺ cells in the small intestine but favoring the colonization of memory-like TCRαβ⁺ CD8αα⁺ IELs in the large intestine.

Given that DN IELs are considered precursors to CD8αα⁺ IELs, we sought to determine whether the distinct Prdm1⁺ and Tcf7⁺ transcriptional profiles observed in mature CD8αα⁺ IELs are already established at the DN stage. Therefore, using our dataset of coreceptor-negative DN IELs, we next examined whether these transcriptional differences could also be detected prior to CD8αα expression. Gene expression analysis across both TCRαβ⁺ and TCRγδ⁺ lineages revealed that DN subsets share many transcription factors with their corresponding CD8αα⁺ counterparts. Specifically, both Tcf7 ⁺ DN and CD8αα⁺ IELs expressed Tcf7 and Zfp706, while Prdm1⁺ subsets shared expression of Prdm1 and Runx1 ( Figures 2a–d , 3a–d ). Furthermore, we found a strong overlap in transcription factor expression between DN and CD8αα⁺ subsets within the TCRγδ⁺ lineage ( Figures 2a–d , 3a–d ). For example, Zmat4 was highly enriched in the Prdm1 ⁺ cluster; Zeb2, Runx2, Stat1, and Bcl11b were enriched in the Zeb2⁺ cluster; and Tcf7, Id3, Aff3, Satb1, Lef1, and Nfkb1 were enriched in the Tcf7 ⁺ cluster ( Figures 2d , 3d ). In addition, we identified several genes, which were specific to Tcf7⁺ DN subsets, including TFs associated with the cell cycle such as Hnrnpk, and Ybx1 (59–61) ( Figures 3a–d ), along with increased expression of ribosomal genes, a pattern also observed among the CD8αα⁺ subclusters ( Supplementary Table 1 ). These findings suggest that DN cells represent distinct populations that nonetheless relate to CD8αα⁺ cells, supporting the possibility of developmental continuum between these subsets.

We also examined the expression level of surface protein coding genes which were differentially expressed between CD8αα⁺ subclusters. Similarly to the CD8αα⁺ IELs, the Tcf7 ⁺ DN clusters had greater mRNA expression of Cd160, while the Prdm1 ⁺ DN clusters showed higher expression of Il2rb ( Figures 3e, f ). These markers also distinguished DN IELs at the protein level by surface expression. However, in contrast to the CD8αα⁺ IELs, which contained a greater proportion of CD122^hi cells ( Figure 2g ), approximately 70% of DN IELs were CD160⁺ and 20% were CD122^hi ( Figure 3g ). Based on Tcf7-GFP reporter expression, we observed a clear trend (although not always statistically significant) indicating that CD160⁺ DN cells expressed higher levels of Tcf7 and contained a higher frequency of GFP+ cells ( Supplementary Figures 5b, d ). These findings align with the predicted distribution from the sequencing data ( Supplementary Figures 3e, f ). Of note, unlike natural IELs, induced IELs were more uniform and did not clearly segregate into memory-like and effector-like phenotype based on expression of CD122 and CD160 ( Supplementary Figure 5e ). An exception was observed for TCRαβ⁺ CD8αα⁺ CD4⁺ (DP) cells, where the majority of cells displayed an effector-like CD160^-CD122^int phenotype.

Interestingly, unlike CD8αα⁺ IELs, the proportion of CD122⁺ and CD160⁺ cells along the intestine remained stable among both TCRαβ⁺ and TCRγδ⁺ DN IELs ( Figures 3h, i ). This suggests that CD8αα⁺ IEL subsets are uniquely sensitive to yet unidentified microenvironmental cues that may drive their “gradient-like” distribution along the small intestine. Notably, this unique sensitivity appears to be specific to the small intestine, as both the CD8αα⁺ and DN IELs in the large intestine (cecum and colon) are predominantly comprised of CD122^lowCD160^int cells ( Figures 2h, i , 3h, i ).

The characteristics of CD8αα⁺ IELs have been suggested to vary between TCRαβ⁺ and TCRγδ⁺ cells, based on functional studies (62). However, the clustering profiles of the small intestinal CD8αα⁺ IELs appeared strikingly similar between the TCRαβ⁺ and TCRγδ⁺ datasets ( Figures 2 , 3 ), suggesting the presence of subsets with shared phenotypes across both cell lineages. To further investigate this, we examined the similarity of IEL subsets between the two datasets.

First, we identified the gene signatures representing the CD8αα⁺ Prdm1 ⁺ (effector-like) and CD8αα⁺ Tcf7 ⁺ (memory-like) clusters in the total TCRαβ⁺ dataset. We calculated the differentially expressed genes comparing each subcluster to the rest of the TCRαβ⁺ cells and sorted the genes by their statistical significance ( Figure 4a ). The top 10 representative genes of the TCRαβ⁺ CD8αα⁺ Prdm1 ⁺ cluster included genes associated with CD8 T cell activation such as Fcer1g, Cd244a (encoding 2B4), and Clnk (a SLP-76 member), or genes associated with cell-cell contacts such as Fgl2, Frmd5, and Osbpl3 (63–68). The top 10 genes representing the TCRαβ⁺ CD8αα⁺ Tcf7 ⁺ cluster included genes associated with memory or stem-like T cells such as Id3, Batf3, Kit, and AW112010 (42), and genes associated with intestinal resident T cells such as Cd160 and Xcl1. Next, we projected the signature of each TCRαβ⁺ CD8αα⁺ subcluster onto the TCRγδ⁺ dataset to evaluate whether similarly annotated clusters showed similar enrichment patterns. Indeed, the TCRαβ⁺ CD8αα⁺ Prdm1 ⁺ signature was significantly enriched in the TCRγδ⁺ CD8αα⁺ Prdm1 ⁺ clusters ( Figure 4b ), indicating strong similarities between effector-like clusters of the TCRαβ⁺ and TCRγδ⁺ datasets. In parallel, the TCRαβ⁺ CD8αα⁺ Tcf7 ⁺ signature was highly enriched in the TCRγδ⁺ CD8αα⁺ Tcf7 ⁺ clusters, highlighting a similar correspondence between memory-like populations. Together, these findings underscore striking parallels between IEL subsets across the TCRαβ⁺ and TCRγδ⁺ lineages.

Although these findings validate the overlap between corresponding clusters, it is possible that additional parallels exist between the two datasets that were not detected using a signature projection approach. For example, neither of the signatures showed enrichment in the TCRγδ⁺ CD8αα⁺ Zeb2⁺ cluster. To further evaluate the similarities between TCRαβ⁺ and TCRγδ⁺ IEL clusters, we applied a second approach, where we integrated populations from the two datasets and re-clustered them with low resolution ( Figures 4c, d ). In addition to the CD8αα⁺ and DN IELs analyzed in Figures 2 , 3 , we also included the TCRαβ⁺ CD8αβ⁺ clusters as they represent a major portion of TCRαβ⁺ IELs ( Supplementary Figure 2a ). As expected, the integrated UMAP analysis showed that the “Tcf7 ⁺”, “Prdm1 ⁺” and “Mki67 ⁺” clusters from TCRαβ⁺ and TCRγδ⁺ cells grouped together ( Figure 4d ). Similar to the previous observation ( Figure 2d ), the TCRγδ⁺CD8αα⁺ Zeb2⁺ cluster grouped together with a portion of the TCRαβ⁺ CD8αβ⁺ Prdm1 ⁺ cells, despite these clusters representing separate lineages of IELs. Altogether, our results demonstrate that many of the TCRαβ⁺ and TCRγδ⁺ cells have similar transcriptional profiles, suggesting a potential convergence of their functional roles in intestinal epithelial immunity.

Next, we investigated precursor-progeny relationships among natural IEL subsets across TCRαβ⁺ and TCRγδ⁺ lineages. Given that conventional memory T cells can differentiate into effector T cells, we asked whether a similar relationship exists between memory-like IEL clusters and effector-like IEL clusters, and whether such transitions could be inferred from our scRNAseq datasets. Since thymic IEL-precursors are phenotypically DN cells, we also explored the possibility of a developmental trajectory from DN cells to CD8αα⁺ IELs. Upon re-clustering of TCRαβ⁺ CD8αα⁺ and DN cells, we observed that only the Tcf7 ⁺ cluster contained a substantial proportion of DN cells ( Supplementary Figures 3c, e ), suggesting a single precursor population for TCRαβ⁺CD8αα⁺ IELs. In contrast, the TCRγδ⁺ dataset contained two clusters, Tcf7 ⁺ and Zeb2 ⁺, that were both enriched for DN cells ( Supplementary Figures 3d, f ), indicating the potential existence of two distinct precursor populations for TCRγδ⁺CD8αα⁺ IELs. To test whether these hypotheses are supported by computational predictions, we performed RNA velocity analysis on combined datasets of TCRαβ⁺ CD8αα⁺ and DN cells, as well as TCRγδ⁺ CD8αα⁺ and DN cells ( Figures 5a, b ).

The TCRαβ⁺ dataset contained distinct trajectories within both the Tcf7 ⁺ and Prdm1 ⁺ clusters. Notably, Tcf7 ⁺ DN cells were positioned at the beginning of the trajectory, clustering together at the earliest point, and therefore appeared as sole precursors to more differentiated Tcf7 ⁺ CD8αα⁺ cells ( Supplementary Figure 3e , Figure 5a ). Along this trajectory, cells within the Tcf7 ⁺ cluster progressively downregulated the expression of transcription factor Lef1 and the activating receptors Klrk1 and Klrc2, while upregulating the memory-associated transcription factor Batf3 (69, 70). In contrast, the trajectory within the Prdm1 ⁺ cluster started near the Ki67 ⁺ population and was characterized by increased expression of Zbtb16, which encodes the innate-like T cell transcription factor PLZF (71, 72) ( Supplementary Figures 6a–c ). These findings indicate the presence of phenotypic gradients within TCRαβ⁺ IEL compartment. Moreover, the directionality of the RNA velocity vectors suggested a developmental transition from the Tcf7 ⁺ cluster to Prdm1 ⁺ cluster, with the Mki67 ⁺ population acting as an intermediate proliferating stage. Thus, our results support the notion that DN Tcf7⁺ cells function as precursors to natural TCRαβ⁺ CD8αα⁺ IELs ( Figure 5c ).

The TCRγδ⁺ dataset also revealed lineage trajectories originating from DN cells ( Supplementary Figure 3f , Figure 5b ). Within the Tcf7 ⁺ cluster, two potential trajectory directions were observed. The first followed a linear progression from DN cells to CD8αα⁺ cells, marked by upregulation of memory-associated genes Xcl1 and Batf3, and downregulation of genes involved in early T cell development, including Rgs10 and Eya2 (73, 74) ( Supplementary Figures 6d–f ). The second trajectory, directed toward cluster 9 ( Supplementary Figure 6d ) and included contributions from both the DN and CD8αα⁺ cells, but lacked clearly distinguishable features. In contrast, the trajectory within the effector-like clusters appeared to originate from the DN cells in the Zeb2⁺ cluster and progress toward CD8αα⁺ cells in the Prdm1 ⁺ cluster, and was characterized by high expression of effector-associated genes such as Tyrobp, Fcgr3, and Ccrl2 ( Supplementary Figures 6d–f ). To better characterize precursor-progeny relationships among TCRγδ⁺ cells, we analyzed the distribution of Vγ chains within the DN and CD8αα⁺ cell populations ( Supplementary Figures 7a–c ). Vγ7⁺ cells, the predominant intestinal TCRγδ⁺ subset (75, 76), were primarily located in the Prdm1 ⁺ cluster and in a subset of Tcf7 ⁺ cluster within CD8αα⁺ cells. Interestingly, they were largely absent from the DN cells within the Tcf7 ⁺ cluster but were present in the DN cells within the Zeb2⁺ cluster. In contrast, Vγ1⁺ cells were enriched in the Tcf7 ⁺ cluster in both CD8αα⁺ and DN populations. Vγ4⁺ cells were relatively evenly distributed across clusters, while Vγ5⁺ and Vγ6⁺ cells were virtually absent. These findings suggest that Vγ7⁺ cells have largely completed the transition from the DN to the CD8αα⁺ state, consistent with their relatively early seeding of the gut (77). In contrast, the broader distribution of Vγ1⁺ and Vγ4⁺ cells across both DN and CD8αα⁺ clusters may reflect ongoing differentiation and migration from thymic precursors into the intestinal IEL compartment (78, 79). We also cannot rule out a possibility that similar to what was observed in the TCRαβ⁺ dataset, a transition from the Tcf7 ⁺ cluster to the Prdm1 ⁺ cluster may occur in TCRγδ⁺ lineage as well, via an intermediate Mki67 ⁺ proliferative population. Collectively, our data support the hypothesis that TCRγδ⁺ IELs arise from two independent precursor subsets ( Figure 5c ).

We also considered an alternative hypothesis that the memory-like and effector-like IELs arise from separate precursors with their fates predetermined prior to maturation in the intestine. In mice, TCRαβ⁺ CD8αα⁺ IELs are known to develop from two thymic precursor populations that seed the intestine at different developmental stages and express oligoclonal TCR variable chains (21, 80–82). To test whether memory-like or effector-like IELs exhibit distinct TCR Vα-chain usage, we measured the expression of Vα2 and Vα3.2 on the surface of TCRαβ⁺ CD8αα⁺ IELs from five sections of the intestine. However, we did not observe consistent enrichment of either Vα-chain among memory-like or effector-like subsets ( Supplementary Figure 8 ), arguing against the idea that these populations originated from separate precursors. Altogether, these findings support a model in which memory-like and effector-like CD8αα⁺ IELs represent cells at different stages of shared differentiation pathways, enabling a spectrum of immune responses across the intestinal mucosa.

C57BL/6J (stock no. 000664) mice were obtained from the Jackson Laboratory. All breedings were maintained at the University of Massachusetts, Amherst. This study was performed in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All animals were handled according to approved institutional animal care and use committee (IACUC) protocols of the University of Massachusetts.

The small intestine and colon were first removed from the rest of the gastrointestinal tract, onto collection media (RPMI supplemented with 25mM HEPES, 1% L-glutamine, 1% penicillin/streptomycin, 50µM β-mercaptoethanol, and 3% FBS). Peyer’s patches lining the small intestine were removed. Tissues were cleaned by flushing out feces with collection media and rinsing in PBS. Tissue fragments were agitated at 37°C for 20 minutes in collection media containing 5mM EDTA and 1mM DTT, then further shaken in serum-free collection media containing 2mM EDTA. The suspension was washed several times in collection media, and the IELs were collected as cells that passed through a 70µm filter. Finally, IELs were resuspended in collection media containing 10% FBS.

Flow cytometry data were acquired on BD LSR Fortessa. The following monoclonal antibodies from BioLegend were used: CD45.2 (104), CD45 (30-F11), TCRβ (H57-597), TCRδ (GL3), CD8α (53-6.7), CD8β (YTS156.7.7), Biotin-CD122 (5H3) and Streptavidin- AF647. The monoclonal CD4 (GK1.5) antibody and Brilliant Stain Buffer were obtained from BD Biosciences. The monoclonal CD160 (CNX46-3) antibody was obtained from eBiosciences.

Live cells were treated with anti-CD16/32 Fc block (2.4G2, BD Pharmingen) prior to staining with antibodies against surface markers. Staining for surface proteins was performed at 4 °C for 40 min, and FACS buffer (PBS + 0.5% BSA + 0.01% sodium azide) was used for washes. Data from BD LSR Fortessa were analyzed in FlowJo™ v10.9.0 Software. IEL populations were analyzed as shown in Supplementary Figure 9 .

Cells were stained for CD45.2, TCRβ, TCRδ, and sorted using the BD FACSAria Fusion instrument (BD Biosciences). After sorting, cells were counted using a Cellometer K2 cell counter (Nexcelom Bioscience) and by manual counting via hemocytometer. Single cell gene expression profiling was performed using the Chromium Next GEM Single Cell 3′ v3.1 (Dual Index) kit. Each cell suspension was loaded onto a well of Chip G on the 10x Genomics Chromium Controller System following the manufacturer's user manual (10x Genomics). Barcoding and cDNA synthesis were performed according to the manufacturer's instructions. Qualitative analysis of cDNA was performed using the 2100 Agilent Bioanalyzer High Sensitivity assay. The cDNA libraries were constructed using the 10x Chromium Single cell 3’ Library Kit v3.1 (dual index) according to the manufacturer’s protocol. Quality assessment of final libraries was done on Qubit fluorometer using a DNA High Sensitivity assay (Thermo Scientific) and a 2100 Agilent Bioanalyzer High Sensitivity assay (Agilent Technologies). Libraries were sequenced on an Illumina NextSeq 500 using the NextSeq 500/550 Mid Output Kit v2.5 (150 Cycles) sequencing kit, with the following read length: 28 bp Read1 for the 10x cell barcode and UMI, 90 bp Read 2 for the insert, and 10 bp I7 and I5 for the sample index. Phix (Illumina) was spiked in at 1% as per kit manual recommendation (10x Genomics).

The 10x Could Analysis Cell Ranger Pipeline (Cellranger version 7.1.0) was used to align reads and generate feature-barcode matrices. Reads were aligned to Mus musculus reference genome (Mouse GRCm39). The aggr pipeline was used to combine data from multiple samples into an experiment-wide feature-barcode matrix and analysis. The 10x Genomics Loupe Browser was used for visualization, initial quality assessment, and filtering of single cell gene expression data. Single Cell Gene Expression was performed at the Genomics Resource Laboratory, University of Massachusetts Amherst, MA.

Data analyses were performed using the Seurat package (version 4.3) in the R software version 4.2.1. The TCRαβ⁺ and TCRγδ⁺ datasets were analyzed individually but with identical procedures. The Seurat object was first generated by keeping all genes expressed by at least 3 cells. Cells were kept if they contained at least 100 unique features, at least 100 reads, and less than 10% mitochondrial genes. Data were normalized, scaled, then the top 2000 variable genes were used for the Principal Component Analysis (PCA) and generation of Uniform Manifold Approximation and Projection (UMAP) plots. Immune cells were filtered by expression of Ptprc > 0.5, then for T cells by keeping cells with at least 0.05% expression of Cd3e/Cd3g/Cd3d. Cluster identities were interpreted by calculating unique cluster marker genes and analyzing the distribution of known marker genes. Clusters with similar transcriptional profiles and localization on the UMAP plot were merged into groups. Poorly defined clusters were removed from further analysis.

Re-clustering of “CD8αα⁺” annotated cells from Figure 1 was performed by subsetting the “CD8αα⁺ Tcf7 ⁺”, “CD8αα⁺ Prdm1 ⁺”, and “Mki67 ⁺” clusters for the TCRαβ⁺ dataset, or the “CD8αα⁺ Tcf7 ⁺”, “CD8αα⁺ Prdm1 ⁺”, “CD8αα⁺ Zeb2 ⁺”, and “Mki67 ⁺” clusters for the TCRγδ⁺ dataset. Cells expressing Cd4 or Cd8b1 were removed, and the remaining cells were processed by the normalization and clustering steps as above. After examining the expression of Cd8a, cells were further separated into the CD8αα⁺ cells which expressed Cd8a, or DN cells which did not express Cd8a. Differentially expressed genes between clusters were calculated by the Wilcoxon ranked test. Pathway enrichment analyses were performed with Metascape, using the differentially expressed genes filtered by significance (adjusted p-value < 0.05) and Average Log₂ Fold Change > ± 0.5. The list of transcription factors in the mouse genome were obtained from the FANTOM5 database. The most differentially expressed transcription factors between clusters were selected to show as bubble plots ( Figures 2b, d , 3b, d ) by significance (adjusted p-value < 0.05) and proportion (difference in the percentage of cells expressing the gene > 20% for Figures 2b, d , Figure 3b , > 30% for Figure 3d ). RNA velocity analysis was performed using scVelo and Velocyto.

Integration of the TCRαβ⁺ and TCRγδ⁺ datasets was performed after preliminary analysis of individual datasets with the SelectIntegrationFeatures(), FindIntegrationAnchors(), and IntegrateData() Seurat functions. The integrated dataset was preprocessed following standard methods, and clusters were defined at low resolution (0.1) to determine which cell types colocalize on the UMAP plot.

Data statistical analysis was performed with Prism 9 (GraphPad software). P-values were determined using a two-tailed paired t-test, or one-way ANOVA.