Introduction

Like αβ T cells, γδ T cells can be broadly grouped into innate-like IFN-γ-producing γδT1, IL-4-producing γδT2, or IL-17-producing γδT17 subsets that are functionally programmed during thymic development^{1, 2, 3, 4}, or inducible subsets that differentiate in the periphery⁵. These innate-like γδ T cells are highly enriched in mucosal tissues such as the skin, gut, and lung and play critical roles in the host response to pathogens, tumor surveillance, autoimmunity, and tissue homeostasis. Accumulating data suggest that the relative balance of these functional γδ T cell subsets can have important implications in disease susceptibility and progression. This is perhaps best exemplified in the setting of cancer where IFN-γ-producing γδT1 have been shown to play a critical role in tumor elimination^{6, 7, 8}, while γδT17 subsets promote tumor growth and metastasis^{9, 10, 11, 12}. The thymic developmental programs, therefore, that regulate the balance of these functional subsets in the periphery may have considerable impact on disease outcome. Unlike the case for their αβ T cell counterparts, however, the mechanism that determine whether a developing γδ T cell acquires an innate-like γδT1 or γδT17 phenotype during thymic development remain unclear.

A significant body of data supports a model in which weak or strong TCR signal strength might contribute to the development of innate-like γδT17 (weak signal) and γδT1 (strong signal) subsets^{4, 13, 14, 15, 16}, and a strong TCR signal has been suggested to be necessary for the development of the IFN-γ and IL-4 producing γδNKT cell subset¹⁷. However, there is also evidence of a TCR-independent model in which pre-existing hard-wired transcriptional programs regulate the development of some innate-like γδΤ17 subsets¹⁸, and accumulating data suggest that γδT17 programming is associated with the use of distinct signaling pathways^{19, 20, 21}.

SLAM family receptors comprise a family of nine receptors, SLAMF1 (SLAM; CD150), SLAMF2 (CD48), SLAMF3 (Ly9; CD229), SLAMF4 (2B4; CD244), SLAMF5 (CD84), SLAMF6 (Ly108; CD352), SLAMF7 (CRACC; CD319), SLAMF8 (BLAME), and SLAMF9 (SF2001; CD84H), that are expressed primarily on hematopoietic cells and which play diverse roles in immune development and function²². Most SLAM family receptors interact in a homophilic manner, the exception to this rule being SLAMF4, which serves as the ligand for SLAMF2. SLAM family receptor signals can be both activating or inhibitory, depending on the recruitment of the cytosolic signaling adapter proteins SAP and EAT-2, and inhibitory SHP-1, SHP-2, or SHIP phosphatases to immunoreceptor tyrosine switch motifs in the SLAM family receptor cytoplasmic domain²². Although SLAM/SAP signaling has long been known to play a role in innate-like NKT²³ and γδNKT¹⁷ cell development, recent findings from multiple groups demonstrating that this pathway is also involved in the development of innate-like γδT1, γδT17²⁴ and innate-like MAIT cells^{25, 26} indicates that SLAM/SAP signaling plays a fundamental role in the development of a majority of the innate-like T cell lineages.

In this study, we set out to define the stage of thymic development during which SAP mediates its effect on γδ T cell developmental programming and to define SAP-dependent changes in the γδ T cells transcriptional program during their development. Finally, since SAP-dependent innate-like T cell subsets are often characterized by highly restricted TCRs (e.g., Vγ1⁺/Vδ6.3⁺ γδNKT²⁷), we wished to identify whether specific TCR clonotypes were also associated with SAP-dependent γδT1 and γδT17 T cells. Therefore, we utilized a single-cell proteogenomics approach with V(D)J profiling that allowed us to simultaneously compare cell surface protein expression, gene expression, and TCR clonotypes between B6 and B6.SAP^−/− γδ T cells from embryonic, neonatal, and adult thymus, as well as from peripheral lung γδ T cells. Collectively, these data provide a comprehensive map of the SAP-dependent γδ T cell subsets and their developmental checkpoints in perinatal and adult thymus.

Results

SLAM family receptor expression marks transcriptionally distinct developmental stages among E17 γδ T cells

γδT17 developmental programming is thought to occur primarily during embryonic/neonatal thymic development^{2, 18, 28}, and is tightly linked to the expression of discrete TCR gamma and delta pairings²⁹. To identify the specific developmental stages of thymic development during which SLAM/SAP signaling pathway exerts its effects, we employed a single-cell RNAseq approach coupled with feature barcoding (i.e., CITEseq³⁰) and V(D)J profiling to compare embryonic day 17 (E17) B6 and B6.SAP^−/− thymic γδ T cells (Fig. 1A). FACS-sorted thymic γδ T cells (Suppl. Fig. 1) from individual embryos (n = 4 B6 and 4 B6.SAP^−/− mice) were identified using hashtags, and Vg1⁺ and Vg4⁺ gd T cells were identified using barcoded antibodies (ADTs). To help define the developmental trajectories of the E17 gd T cells, we also stained cells with barcoded antibodies specific for CD24, CD73, CD44, and CD45RB, which define distinct developmental pathways of γδT17 and γδT1 cells in the thymus^{4, 31, 32}, as well as for SLAMF1 and SLAMF6. After data filtering and QC, we verified that clusters were evenly distributed among the biological replicates (Suppl. Fig. 1). The resulting dataset represented a total of 5278 thymic γδ T cells from B6 (n = 2597) and B6.SAP^−/− (n = 2681) γδ T cells. Data were visualized using uniform manifold approximation and projection (UMAP) for dimensionality reduction and a graph-based clustering approach³³.

Figure 1.

First, we annotated the B6 dataset using a combination of cell surface protein expression, gene expression, and TCR clonotype sequence. This analysis revealed 11 B6 E17 γδ T cell clusters, mostly distributed between two major branches (Fig. 1B). We observed the presence of CD24^low73^low44^high45RB^lowSLAMF1^pos (c8) and CD24^low73^high44^high45RB^highSLAMF6^low/pos (c0, c9) clusters, consistent with mature γδT17 phenotype and γδT1 phenotypes, respectively (Fig. 1C, Suppl. Fig. 1). In support of this annotation, the c8 cluster was enriched in Il17a, Il17f, Il23r, Rorc, and TRGV6/TRDV4 (IMGT nomenclature; encoding Vγ6/Vδ1) transcripts with limited CDR3 diversity, while the c0/c9 cluster was enriched in Ifng, Tbx21, Klrd1, Il2rb, Xcl1, and TRGV5/TRDV4 (encoding Vγ5/Vδ1) transcripts also with little to no diversity (Figs. 1D – F, Suppl. Fig. 2, Suppl. Tables 1 and 2)^{34, 35}. Consistent with previous observations²⁴, we noted that the c0/c9 cluster was highly enriched in Cd244a (Slamf4), and was enriched in transcripts from another SLAM family member, Slamf7 (Figs. 1D - E). Flow cytometric analysis confirmed that SLAMF7 was expressed primarily on mature γδ T cells and that its expression was mostly confined to γδ T cells expressing the canonical CD44⁺CD45RB⁺ markers of γδT1 cells (Suppl. Fig. 1)^{4, 36}.

We identified a small CD24^low/int73^high44^high45RB^pos c3 cluster characterized by enrichment in Cd200 (whose expression is upregulated upon TCR ligation³⁷), Il2rb, Nr4a1, Cd244a, Xcl1, and Tnfrsf9 (encoding 4-1bb receptor). A high frequency of invariant TRGV5/TRDV4 clonotypes and decreased TCR clonotype diversity (Suppl. Fig. 2, Suppl. Table 2) suggested that c3 represented an immature γδT1 population. The CD24^low/posCD73^low/pos44^low/posCD45RB^low/pos c4 cluster was enriched in Ms4a4b and Ms4a6b (encoding CD20 homologues that modulate T cell activation^{38, 39}), S1pr1, Il6ra, Cd8a, Cd8b1, and Ccr9, suggesting that the c4 cluster represents the previously described emigrating thymic Ccr9⁺ S1pr1⁺ γδ T cells⁴⁰ (Fig. 1C and 1E). Consistent with the c4 cluster representing naïve emigrating thymic γδ T cells, the diversity of the c4 TCR clonotypes was the highest among the CD24^low/int (c0/c9, c8, c3) clusters (Suppl. Fig. 2).

We identified several immature CD24^highCD73^neg ADT clusters (c1, c2, c7, c5, c6, and c10) that also expressed SLAMF1 and/or SLAMF6 (Fig. 1B and Suppl. Fig 1). The c1 and c2 clusters were enriched in γδT17-associated genes such as Blk, Maf, Sox13, Ccr2, Tcf12, and Rorc^{32, 41, 42, 43, 44}, as well as Sh2d1a (encoding SAP) and Ccr9¹⁶ (Fig. 1E, Suppl. Table 1). Within these two clusters, the CD44^low/posCD45RB^low c2 cluster was especially enriched in genes previously implicated in γδ T cell development such as Etv5, Tox, Sox4, Tcf7, Zbtb16 ^{44, 45, 46, 47, 48}, Igfbp4 and Igf1r which were recently implicated in Th17 differentiation⁴⁹, Tnfsf11 (RANKL) which regulates Aire⁺ mTEC development ⁵⁰, as well as Themis and Tox2 (Fig. 1E, Suppl. Fig. 1). A high level of Gzma expression identified c2 as the previously described Gzma⁺ γδ T cell subset⁴⁰. Interestingly, in addition to expressing γδT17-associated genes, the c2 cluster was also enriched in γδT2-associated Gata3 and Il4 transcripts (Fig. 1D – E, Suppl. Table 1). The CD44^low/posCD45RB^neg c1 cluster, in contrast, exhibited significantly higher levels of Il17a, Il17f, Rorc, Maf ⁴², as well as Ifngr1, consistent with the c1 cluster representing immature γδT17 cells (Fig. 1D – E, 1G, and Suppl. Fig. 1).

Analysis of the c1 and c2 TCR transcripts revealed that the c2 cluster was notably enriched in TRGV4-expressing cells and that 87% of Vγ4 T cells utilized just 3 TRDV chains (50% TRDV5, 24% TRDV2-2, and 12% TRAV13-4/DV7) (Fig. 1F, Suppl. Table 2). Clonotype analysis revealed a high frequency of semi-invariant TRGV4/TRDV5 and TRGV4/TRDV2-2 clonotypes previously associated with γδT17^{51, 52, 53} (Suppl. Fig. 2, Suppl. Table 2). In contrast, the c1 cluster possessed comparatively fewer of these TRGV4 clonotypes, and was enriched in TRGV6, TRGV5, and TRGV1 paired with TRDV4 (Suppl. Fig. 2, Suppl. Table 2). Indeed, we noted that the Vγ1 repertoire, which is diverse in adult tissues, was dominated by TRGV1/TRDV4 clonotypes (35% of all TRGV1) in E17 thymus, and just two closely related semi-invariant TRGV1/TRDV4 clonotypes comprised 17.4% of all TRGV1. The location of these clonotypes in the immature γδT17-related clusters (Suppl. Fig. 2) suggested that they represent the dominant IL-17-producing Vγ1 clonotypes at E17, which we verified using flow cytometry (Suppl. Fig. 3).

Since the IL-4-producing Vγ1/Vδ6.3 γδNKT cell subset develops in fetal thymus¹, we asked whether the Il4-enriched γδ T cells in cluster 2 (Fig. 1D) were immature Vγ1 γδNKT cells with previously identified canonical Vδ6.3 sequences²⁷. Interestingly, while we did identify a small number of Vγ1 T cells with these canonical sequences (encoded by TRAV15-1/DV6-1 or TRAV15N-1) in our dataset, these were not present in the Il4-enriched c2 cluster.

The CD44^posCD45RB^pos c5 cluster was enriched in Il2ra, Scin, Bex6, Cd5, Slamf6, Cd28, Hivep3, and cell-cycle/metabolism genes, consistent with γδ T cells that have recently developed from DN2/DN3^{40, 54, 55}. Within the c5 and c7 clusters, we noted that while there was expression of γδT17-associated genes such as Blk, Maf, Sox13, and Etv5, there was relatively little Rorc expression. Rorc expression became apparent in the c2 cluster (Rorc^low) and was significantly increased in the c1 cluster (Suppl. Fig. 1). Flow cytometric analysis confirmed BLK and PLZF expression in CD25^high E17 Vγ4 T cells, while RORγt expression was primarily observed in the CD25^neg population, suggesting that BLK is expressed soon after γδ T cell selection, followed by RORγt (Suppl. Fig. 3). These data were consistent with MAF-directed expression of Rorc^{40, 42} and suggested a c5 – c7 – c2 – c10 / c1 – c8 γδT17 developmental trajectory, which was supported by trajectory inference analysis (Fig. 1H).

Finally, analysis of the CD24^high73^neg44^neg45RB^intSLAMF1^posSLAMF6^pos c6 population revealed that it was highly enriched in Cd8a, Cd8b1, Cd4, Rag1, Rag2, Ptcra, Arpp21, and Dgkeos, consistent with cells developing along the αβ T cell lineage (Figs. 1E, 1G, Suppl. Fig. 1, and Suppl. Table 1)⁵⁶. Analysis of the Vγ1 and Vγ4 ADTs confirmed that c6 cells were not simply contaminating αβ T cells as both Vγ1- and Vγ4-expressing cells were present in the c6 cluster (Suppl. Fig. 3) and analysis of the c6 TCR clonotypes revealed the presence of TRGV1, GV4, GV5, GV6, and GV7 clonotypes, none of which appeared to be specific to the c6 cluster (Suppl. Table 2).

Since trajectory inference analysis suggested a strong relationship between the Il2ra⁺ c5 and c6 clusters (Fig. 1H), we conducted differential gene expression analysis between the c5 and c6 from B6 mice. Interestingly, we found that the c6 cluster was highly enriched in Rorc, suggestive of a γδT17 phenotype, but was mostly lacking in the expression of Maf, Blk, Ccr2, and other canonical γδT17-associated transcripts (Fig. 1E, 1G, and Suppl. Fig. 1). Indeed, flow cytometric analysis of E17 TCRβ^negTCRδ^pos γδ T cells confirmed that CD4CD8 DP γδ T cells exhibited an immature CD44^lowCD25^lowMAF^negBLK^negRORγt^pos phenotype (Suppl. Fig. 3), and that RORγt expression in the absence of BLK or MAF is a characteristic of Lin^neg DN thymocytes as they progress from CD44^low25^high DN3 to CD44^low25^low DN4 (Suppl. Fig. 3). Together, these data suggested that the c6 cluster γδ T cells represented previously described CD4CD8 DP γδ T cells that were developing along an alternative αβ T cell pathway^{57, 58, 59}.

SAP regulates γδ T cell clonotype diversion into both the γδT17 and the αβ T cell developmental pathways

To identify SAP-dependent developmental checkpoints, we next compared the hashtag-separated B6 and B6.SAP^−/− E17 γδ T cell datasets. A comparison of the cluster frequencies between B6 and B6.SAP^−/− γδ T cells revealed significantly decreased frequencies of the CD24^high73^neg44^int45RB^intSLAMF1^posSLAMF6^pos c2 and c1 clusters in B6.SAP^−/− mice, and a striking increase in the frequency of the CD24^high73^neg44^neg45RB^negSLAMF1^highSLAMF6^high c6 cluster (Fig. 2A). A comparison of differentially expressed genes among the immature (c1, c2, c7, c5, c6, and c10) clusters using pseudobulk analysis revealed 234 differentially expressed genes with a p_adj threshold < 0.005 and log₂ fold change greater than 0.5, with 148 genes exhibiting decreased expression in B6.SAP^−/− γδ T cells, and 86 genes exhibiting increased expression (Fig. 2B, Suppl. Table 3). Consistent with the decreased frequencies of the c1/c2 clusters in B6.SAP^−/− mice, a significant fraction of the differentially expressed genes that were down-regulated in B6.SAP^−/− mice reflected γδT17-associated genes (e.g., Maf, Blk, Ccr2, Sox13, Etv5, etc.) that were present in the c2 and c1 clusters, as well as a number of c1/c2-enriched genes (e.g., Cpa3, Tox2, Pdcd1, Tmem121, Bex6, and Scin) whose importance in γδ T cell development is unclear (Fig. 2B-C, Suppl. Table 2). This was not solely due to the decreased frequency of the c2 and c1 clusters in B6.SAP^−/− mice, as a closer analysis revealed a decreased number of transcripts for many signature genes (e.g., Blk, Etv5, Ccr2) prevalent in the c2 cluster and to a lesser extent the c7 cluster (Suppl. Fig. 4). Given the decreased expression of Blk, we also examined the expression of Src tyrosine kinase family members Lck and Fyn, which have previously been implicated in T cell development⁶⁰. Unlike Blk, which was mostly restricted to γδT17 clusters, Lck and Fyn expression was more broadly distributed (Suppl. Fig. 1). Using pseudobulk analysis to compare their expression between B6 and B6.SAP^−/− immature E.17 γδ T cells, we found that neither Lck nor Fyn expression was appreciably altered (Suppl. Fig. 4), suggesting that Blk expression was uniquely dependent on SAP. We did note, however, that the magnitude of Lck differential expression was close to the 0.5 log₂ FC cut-off and that its expression was increased in B6.SAP^−/− γδ T cells (Suppl. Fig. 4). Indeed, a large number of genes whose expression was increased in B6.SAP^−/− mice were those whose expression was mostly confined to the c6 cluster (e.g., Cd4, Cd8b1, Ptcra, Rag2) or whose expression was highest in c6 (e.g., Tcf7, Thy1, Lck) (Fig. 2B-C, Suppl. Fig. 4). Together, these data identified the Gzma⁺Etv5⁺Sox13⁺Blk⁺Tox2⁺c2 and Cd4⁺Cd8b1⁺Ptcra⁺Rag2⁺c6 clusters as SAP-dependent developmental checkpoints.

Figure. 2.

Next, we asked whether these SAP-dependent alterations in γδ T cell development were associated with any changes in the TCR repertoire. A comparison of all γδ TCR clonotypes between E17 B6 and B6.SAP^−/− γδ T cells revealed no significant changes in overall clonotype frequency (Fig. 2D). However, a closer examination of TCR repertoire revealed significant changes in the distribution of specific TCR clonotypes. For example, within the TRGV4 repertoire, which represents the dominant TCR among immature clusters at E17, we observed a significant decrease in TRGV4/TRDV5 clonotypes in the c2 cluster and a concomitant increase of these clonotypes in the c6 cluster (Fig. 2E). Among the TRGV4 clonotypes, this appeared to be specific to the TRDV5 clonotypes, as we observed no significant differences c2 and c6 frequencies among TRGV4/TRDV2-2 and TRDV7 clonotypes (Fig. 2E). We did note, however, a significant increase in the TRGV4/TRDV7 clonotypes in the Il2ra⁺ c5 cluster (Fig. 2E). Within the TRGV1 repertoire, we noted a significant increase in TRGV1/TRDV4 clonotypes within the c6 cluster (Fig. 2E). Together, these data suggested that SAP operated at an early stage of γδ T cell development during or soon after the down-regulation of Il2ra, and that it positively regulated entry into the γδT17 pathway at the c5/c7 to c2 transition, while at the same time inhibiting entry into the CD4⁺CD8⁺ αβ T cell-like c6 cluster. The net result of these changes appeared to be that specific TCR clonotypes such as TRGV4/TRDV5 were diverted from the γδT17 pathway into the αβ T cell-like pathway.

SAP deficiency is associated with increased numbers of immature CD4⁺CD8⁺RORγt⁺ thymic γδ T cells

Next, we independently evaluated these observations using flow cytometry. This analysis revealed that while SAP deficiency did not affect the overall number of E17 γδ T cells (Suppl. Fig. 5), it did result in a significant decrease of immature CD24^posBLK^posMAF^posRORγt^pos γδ T cells in E17 B6.SAP^−/− thymus (Fig. 3A-B). We found that this decrease was most pronounced in the CD44^pos population after downregulation of CD25 (Fig 3A, Suppl. Fig. 5) and that this population exhibited a high level of PLZF expression, suggesting it corresponded to the CD44⁺Blk⁺Maf⁺Zbtb16⁺Rorc⁺ c1 cluster (Suppl. Fig. 5). In addition, UMAP analysis of E17 Vγ4 T cells flow cytometric data revealed the presence of a SAP-dependent CD24^posBLK^posPLZF^posRORγt^neg population in B6.SAP^−/− thymus (Suppl. Fig. 5) consistent with the phenotype of the Rorc^low c2 cluster. In contrast, we found no evidence of a SAP-dependent decrease in the CD25^pos population frequency, consistent with the notion that SAP exerts its effect during or after this stage of development (Fig. 3A). We also noted that in addition to its effect on the numbers of immature thymic γδT17, SAP deficiency resulted in the reduced expression levels of some of these markers, most notably BLK (Fig. 3C), a finding that was consistent with the decreased number of Blk transcripts observed in Fig. 2. Interestingly, as PLZF induction is often associated with SLAM/SAP signaling, we noted that SAP did not affect overall PLZF expression levels in E17 thymic γδ T cells (Suppl. Fig. 5).

Figure 3.

In contrast to the decreased numbers of CD24^posBLK^posMAF^posPLZF^pos RORγt^pos γδ T cells, we observed a significant increase in the frequency and number of CD24^posBLK^negMAF^negPLZF^negRORγt^pos γδ T cells (Figs. 3A-B), whose phenotype was consistent with the αβ T cell-like c6 cluster in Fig. 2. Consistent with our finding that the c6 cluster was enriched in CD4 and CD8, we found that the frequency and number of CD4 and CD8-expressing CD24^pos TCRδ^posTCRβ^neg BLK^negMAF^negRORγt^pos thymic γδ T cells was significantly increased in B6.SAP^−/− mice (Fig. 3D). Together, these analyses confirmed our single-cell CITEseq analysis, and suggested that SAP regulated both the throughput of immature γδ T cells in the γδT17 cell pathway as well as the αβ T cell-like pathway. More specifically, the data suggested that SAP promoted progression to the CD24^pos BLK^posMAF^posPLZF^posRORγt^low c2 cluster that is enriched in a number of genes (e.g., Sox13, Etv5, Blk, Maf) known to regulate γδT17 development, but that it inhibited entry into the CD24^pos44^negBLK^negMAF^negRORγt^pos αβ−T cell-like pathway (c6 cluster).

Identification of SAP-dependent developmental checkpoints during neonatal and adult thymus γδ T cell developmental programming

The neonatal period is a time of rapid expansion of numerous γδ T cell subsets in the thymus including Vγ4 γδT17 and γδNKT. Therefore, we conducted a comparison between B6 and B6.SAP^−/− neonatal and adult thymic γδ T cells, to identify novel SAP-dependent developmental checkpoints and to track the outcome of the SAP-dependent alterations in γδ T cell development observed in E17 thymus. Neonate datasets revealed a similar branched clustering pattern as in E17 with some notable differences. Similar to E17, we observed immature CD24^pos73^lowSLAMF1^posSLAMF6^pos clusters (c0, c1, c2, c3, c9, c12, c14), CD24^pos73^low/posSLAMF1^posSLAMF6^pos (c4, c5, c10, c15) S1pr1-enriched “naïve” clusters, mature CD24^neg73^lowSLAMF1^posSLAMF6^neg (c6, c7) gdT17 clusters, and a mature CD24^neg73^posSLAMF1^negSLAMF6^low/^posSlamf7^high (c8) gdT1 cluster (Fig. 4A-B, Suppl. Fig. 6). In general, we noted good correlation with cluster-specific gene expression between the E17 and D9 (Day 9) datasets. A notable exception was a striking decrease in Zbtb16 expression among all immature neonatal γδ T cells, which is consistent with previous reports^{16, 40} (Suppl. Fig. 6, Suppl. Table 3).

Figure 4.

In contrast to the E17 dataset, the neonatal thymic γδ T cell dataset contained a CD24^neg73^posSLAMF1^negSLAMF6^highSlamf7^neg (c13) γδT2 cluster enriched in Il4, Il13, Zbtb16, Cd40lig, and Nos1 (Suppl. Fig. 6, Suppl. Table 1), suggesting it represented IL-4-producing γδNKT cells. Indeed, examination of the γδT2 (c13) and γδT1 (c8) cluster TCR repertoires revealed a significant enrichment of canonical TRGV1 and TRAV15N-1 or TRAV15-1/DV6-1 sequences characteristic of γδNKT cells (Fig. 4C, Fig. 5)²⁷. Analysis of the differentially expressed genes in the γδT1 (c8) and γδT2 (c13) clusters revealed that Zbtb16, Icos, and Slamf6 were markers of γδT2 cells while Tbx21, Klrb1c (NK1.1), and Slamf7 were markers of γδT1 cells in neonatal thymus (Suppl. Fig. 6). Flow cytometric analysis confirmed that PLZF^hi γδT2 Vγ1Vδ6.3 cells co-express SLAMF6 and ICOS while Tbet^hi γδT1 Vγ1Vδ6.3 cells are NK1.1⁺ and are either SLAMF6⁻SLAMF7⁺ or SLAMF6⁺SLAMF7⁺ (Suppl. Fig. 7).

Figure 5.

A comparison of B6 and B6.SAP^−/− neonatal gd T cells revealed multiple SAP-dependent changes in clustering. Among the immature γδ T cells, we observed a decreased frequency of the immature CD24^high73^negSLAMF1^posSLAMF6^pos c1 cluster (Fig. 4D) that was enriched in Blk, Etv5, Sox13, Rorc, and Maf (Suppl. Fig. 6), suggesting it was analogous to the SAP-dependent c2 cluster in E17 γδ T cells. Similar to our observations in E17 γδ T cells, flow cytometric analysis confirmed both a decrease in BLK^posRORγt^pos and BLK^posRORγt^neg γδ T cell populations, as well as a decreased level of BLK expression (Fig. 4E-F). In addition, we observed an increased frequency of immature CD24^high73^negSLAMF1^posSLAMF6^pos cluster with characteristics of αβ T cells (c14) in B6.SAP^−/− D9 thymic γδ T cells (Fig. 4D, Suppl. Fig. 6). While this increase did not reach the level of significance, presumably due to the very low number of c14 cells in the dataset, flow cytometric analysis revealed a significant increase in immature CD24^pos73^negBLK^negRORgt^pos gd T cells, a phenotype consistent with gd T cells that have been diverted to the ab developmental pathway (Fig. 4E). Last, we observed an increased frequency of immature CD24^pos73^low/pos SLAMF1^posSLAMF6^pos S1pr1-enriched c4 and c5 clusters, but not the closely related Cd200-enriched c10 cluster (Fig. 4D). While we did confirm an increased number of CD200^negS1pr1^pos γδ T cells in neonatal thymus by flow cytometry, this increase was relatively small and its significance is unclear (Suppl. Fig. 7). Independent confirmation of these SAP-dependent changes in neonatal γδ T cell gene expression was obtained by conducting bulk RNAseq analysis of FACSs-sorted CD24^high neonatal thymic γδ T cells from B6 and B6.SAP^−/− mice (Suppl. Table 4).

Among the mature neonatal γδ T cells, we observed a decreased frequency of the CD24^low73^negSLAMF1^posF6^low γδT17 c6 cluster, which was consistent with our previous observation of decreased IL-17 production in B6.SAP^−/− neonatal γδ T cells²⁴. Interestingly, we observed a striking loss of the CD24^low73^highSLAMF1^lowF6^highF7^neg γδT2 c13 cluster in B6.SAP^−/− mice, as well as a decreased frequency of the mature CD24^low73^highSLAMF1^lowF6^lowF7^high γδT1 c8 cluster (Fig. 4A - D). These findings were consistent with a SAP-dependent loss of both the IL-4- and IFN-γ-producing γδNKT as well as non-γδNKT γδT1 cells (discussed further below). Unfortunately, while we were able to define mature SAP-dependent γδNKT subsets, we found no evidence of immature γδNKT cells in our dataset, suggesting that they develop during a brief developmental window between E17 and D9 after birth.

A comparison of B6 adult thymic γδ T cells between B6 and B6.SAP^−/− mice recapitulated many of our observations of the neonatal thymus, although to a somewhat lesser degree. For example, we observed a small but significant decrease in the frequency of the CD24^highCD73^low adult thymic c1 and c2 clusters in B6.SAP^−/− mice that was coupled with an increased frequency of the CD24^posCD73^+/− c0 cluster (Suppl. Fig. 8). Like the SAP-dependent E17 c2 and D9 c1 clusters, the adult thymus c1 and c2 clusters were enriched in γδT17-associated Blk, Etv5, Sox13, Maf, and Rorc (mostly in c2) indicating that SAP regulates this γδ T cell thymic developmental checkpoint at all ages. Likewise, both the adult thymus c0 and the neonatal c4 and c5 clusters that were enriched in Ms4a4b, Ms4a6b, Dgka, Klf2, and S1pr1 exhibited increased frequencies in B6.SAP^−/− mice. Just as we observed in neonatal thymic γδ T cells, we noted a near complete loss of mature γδNKT cells (Zbtb16, Icos, and Il4-enriched γδT2 c8 cells and Tbx21, Slamf7, and Il2rb-enriched γδT1 c11 cells) as well as a decreased frequency of mature non-γδNKT γδT1 c10 cluster cells in B6.SAP^−/− adult thymus (Suppl. Fig. 8). Unlike the E17 and neonatal thymic datasets, however, the frequency of γδ T cells with an αβ T cell-like signature in the adult thymic γδ T cell dataset was minimal. Collectively, these data supported the presence of multiple SAP-dependent developmental checkpoints during neonatal and adult γδ T cell development—a reduction in the CD24^posCD73^neg SLAMF1⁺SLAMF6⁺ cluster enriched in γδT17 development-associated transcripts, an increased frequency of immature S1pr1-enriched γδ T cells, and decreased frequencies of both γδNKT cells (γδT1 and γδT2) as well as non-γδNKT γδT1 cells.

SAP shapes the thymic γδ TCR repertoire

Given our observation of altered γδ TCR clonotype distribution among E17 γδ T cells and the substantial changes in clustering we observed in SAP-deficient mice, we were next interested in determining the effect of SAP deficiency on the neonatal and adult γδ TCR repertoires. A comparison of frequencies of the top 100 TCR clonotypes between B6 and B6.SAP^−/− mice revealed significant changes in the overall frequency of several TRGV1 and TRGV4 clonotypes (Fig. 5A). Consistent with a SAP-dependent loss of γδNKT, we observed a profound loss of TRGV1 clonotypes utilizing the canonical TRAV15N-1 and TRAV15-1/TRDV6-1 γδNKT cell TCRs (Fig. 5B). Loss of these clonotypes was directly tied to the significant decrease of the γδT1 c8 and γδT2 c13 clusters, in which nearly all these clonotypes resided.

Within the B6 γδT2 c13 cluster the vast majority (∼86%) of TCR sequences were highly restricted canonical γδNKT cell clonotypes (Suppl. Table 2) and nearly all of these clonotypes were also present in the γδT1 c8 cluster consistent with the ability of γδNKT cells to produce both IFN-γ and IL-4. In contrast, while the B6 γδT1 c8 cluster contained a large proportion of canonical γδNKT cell clonotypes (∼46% of c8 TCRs), it also contained TCRs with a diverse pool of clonotypes that utilized TRGV1 (∼70% of c8 TCRs), TRGV4, TRGV5, TRGV6, and TRGV7 that were present only in the c8 γδT1 cluster. Therefore, in neonates, the thymic γδT1 population is a mixture of γδNKT and “non-γδNKT” γδT1 clonotypes, and each of these clonotype pools is SAP-dependent. As noted above, these findings were largely reproduced in the adult thymus where the γδNKT clonotypes associated with the c8 and c11 clusters were virtually absent in B6.SAP^−/− mice, and where the diverse pool of γδT1 c10 clonotypes was significantly diminished (Suppl. Table 2).

Interestingly, a closer examination of the neonatal thymic TRGV4 clonotypes revealed a significant decrease in the frequencies of γδT17-associated clonotypes with a germline-encoded TRDV5⁵² (Fig. 5B), some of which were identical to those that were redistributed to the αβ T cell-like c6 cluster in E17 γδ T cells (Fig. 2E). When these clonotypes were mapped to their D9 clusters, we noted a significant decrease in the frequency of TRGV4/TRDV5 clonotypes in the mature γδT17 c6 cluster (Fig 5C). We note that while we observed a trend toward an increased frequency of TRGV4/TRDV5 and TRGV4/TRDV2 clonotypes in the S1pr1⁺ c4/c5 clusters, analysis of these sequences revealed that they largely exhibited diverse CDR3 rather than the restricted CDR3 observed in the γδT17 c6 cluster (Suppl. Table 2). Taken together, these data indicated that the immature embryonic day 17 γδ TCR clonotypes that were redirected to the αβ T cell pathway in B6.SAP^−/− mice were underrepresented among the mature B6.SAP^−/− neonatal γδ T cells.

Restricted TCR repertoire in innate-like γδT1 as well as γδT17 in the periphery

We previously demonstrated that SAP-deficiency results in a significant impairment in γδ T cell function in the periphery²⁴. Specifically, we observed that γδ T cell IL-17 production is impaired in the lungs of young mice, but that this had largely corrected itself as the mouse reached maturity. In contrast, we observed that both Vγ1 and Vγ4 IFN-γ production were significantly impaired in B6.SAP^−/− adult lung and spleen. To better understand how SAP-dependent alterations in thymic development ultimately affected the peripheral γδ T cell compartment, we compared B6 and B6.SAP^−/− lungs using single-cell RNAseq with V(D)J profiling.

Our initial annotation of the B6 lung γδ T cell transcriptome revealed the presence of 11 clusters of the lung γδ T cells that we annotated based on cluster-specific gene expression and TRGV and TRDV usage (Figure 6A). Based on their expression of signature cytokines and transcription factors, we identified 2 γδT17 clusters (c1 and c4), 4 γδT1 clusters (c2, c3, c7, and c9), 1 γδTNKT-like cluster (c10) with very few cells, as well as a Cd24-enriched cluster (c5) (Fig. 6B and 6C, Suppl. Table 5). Consistent with previous observations, we noted that the c1 and c4 γδT17 clusters were enriched in Slamf1 and that the c2/c3/c7/c9 γδT1 clusters were enriched in Slamf7 and Slamf6 (Fig. 6B). Indeed, flow cytometric analysis revealed that SLAMF7 expression was restricted to the innate-like lung CD44⁺CD45RB⁺ γδT1 cell population, and that γδT1 cells are SLAMF1⁻F6⁺F7⁺, while the CD44⁺CD45RB⁻ γδT17 subsets are SLAMF1⁺F6⁻F7⁻ (Fig. 6D). Further annotation using V(D)J profiling revealed that the c1 cluster corresponded to Vγ4 γδT17 (primarily TRGV4/TRDV5 and TRGV4/TRDV2-2) while the c4 cluster corresponded to lung Vγ6 (TRGV6/TRDV4) (Fig. 6E, Suppl. Table 6). The c10 cluster was enriched in Zbtb16, Il4, and Icos (Fig. 6B-C) as well as TRGV1/TRAV15N-1 and TRGV1/TRDV6 sequences (Fig. 6E, Suppl. Table 6) consistent with γδNKT cells. We did note, however, that the TRAV15N-1 and TRDV6 sequences in the c10 cluster were, for the most part, non-canonical γδNKT cell clonotypes, suggesting a degree of flexibility in the TRAV15N-1 and TRDV6 CDR3 sequences that can be utilized by γδNKT cells.

Figure 6.

Interestingly, while the γδT1 c2/c3/c7/c9 clusters contained primarily Vγ1 and Vγ4 T cells as expected, we noted that Vγ4 γδT1 cells exhibited a striking preference in the usage of TRAV13-4(DV7), and that TRGV4/TRDV7 pairings accounted for 80% of all Vγ4 clonotypes in the γδT1 clusters (Fig. 6E, Suppl. Table 6). Analysis of the γδT1 TRAV13-4/DV7 CDR3 sequences revealed a high level of diversity that showed no evidence of clonal expansions (Fig. 6F). This was in contrast to the highly restricted Vγ6 repertoire in the c4 γδT17 cluster, as well as the restricted Vγ4 repertoire in the c2 γδT17 cluster that was dominated by semi-invariant TRDV5 and TRDV2-2 sequences (Fig. 6F). As no Vδ7-specific antibody was available, we independently confirmed this observation using single-cell PCR to assess the lung Vγ4 CDR3δ sequences from 13 individual mice. These data demonstrated that, other than TRDV2-2 and TRDV5, TRDV7 is the only other TCR delta chain used at an appreciable frequency by Vγ4 T cells in the B6 mouse lung and that it exhibits significant diversity in its CDR3 (Suppl. Fig. 9). Vγ1 γδT1 cells, in contrast, exhibited a more mild pairing preference where 44% of the γδT1 clonotypes utilized TRAV15N-1, TRAV15-1/DV6-1, or TRAV15D-2/DV6D-2 and 13% utilized TRAV13-4/DV7 (Suppl. Table 6). Together these data formed a comprehensive map of B6 lung γδ T cell transcriptomes and revealed that lung innate-like γδT1, like γδT17 and γδNKT, are all characterized by preferential TCR γ- and δ-chain pairings. This was especially notable for Vγ4 and our data suggested that Vγ4 γδT1 exhibited a distinct preference for the utilization of TRAV13-4(DV7) chains with diverse CDR3.

SAP regulates peripheral Vγ4/Vδ7 γδT1 cells and shapes the lung γδT17 repertoire

Next, we assessed the effect of SAP deficiency on the lung γδ TCR repertoire by comparing B6 and B6.SAP^−/− lung γδ T cells (Fig. 7A). Consistent with our previous observation that lung γδ T cell IFN-γ production is deficient in SAP-deficient mice²⁴, we observed a decrease in the frequency of the γδT1 c2, c3, and c7 clusters (Fig. 7B). Moreover, we observed a striking decrease in the frequency of TRGV4/TRDV7-expressing cells in the lungs of B6.SAP^−/− mice (Figs. 7C-D, Suppl. Table 6), and we noted that we found no evidence of a compensatory increase in the frequency of other TRDV-expressing lung Vγ4 γδT1 cells in the B6.SAP^−/− lung. We independently confirmed this SAP-dependent decrease in TRAV13-4/TRDV7 in both lung and spleen using qPCR (Suppl. Fig. 10), and we demonstrated a significant decrease in the number of CD44⁺SLAMF7⁺ Vγ4 T cells in both the lungs and spleen of B6.SAP^−/− mice (Fig. 7E, Suppl. Fig. 10). The decreased lung γδT1 did not appear to be due to a SAP-dependent decrease in homeostatic proliferation, nor to increased apoptosis, since we observed no significant change in lung γδ T cell BrdU incorporation or Annexin V staining in B6.SAP^−/− mice (Suppl. Fig. 10).

Figure 7.

In addition to the decrease in lung γδT1, we also noted a loss of the c10 γδT2 cluster consistent with the enrichment in SAP-dependent γδNKT TCRs in this cluster, and we observed an increased frequency of the Vγ4 γδT17 c1 cluster (Fig. 7B). Given our previous observation of a SAP-dependent decrease in lung γδ T cell IL-17 production in young, but not adult mice²⁴, and our findings here of SAP-dependent alterations in γδT17 clonotypes, we were interested in determining whether the adult lung γδT17 repertoire was altered in SAP-deficient mice. Indeed, a comparison of the lung Vγ4 c1 γδT17 CDR3 sequences between B6 and B6.SAP^−/− mice revealed significant alterations in the γδT17 c1 cluster TCR repertoire. Specifically, we noted a decreased frequency of TRGV4 chains that exhibited a germline-encoded CDR3γ (CSYGYSS) and an increased frequency of TRGV4 chains with a CDR3γ exhibiting evidence of N/P additions (CSYG(X)YSS). These changes were associated with a significant shift in the usage of TRDV2-2 CDR3δ sequence usage (Fig. 7F) and an altered c1 TCR clonotype frequency in B6.SAP^−/− mice (Fig. 7G). Although the number of germline-encoded invariant TRGV4/TRDV5 clone (clonotype 1 in Fig. 6B) was diminished in B6.SAP^−/− lung γδ T cells (Fig. 7G, Suppl. Table 6), the low numbers of cells available for analysis precluded us from making strong conclusions regarding the loss of any one particular clonotype. Taken together, these data suggest that Vγ4 γδT1 preferentially utilize the TRAV13-4(DV7) chain in the lung and that the SAP-dependent decrease in peripheral γδ T cell IFN-γ production is not due to dysregulated cytokine production, but rather due to a specific loss of innate-like γδT1 subsets. This finding, together with our finding that the adult lung γδT17 TCR repertoire is altered in SAP-deficient mice suggests that SLAM/SAP signaling plays a role in shaping the γδ TCR repertoire.

Discussion

A common denominator in most innate-like T cell developmental programming is a dependence on the SLAM/SAP signaling pathway^{17, 24, 25, 61}. Within the γδ T cell lineage, SAP plays a critical role in the development of IFN-γ and IL-4-producing γδNKT cells²⁷, as well as in γδT1 and γδT17²⁴, but exactly when during γδ T cell development SAP functions, and how it exerts these effects is unknown. Here, we provide evidence that SAP mediates its effects on immature, uncommitted thymic γδ T cells and that it works to simultaneously promote progression through the γδT17 developmental pathway while inhibiting the diversion of γδ T cells into the αβ T cell pathway. The redirection of γδ T cell clonotypes into the αβ T cell pathway at E17 is associated with a decreased frequency of these mature clonotypes in neonatal thymus, and an altered γδT17 TCR repertoire in the periphery. Finally, we identify Vγ4/Vδ7 T cells as a novel, SAP-dependent IFN-γ-producing Vγ4 γδT1 subset.

The developmental pathway(s) that gives rise to γδ T cell functional subsets remains unclear. γδ T cells develop at the DN2/DN3 stage after receiving a sufficiently strong γδTCR-mediated signal⁶². Numerous reports support a model where, within the continuum of signals that can drive γδ T cell development, it is the relatively weaker TCR signals that promote γδT17 programming while the stronger TCR signals foster γδNKT and γδT1 programming^{4, 14, 17, 20, 63, 64}. Despite considerable effort, the precise thymic developmental stages at which γδ T cell lineage diversification occurs remains unclear. Recent studies examining this issue have demonstrated that effector programming can occur at a CD44^negCD45RB^neg stage^{4, 20}, and that c-Maf-directed γδT17 programming occurs at a CD45RB^low stage⁴². In addition, by altering TCR signal strength in a transgenic model, Chen et al., defined a Ccr9-enriched cluster as a point of effector lineage diversification¹⁶. Our findings here suggest that SAP influences γδT17 developmental programming during the immature uncommitted CD24^hiCD73^neg stage, either during or soon after Il2ra (CD25) downregulation. The SAP-dependent c2 cluster we define here was highly enriched in numerous genes (e.g., Sox13⁴⁶, Blk⁴¹, Maf⁴², Etv5⁴⁵, Sh2d1a²⁴) previously associated with γδT17 development as well as Ccr9, and the cell surface phenotype of c2 γδ T cells is SLAMF1^posSLAMF6^posCD44^low/highCD45RB^low. Altogether, these data suggest that the SAP-dependent c2 cluster is identical to these previously defined developmental checkpoints. Given that these previous studies have demonstrated that this developmental checkpoint is influenced by TCR signaling, it is interesting to note that SLAM/SAP signaling has previously been demonstrated to foster semi-invariant NKT cell differentiation through its ability to regulate TCR signal strength⁶⁵.

To that end, we noted that the SAP-dependent c2 cluster in E.17 γδ T cells was enriched in Blk, Lck, and Fyn, Src tyrosine kinase family members that have all been implicated in γδ T cell development^{41, 60}. Decreased expression of Blk, which is known to promote γδT17 development⁴¹, in B6.SAP^−/− γδ T cells suggests one possible mechanism through which SAP may regulate γδT17 development via TCR signaling. The mechanism, however, through which SAP regulates Blk gene expression remains unclear. Since SAP is known to regulate conventional αβ T cell signaling through its ability to recruit FYN^{66, 67} and LCK^{68, 69} to SLAM family members, it seems reasonable to hypothesize that SAP could regulate TCR signaling by recruiting one or more of these Src kinase family members to the γδ TCR signalosome.

Our finding that SAP inhibited the number of γδ T cells that enter the αβ T cell developmental pathway is reminiscent of previous work using transgenic TCR models demonstrating that a reduction of TCR signal strength results in a redirection of γδ T cells to the αβ T cell lineage^{59, 70}. Fahl et al. recently demonstrated that both entry into the γδ T cell lineage as well as γδ T cell effector programming involved repression of TCF1 expression by TCR-driven increases in Id3⁴⁸. Consistent with these results, we found that Tcf7 (encoding TCF1) transcript levels were highest in the αβ T cell-like c6 cluster, and we note that Tcf7 transcripts levels were noticeably higher in the Il2ra-enriched post-selection c5 cluster in B6.SAP^−/− γδ T cells. Id3 expression, however, did not appear to be appreciably altered in SAP-deficient mice. Whether SAP-mediated signals could alter γδ T cell Tcf7 expression through a TCR-independent mechanism remains an open question.

In addition to the considerable evidence supporting a role for TCR signal strength in γδ T cell developmental programming, there is increasing support for a model in which γδT17 programming occurs earlier in development, prior to expression of the γδ TCR. Indeed, evidence that some DN1 thymocytes possess a SOX13^hi IL-17-like phenotype and that Vγ4 γδT17 arise from DN1 precursors^{18, 71} indicates the existence of such a TCR-independent developmental programming stage. In consideration of such a model, it may be possible that SLAM/SAP-mediated signals could differentially regulate the survival and/or fitness of the functionally pre-programmed γδ T cell subsets. Such a function would be consistent with the role of SLAMF6 and SAP in promoting reactivation-induced cell death in conventional αβ T cells^{69, 72}.

While our data indicate that SAP regulates the diversion of at least some TCR clonotypes into the αβ TCR pathway, the fate of those cells is unclear. Certainly, the fact that the clonotypes that were diverted to the αβ T cell pathway in E17 thymus are under-represented in the mature D9 thymus γδT17 clusters suggests that this diversion has significant consequences on the γδ T cell receptor repertoire, but whether these diverted γδ T cells contribute to the αβ T cell compartment is unclear. Recently, a population of hybrid αβ-γδ T cells was described in both mouse and human⁷³ that responds to peptide/MHC stimulation like conventional αβ T cells, as well as to IL-1 and IL-23 which induce both IFN-γ and IL-17 production⁷³. In the mouse, at least some of these cells are Vγ4/Vδ5, Vγ4/Vδ4, or Vγ4Vδ7 and are present in the embryonic thymus, raising the possibility that γδ T cells diverted to the αβ T cell pathway in B6.SAP^−/− mice contribute to this αβ-γδ T cell population. While we have not observed a significant increase in the overall numbers of αβ-γδ T cells in B6.SAP^−/− thymus (not shown), a direct comparison of B6 and B6.SAP^−/− αβ-γδ TCRs may be necessary to see such an effect.

Our previous observations that SAP deletion did not completely abolish the development of IL-17- and IFN-γ producing γδ T cells²⁴ suggested the possibility that SAP may be critical for the development of additional γδ T cell clonotypes other than the Vγ1Vδ6.3 subset. Therefore, we compared the SAP-dependent and SAP-independent γδTCR repertoires using paired TCR clonotype analysis, which allowed us to track clonotype fate during development. Our findings not only confirmed previous results demonstrating highly restricted γδT17^{51, 52}, and γδNKT²⁷ TCR repertoires, they also revealed that the vast majority of Vγ4 γδT1 subsets in the periphery preferentially utilize TRGV4/TRAV13-4(DV7) clonotypes, thereby further highlighting the strong links between γδ T cell function and the utilization of specific TCRs²⁹. The reason for this preferential TRDV usage in Vγ4/Vδ7 γδT1 is unclear. While the TRAV13-4(DV7) CDR3 sequences were diverse and did not suggest evidence of clonal expansion, it may be possible that TRDV-encoded sequence motifs could be required to interact with a ligand, similar to the way that TRGV-encoded sequences mediate recognition of butyrophilins^{74, 75}.

Exactly when during development SAP exerts its effects on Vγ4/Vδ7 T cells is still somewhat unclear. We found comparatively low numbers of TRGV4/TRDV7 γδ T cells in embryonic, neonatal, or adult thymus, and those that were identified were generally found in the immature clusters. The increased frequency in B6.SAP^−/− E17 thymus of TRGV4/TRDV7 γδ T cells in the Il2ra⁺ c5 cluster coupled with their decreased frequency in the emigrating S1pr1⁺ c4 cluster suggests a model where SAP promotes progression through the early stages of development through thymic egress. SAP deficiency, therefore, could result in lower Vγ4/Vδ7 output. However, this model is at odds with our finding that the frequency of S1pr1⁺ clusters in neonatal thymus is increased in B6.SAP^−/−, and that these clusters contained increased numbers of TRGV4/TRDV7 γδ T cells. A possible explanation could be found in a previous report demonstrating that a Vγ4/Vδ7 subset increases in frequency in the periphery of B6 mice over time and undergoes extrathymic expansion⁷⁶. Our finding here that that the Vγ4/Vδ7 T cells represent a SAP-dependent innate-like γδT1 subset raise the possibility that SAP regulates the homeostatic expansion of this subset after it leaves the thymus. The mechanism, however, through which SAP is regulating the number of Vγ4/Vδ7 T cells remains unclear, as we have so far observed no evidence of changes in homeostatic proliferation or apoptosis of these cells in the lungs.

We previously demonstrated that the SAP-dependent impairment in γδT17 function in the neonatal periphery largely corrects itself as mice mature²⁴. This suggested the possibility that SAP-dependent γδT17 clonotypes originating in the embryonic thymus were being replaced over time by SAP-independent γδT17 clonotypes. Although our data certainly point to a SAP-dependent decrease in a germline-encoded Vγ4/Vδ5 clonotype in neonatal thymus, the low numbers of this clonotype in the lung makes it difficult to draw definitive conclusions in the periphery. Still, we did observe a role for SAP in shaping the peripheral γδT17 repertoire, as we noted a decreased frequency of germline-encoded TRGV4 sequences (so-called “GYSS”⁷⁷) coupled with an increased frequency of TRGV4 sequences with more diverse CDR3 (“GXYSS”) in the periphery. Although the potential significance of this change on γδ T cell function is unclear, it does support a model in which SAP regulates γδT17 developmental programming during embryonic thymic development, where the low expression of terminal deoxynucleotide transferase results in a high frequency of germline-encoded V-J junctions.

In summary, our findings reveal the presence of multiple SAP-dependent developmental checkpoint during the very early stages of γδ T cell development in the thymus, indicating a more complex role for this signaling pathway in γδ T cell developmental programming than previously appreciated. Indeed, our data indicate that SAP functions to regulate entry into the γδ and αβ T cell developmental pathways, and we make the striking observation that nearly all Vγ4 γδT1 cells in the periphery preferentially use Vδ7 chains and that it is a decrease in the Vγ4/Vδ7 subset that, in part, underlies the SAP-dependent decrease in γδ T cell IFN-γ production. Altogether, these findings suggest that SAP regulates the balance of innate-like γδT1/γδT2/γδT17 subsets during thymic development and in the periphery.

Methods

Mice

C57BL/6J (B6) and B6.Sh2d1a^−/− (SAP^−/−) mice were purchased from the Jackson Laboratory (Bar Harbor, ME) and were housed and bred in a specific pathogen-free facility at the animal facility of the University of Vermont. For experiments using E17 mice, timed-pregnant matings were used. The presence of a vaginal plug indicated day 0 of pregnancy. The flow cytometry experiments included 8-13 weeks old B6 and B6.SAP^−/− mice (age and sex-matched). All experimental procedures on animals were carried out with the approval of the University of Vermont Institutional Animal Care and Use Committee.

Tissue Processing

Single-cell suspensions from the spleen and thymus were obtained by gently passing tissue through nylon mesh. Red blood cells (RBCs) were lysed with Gey’s solution. Single-cell suspensions from the lungs were prepared as described (Benoit et al., 2015). Briefly, dissected lungs were finely minced in a digestion buffer containing DMEM, 1 mg/ml collagenase type IV (Gibco), and 0.2 mg/ml DNase (Sigma-Aldrich). After shaking at 200 rpm at 37°C for 20 min, the lung digests were triturated with a 16-gauge needle, followed by an additional 20 min shaking at 200 rpm at 37°C. Following the final trituration, red blood cells (RBC) were lysed with Gey’s solution. Finally, the digested lung tissues were filtered through a 70-μm filter and washed in PBS/2%FBS. Cells were counted either on a MACSQuant VYB or using a hemacytometer following trypan blue staining.

Flow Cytometry

Cells were incubated with LIVE/DEAD™ Fixable Blue Dead Cell stain (Thermo Fisher Scientific, Grand Island, NY) for 30 min at 4°C in PBS, after which they were resuspended in Fc-Block (BioLegend) containing PBS/2% FCS/0.1% sodium azide buffer. After a 5-10 min incubation on ice, antibodies specific for cell surface markers were added, incubated for 30 min on ice, and washed in PBS/2% FCS/0.1% sodium azide buffer. Brilliant Violet Staining buffer (BD Biosciences) was used in experiments contain more than one Brilliant Violet, Brilliant UV, or SuperBright-conjugated Ab. Cells were fixed in PBS/1% PFA buffer and stored at 4° C prior to data collection.

For nuclear staining, cells were permeabilized with Foxp3 transcription factor fixation/permeabilization buffer (ThermoFisher) following surface staining. Cells were incubated overnight at 4°C, washed, and incubated with 25 micrograms/ml total rat IgG before adding Abs against transcription factors. After staining, cells were washed and resuspended in PBS/2% FCS/0.1% sodium azide buffer. For intracellular cytokine staining, cells were permeabilized with PBS/1% FCS/0.1% sodium azide/0.1% saponin following surface staining and fixation. To reduce unspecific binding, 25 μg/ml total rat IgG was added prior to cytokine staining with fluorescent-labeled antibodies.

Flow cytometry data were collected on a Cytek® Aurora (Cytek Biosciences) spectral cytometer, and data analysis was performed using FlowJo Software (Tree Star). In some cases, visualization of flow cytometry data was conducted by concatenation of individual .fcs files, followed by dimensionality reduction and clustering using UMAP and FlowSOM⁷⁸ in FlowJo. Antibodies used in these experiments are anti-CD4 (RM4-5), anti-CD11b (M1/70), anti-CD11c (N418), anti-CD19 (1D3), anti-CD24 (M1/69), anti-CD25 (PC61), anti-CD27 (LG.3A10), anti-CD44 (IM7), anti-CD45RB (C363-16A), anti-CD73 (TY/11.8), anti-CD150 (TC15-12F12.2), anti-CD196 (29-2L17), anti-CD319 (4G2), anti-cMAF (sym0F1), anti-IFNγ (XMG1.2), anti-IL17a (TC11-18H10.1), anti-IL4 (11B11), anti-PLZF (9E12), anti-T-bet (4B10), anti-TCRβ (H57-597), anti-Vγ5 (536), and anti-Vd6.3 (C504.17c) from BioLegend; and anti-RORγt (Q31-378), anti-TCRδ (GL3), anti-Vγ1 (2.11) from BD biosciences; and anti-CD8a (53-6.7), anti-CD45 (30-F11), anti-CD200 (OX90), anti-CD278 (7E17G9), anti-Ki-67 (SolA15), anti-Ly108(13G3-19D), anti-Vγ4 (UC3-10A6) from eBiosciences; anti-S1P1 (713412) from R&D systems; anti-Blk (polyclonal) from Cell Signaling Technology; and anti-17D1

In vivo proliferation and apoptosis

For in vivo proliferation analysis, each mouse was intraperitoneally injected with 1 mg of 5-bromo-deoxyuridine (BrdU) (Sigma-Aldrich) each day for 3 days. Lung cell suspensions were prepared as above, stained with surface markers, after which they staining with anti-BrdU antibody according to the manufacturer’s instructions (BrdU Flow Kit; BD Biosciences). For apoptotic cell detection, annexin V staining was performed. Briefly, surface staining was followed by dead cell staining with Fixable Viability Dye eFluor™ 780 dye. The cells were washed and resuspended in an annexin-binding buffer (10 mM HEPES, 140 mM NaCl, and 2.5 mM CaCl2, pH 7.4). Finally, 1 μl of Annexin V Conjugate (Thermo Fisher Scientific, Grand Island, NY) was added to each tube, and cells were incubated for 15 minutes at room temperature (RT) prior to data collection.

Bulk RNA Sequencing and data analysis

For the bulk RNA sequencing, thymocytes from 10-day old B6 (n=4 male, n=2 female) and B6.SAP^−/− (n=5 male, n=1 female) mice were prepared by depleting CD4+ cells from thymocytes using the EasySepTM Mouse CD4 Positive Selection Kit II (STEMCELL Technologies). Cells were then incubated with Fixable Viability Dye eFluor™ 780 (Thermo Fisher Scientific, Grand Island, NY) for 30 min at 4°C in PBS, after which they were resuspended in Fc-Block (BioLegend) containing PBS/2% FCS. After washing, cells were stained with antibodies against CD11b (M1/70), CD19 (1D3), CD24 (M1/67), CD45 (30-F11), TCRβ (H57-597), and TCRδ (GL3) for 30 min at 4°C. Following surface staining, CD24^hiTCRδ+ cells were sorted from a CD11b⁻CD19⁻TCRb-population of live CD45⁺ lymphocytes. Immature (CD24^high) thymocytes from both B6 and B6.SAP^−/− mice were sorted directly into 350 microliters of RLT lysis buffer (Qiagen). RNA was extracted using the RNeasy Micro kit (Qiagen) according to manufacturer instructions. After RNA quantification and quality assessment, cDNA libraries were prepared using SMARTer Stranded Total RNA-Seq Kit v3 - Pico Input Mammalian (Takara, Japan) according to the manufacturer’s protocol. Following library clean-up and quantification, the libraries were sequenced (single-end 75 bp) on an Illumina HiSeq1500.

The Fastq files were checked for data quality using FastQC and MultiQC⁷⁹ bioinformatics tools. Cutadapt⁸⁰ was used to filter out low-quality sequences and to remove adapter sequences from all reads. The expression of transcripts was quantified using the mapping-based mode of Salmon⁸¹, which quasi-mapped RNA-seq reads to an index created from a mouse reference transcriptome (GRCm39). The transcript abundance estimates were imported using tximport⁸² and differential gene expression (DE) analysis was performed using DEseq2⁸³ employing the Wald test of the negative binomial model coefficients to determine statistical significance. The results were analyzed based on Log₂ fold change (Log₂FC ≥0.5 or ≤-0.5) and adjusted P-value (padj ≤0.0001) to determine whether a gene was differentially expressed.

Cell sorting and library preparation for CITE-seq and V(D)J enrichment

For thymus single-cell transcriptome analysis, we used embryonic day 17 (E17) B6 (n=4) and B6.Sh2d1a^−/− (SAP^−/−) (n=4) mice, nine day-old (D9) B6 (n=3) and B6.SAP^−/− (n=3) mice, and six-week-old B6 (n=3 male) and B6.SAP^−/− (n=2 male, n=1 female) mice. For single-cell analysis of adult lung γδ T cells, we used 13 week old B6 and B6.SAP^−/− littermates (pooled cells from 2 males and one female mouse in each group).

Analysis of gene expression, surface protein expression, and TCR V(D)J repertoire was performed using the 10X Genomics Chromium Next GEM Single Cell 5’ V(D)J with feature barcoding kit (version 1.1). To sort cells for scRNAseq, thymocytes were stained with anti-TCRδ-PE (GL3) antibody and enriched using anti-PE MicroBeads kit (Miltenyi Biotec). The exception was E17 thymus which did not undergo TCRδ enrichment. Following live/dead staining with Fixable Viability Dye eFluor™ 780 (Thermo Fisher Scientific, Grand Island, NY) for 30 min at 4° C in PBS, cells were resuspended in Fc-Block (BioLegend) containing PBS/2% FCS. Cells were then simultaneously stained with oligonucleotide-conjugated antibodies targeting CD24, CD44, CD45RB, CD73, SLAMF1, SLAMF6, Vγ1, and Vγ4 (TotalSeqC-BioLegend), oligonucleotide-conjugated hashtag antibodies (TotalSeqC-BioLegend), and fluorophore-conjugated antibodies for cell sorting which included the anti-CD11b (M1/70), anti-CD11c (N418), anti-CD19 (1D3), anti-CD45.2 (104) and anti-TCRβ (H57-597) antibodies. Following surface staining of the enriched population, TCRδ⁺ cells were sorted from a CD11b⁻CD11c⁻CD19⁻TCRβ⁻ population of CD45.2⁺ live lymphocytes.

Libraries were prepared according to the manufacturer’s instructions with modifications to the V(D)J enrichment procedure. Equal numbers of live GL3⁺ γδ T cells from B6 and B6.SAP^−/− mice were pooled and captured (between 5000 – 7000 cells) using a 10x Chromium controller. Following cDNA preparation, the gene expression and CSP libraries were created according to the manufacturer’s instructions. To obtain TCR gamma and delta chain sequences, we replaced kit-provided enrichment primers with custom mouse γ and δ chain-specific primers and modified the number of PCR cycles used from manufacturer recommended cycling conditions. In the first enrichment step, an equal amount of cDNA was used as a template for PCR amplification of γ and δ chains separately using a total of 12 and 10 PCR cycles, respectively. The primers used for the first γ-chain enrichment are: Forward Primer 5’-AATGAT ACGGCGA CCACCG AG ATCTACACTCTT TCCCTACAC GACGCTC-3’ (2 μM) and Reverse Outer Primers 5’-GGAA AGAACTT TTCAAGGAS ACAAAG-3’ (1 μM), 5’-CCCTTATG ACTTCAGG AAAGAA CTTT-3’ (0.5 μM). The primers used for the first δ-chain enrichment are: Forward Primer 5’-AATGAT ACGGCGACCACCG AGATCTACACTC TTTCCCTAC ACGACGCTC-3’ (2 μM) and Reverse Outer Primer 5’-CCACAATC TTCTTGGATG ATCTGAG-3’ (0.5 μM). Following PCR product purification, the γ-chain and δ-chain amplicons went through the second round of enrichment separately using a total of 12 and 10 PCR cycles, respectively. The primers used for the second γ-chain enrichment are: Forward Primer 5’-AATGA TACGGCGACCA CCGAGATCT-3’ (0.5 μM) and Reverse Inner Primers 5’-ACAAA GGTATGTCCCA GTCTTATGGA-3’ (0.5 μM), 5’-GGAGACA AAGGTAGGT CCCAGC-3’ (0.5 μM). The primers used for the second δ-chain enrichment are: Forward Primer 5’-AATGATA CGGCGACCACC GAGATCT-3’ (0.5 μM) and Reverse Inner Primer 5’-GTCACCTC TTTAGGGTAG AAATCTT-3’. Following PCR product cleanup, the second round PCR products (γ-chain and δ-chain amplicons) were quantified and combined in equal quantity making the final TCR V(D)J library. Finally, all three libraries (gene expression, cell surface protein expression, and TCR V(D)J) were sequenced (paired-end 26 x 91 bp) on an Illumina HiSeq 1500 sequencer or a NextSeq 2000.

CITE-seq quality control, filtering and demultiplexing

Raw sequence reads from the gene expression and CSP libraries were processed using the Cell Ranger software (v4.0.0)⁸⁴. Reads were mapped to the GRCm38 (mm10) mouse reference genome. To exclude the possibility of TCR transcripts playing a role in UMAP clustering, all TCR genes were omitted from the analysis using the “cellranger reanalyze” argument of Cell Ranger software.

Quality control, filtering, UMAP projection, and clustering analysis were conducted using Seurat (v4.3.0.1)⁸⁵ R package. Briefly, cells with a low number (≤2000 for E17 thymus dataset, ≤1200 for D9 thymus dataset and ≤750 for lung dataset) of unique genes (nFeature_RNA), a frequency of mitochondrial reads > 4% and a frequency of ribosomal protein < than ∼12% were eliminated from further analysis. To reduce cell cycle-related and stress-related heterogeneity in the normalized data, selected cell cycle genes and stress-induced genes ⁸⁶ were regressed out during the data normalization step. Note that in some instances, contaminating clusters expressing macrophage specific marker genes were manually removed using Seurat “CellSelector ()” argument. We used the Seurat function “HTODemux()” to demultiplex single cells back to their sample of origin ⁸⁷.

CITE-Seq and V(D)J enrichment data analysis

The demultiplexed datasets were further analyzed for gene expression, cell surface protein expression following the multimodal analysis pipeline of Seurat⁸⁵. For each dataset, we performed normalization and variance stabilization of the gene expression data using the scTransform⁸⁸. Following identification of top significant principal components (PCs) using the principal component analysis (PCA), we performed a graph-based (UMAP) semi-supervised clustering of cells with similar gene expression patterns using specific number of top PCs (E17 = 25, D9 = 30, W6 = 27, lung = 28) and resolution values (E17 = 1.2, D9 = 1.6, W6 = 0.65, lung =1.2) for each dataset. The cell surface protein expression data (Centered Log Ratio (CLR) normalized) from CSP libraries as well as the gene expression data were then visualized on the clustered cells using the R package scCustomize⁸⁹. The differentially expressed features of each cluster (cluster biomarkers) were identified using FindAllMarkers() function using the Wilcoxon Rank Sum test. For identifying differentially expressed genes (DEGs) between B6 and B6.SAP^−/− samples, we performed a pseudobulk DE (differential expression) analysis of aggregated read counts using the DEseq2 package, which accounts for biological replicates making it a better approach than other specialized single-cell DE analysis methods⁹⁰. Briefly, for cells of a given sample type, we aggregated RNAseq reads across biological replicates identified by individual hashtags. Then transformed the genes-by-cell matrix to a genes-by-replicates matrix and performed DE analysis using the DESeq2⁸³ package employing the Wald test of the negative binomial model coefficients to determine statistical significance. Finally, the trajectory and pseudotime inference analysis was performed using the partition-based graph abstraction (PAGA)⁹¹ analysis, which is a part of the Scanpy⁹² python package.

For the analysis of V(D)J library sequences, Cell Ranger v6.1.2 (“cellranger vdj” argument) was used to map the TCR V(D)J) reads to a custom-made reference containing sequences from the IMGT database⁹³. The Cell Ranger V(D)J output was then processed with custom bashscripts. Briefly, the TRGV and TRAV/TRDV contigs identified by Cell Ranger V(D)J were demultiplexed to individual barcodes after which the dominant TRGV and TRAV/DV transcripts were selected per cell. Finally, the V(D)J information of each cell was integrated with scCITE-seq data and visualized using a combination of Seurat and custom-made R scripts. We note that since about half of Vγ1, Vγ4, and Vγ7 T cells contain functionally rearranged TRGV2 chains, but only a small percentage express Vγ2 on the cell surface⁹⁴, the methods used here cannot adequately assess the true TRGV2 repertoire. Therefore, we report it here for accuracy but have omitted TRGV2 from our analysis.

Quantitative PCR

Lung qPCR samples were prepared by positively selecting leukocytes using CD45 MicroBeads (Miltenyi Biotec). Spleen qPCR samples were prepared by depletion of CD19+ TCRβ+ cells using biotinylated antibodies and Streptavidin MicroBeads (Miltenyi Biotec) according to the manufacturer’s instructions. Total RNA was extracted from approximately 5-10 million cells from each organ using the RNeasy Mini Kit (QIAGEN). SuperScript® IV Reverse Transcriptase (Thermo Fisher Scientific, Grand Island, NY) was used to generate cDNA using 400 ng of RNA and oligodT primer according to the manufacturer’s instructions. qPCR was performed using 1μl of cDNA template in iTaq Universal SYBR Green Supermix (Bio-Rad) containing TCR delta chain specific primer pairs. The qPCR reaction was run on a Bio-Rad CFX96 qPCR machine using the following cycling conditions: 95°C for 30 sec; 40 cycles of [95°C for 5 sec; 60°C for 30 sec; 72°C for 30 sec]. The β2Μ and TRDC genes were used as endogenous housekeeping controls. Primers used were as follows: mb2mforward, 5’-CATGGCTCG CTCGGTGACC-3’, mb2mreverse, 5’-AATGTGAGGC GGGTGGAACTG-3’, TRDCforward, 5’-CTCCGGCCA AACCATCTGTT-3’, mouseDV2, 5’-CCAAGAAG CATACAAGCAGT ATAATG-3’, mouseDV5, 5’-CCCATGA TGCAGATTT TGTTCAAGG-3’, mouseDV7, 5’-GGAAG MCTCGTCAGC CTGTTGT-3’, mouseDCrev, 5’-CCACAATCTT CTTGGAT GATCTGAG-3’. Some primer sequences were adopted from Wei. et. al. (2015)⁵¹.

Vγ4 single cell sorting and TCR Sequencing

The plate-based based γδ TCR repertoire analysis of single-cell sorted Vγ4 cells involved adult lungs from 8-10 weeks old B6 male (n=7) and female mice (n=6). Lung Vγ4 TCR sequences were obtained using a protocol modified from Wei et al., (2015)⁵¹. Briefly, lung cell suspensions were prepared as described above. After washing, cells were resuspended in Fc-Block (BioLegend) containing PBS/2% FCS buffer, and stained with anti-CD11b (M1/70), anti-CD11c (N418), anti-CD19 (1D3), anti-CD45 (30-F11), anti-TCRδ (GL3) and anti-Vγ4 (UC3-10A6) from for 30 min at 4°C in PBS. Fixable Viability Dye eFluor™ 780 (Thermo Fisher Scientific, Grand Island, NY) was used to identify dead cells. Single-cell sorting of lung Vγ4⁺TCRδ⁺ from a CD11b⁻ CD11c⁻CD19⁻TCRβ⁻ population of CD45⁺ live lymphocytes were performed using a FACS Aria III cell sorter (BD Bioscience). Cells were sorted directly into 12.5 μl (per well) of OneStep RT-PCR Buffer mix (QIAGEN) in a 96-well PCR plate (one cell/well).

TCR sequences were then amplified and barcoded according to (Wei et al., 2015). Briefly, the first two nested PCR reactions involved primer-specific amplification of gamma and delta TCR sequences. In both the first and second PCR reactions, the final concentration of each V-region primer is 0.36μM, and each C-region primer is 0.6μM. The first RT-PCR step was performed using the following cycling conditions: 50°C for 30 min; 95°C for 5 min; 25 cycles of [94°C for 30 sec; 62°C for 1min; 72°C for1 min]; 72°C for 7 min; 4°C. Next, 1 μl of the first round PCR product was used as a template for the second nested PCR reaction using HotStartTaq Master Mix Kit (QIAGEN). The cycling conditions for the second PCR reaction were 95°C for 5 min; 25 cycles of [94°C for 30 sec; 64°C for 1min; 72°C for1 min]; 72°C for 7 min; 4°C. In the third PCR reaction, 1 μl of second round PCR product as the template in a new 14 μl reaction using HotStartTaq Master Mix Kit (QIAGEN). Each well was barcoded. The final concentration of 5’ and 3’ barcoding primers was 0.05μM each and paired-end primers were 0.5μM each. The cycling conditions for the third PCR reaction were 95°C for 5 min; 36 cycles of [94°C for 30 sec; 62°C for 1min; 72°C for1 min]; 72°C for 7 min; 4°C. The γ and δ chain products were run on 1.2% agarose gel separately and purified using the Qiagen QIAquick Gel Extraction Kit according to kit instructions. The final products were then analyzed using the 2100 Bioanalyzer (Agilent Technologies). Finally, the γ and δ TCR products were mixed in equal concentration and sequenced on the Illumina MiSeq platform using the Illumina MiSeq Reagent Kit v2 (500-cycles).

The paired-end TCR sequencing reads were analyzed using a custom software pipeline capable of de-multiplexing barcoded sequences to individual cells. The resulting de-multiplexed reads were then analyzed by MiGMAP software, which identifies CDR3 sequences using the IgBlast 1.4.0 (NCBI) V(D)J mapping tool^{95, 96}. The single-cell TCR V(D)J profiling data were then visualized using custom R scripts.

Statistical analysis

Prism (GraphPad Software, San Diego, CA) was used for the statistical analysis of data. Unpaired Student t-tests, one-way ANOVA, or two-way ANOVA were used where appropriate. Tests were considered significant when p-values were less than 0.05. Error bars represent standard deviation (SD) unless otherwise indicated.

Acknowledgements

We thank Marlana Winschel, Cameron Moquin, and Erin Mathieu for help in tissue preparation and flow cytometry. We acknowledge Roxana del Rio Guerra for help in cell sorting. All flow cytometry data were collected at the Harry Hood Bassett Flow Cytometry and Cell Sorting Facility (Larner College of Medicine, University of Vermont) and was supported by NIH grant (S10OD026843). We thank Jessica Hoffman, Kris Finstaad, and the Vermont Integrated Genome Resource (Larner College of Medicine, University of Vermont) for help in scRNAseq, Julie Dragon and Korin Eckstrom for help in bioinformatics analysis. We thank Fred Kolling at the Genomics Shared Resource (Dartmouth Geisel School of Medicine) for help in sequencing. We also thank Willi Born and Rebecca O’Brien for the 17D1 hybridoma and Pablo Pereira for helpful discussion. This work benefitted from data assembled by the ImmGen consortium, and was supported by NIH grants R21AIO66465 and R03AI153902, S10OD026843 (JEB) and an AAI Careers in Immunology Fellowship (SKM, JEB).

Author contributions

S.K.M. and J.E.B. conceived and designed the project; S.K.M., J.E.B., E.A., K.J.H., B.M.H., R.S., and O.D. performed experiments and analyzed data; K.J.H., D.G., J.T.R., and N.S. aided with the single-cell TCR sequencing; D.M. assisted with bulk RNA sequencing; S.K.M. and J.E.B. performed the bioinformatics analysis of bulk RNA-seq and single-cell CITEseq (with immune profiling) data; J.E.B. supervised the project; S.K.M. and J.E.B. prepared the manuscript.

Data Availability

All sequencing data in this manuscript have been deposited in NCBI GEO and have been assigned the accession number GSE262064. All scripts used in these analyses have been made publicly available on Github and can be accessed here: https://github.com/Boyson-Lab/scRNAseq_with_immune_profiling_tutorial

Competing interests

All authors declare no competing interests.

Supplementary Figure Legends

Supplemental Figure 1:

Supplemental Figure 2:

Supplemental Figure 3:

Supplemental Figure 4:

Supplemental Figure 5:

Supplemental Figure 6:

Supplemental Figure 7:

Supplemental Figure 8:

Supplemental Figure 9:

Supplemental Figure 10: