Phased differentiation of γδ T and T CD8 tumor-infiltrating lymphocytes revealed by single-cell transcriptomics of human cancers.

γδ T lymphocytes diverge from conventional T CD8 lymphocytes for ontogeny, homing, and antigen specificity, but whether their differentiation in tumors also deviates was unknown. Using innovative analyses of our original and ~150 published single-cell RNA sequencing datasets validated by phenotyping of human tumors and murine models, here we present the first high-resolution view of human γδ T cell differentiation in cancer. While γδ T lymphocytes prominently encompass TCRVγ9 cells more differentiated than T CD8 in healthy donor's blood, a different scenario is unveiled in tumors. Solid tumors and lymphomas are infiltrated by a majority of TCRVγnon9 γδ T cells which are quantitatively correlated and remarkably aligned with T CD8 for differentiation, exhaustion, gene expression profile, and response to immune checkpoint therapy. This cancer-wide association is critical for developing cancer immunotherapies.

Introduction

The presence and functional differentiation of cytotoxic tumor-infiltrating lymphocytes (TILs) is critical for the clinical outcome of cancer patients. Although T CD8 TILs are associated with favorable prognosis in all types of cancers, little is known about other cytotoxic TILs of γδ T lineage. Human γδ T cells are CD4−CD8− T lymphocytes expressing either TCRVδ2 or TCRVδ1/3 chains prominently paired with TCRVγ9 or TCRVγ2/3/4/5/8 (TCRVγnon9), respectively. On the one hand like T CD8 cells, circulating γδ T lymphocytes differentiate into naïve (Tn), central memory (Tcm), effector memory (Tem), and terminally differentiated effector memory CD45Ra+ (Temra) cells.[1] Similarly, both can infiltrate tumors and mediate cytotoxic functions important for antitumor immunity. On the other hand, γδ T lymphocytes differ from conventional T CD8 by tissue distribution, relative abundance, and a more pleiotropic immunobiology.[2] Unlike T CD8, γδ T cells react to TCR-independent and TCR-dependent stimuli that differ notably from major histocompatibility complex (MHC)-restricted peptides. Peripheral TCRVγ9 γδ T lymphocytes are specific for phosphoantigen-sensing butyrophilins, while tissular TCRVγnon9 γδ T cells recognize more diversified antigens.

Whether T CD8 and γδ T lymphocytes infiltrating tumors are somehow associated remains unclear, but it is important to develop and monitor γδ-targeted cancer immunotherapy.[3] By deconvoluting 10,000 tumors from 55 cancer types, we and others have previously found that the abundance of TCRVγ9 TILs is unrelated to that of TCRαβ TILs, varies considerably across individuals and cancers, and is associated with better outcomes in cancers.[4-6] Likewise, the presence, function, and potential of TCRVγnon9 cells which represent another prominent subset of γδ TILs are demonstrated in different human tumors.[7-9] So far, however, neither the general landscape of both subsets of cytolytic γδ TILs nor their putative link with the T CD8 had been characterized both at the single cell level and at a cancer-wide scale.

To address this issue, the characterization of single-cell transcriptomes from γδ TILs requires to identify exhaustively and selectively these cells in multiple scRNAseq datasets from different human tumors. Although theoretically only the γδ T lymphocytes express a cell surface γδTCR, other cell subsets do actually express some level of TCRγ or TCRδ mRNAs. Both TCRαβ T CD8 lymphocytes and NK cells also express either of the TRDC, TRGC1, or TRGC2 genes, imposing the use of a multigene signature rather than single genes to reliably identify γδ T lymphocytes only.[10] Here, we incepted a ‘score and gate’ strategy that significantly improves the digital purification of γδ T lymphocytes from scRNAseq datasets, and developed a cell pseudotime-based visualization of differentiation trajectories. With the γδ T and the T CD8 lymphocytes extracted from ~150 human samples of PBMC and tumor biopsies, the pseudotimed differentiation trajectories of these lineages revealed their strikingly contrasting hallmarks in blood and cancer.

Results

Circulating γδ T are more differentiated than T CD8 lymphocytes in healthy individual’s blood

From seven scRNAseq datasets of healthy donor’s PBMC and γδ T lymphocytes cell-sorted from PBMC,[10] (n = 6,122) γδ T lymphocytes were digitally extracted by ‘score and gate’ (Supplemental Data Figure S1) (Methods). Then, a differentiation trajectory of these extracted cells was computed, and the trajectory’s principal components were projected onto the linearized pseudotime axis to visualize a pseudotimed trajectory (Methods). The pseudotimed differentiation trajectory of peripheral blood γδ T lymphocytes comprised 14 milestones (Figure 1a) in which TCRVγ9 were prominent over TCRVγnon9 cells, as expected (not shown). This trajectory was anchored by expression of Tn genes including the lymph node homing receptor CCR7[11] and the TCF1 differentiation transcription factor-encoding gene TCF7,[12] which peaked in the M1, M2 milestones. Further M3-M7 displayed intermediate levels of Tn and Tem gene expression, a defining hallmark of Tcm lymphocytes.[13,14] The Tcm M6 and M7 upregulated both Tn (CD27, CD28, SELL, IL7R) and tissue migration genes (CD44, CXCR3, CXCR4). The Tem milestones M10-M11 included activated cells (IFNG, CCL5), and cytotoxic Tem cells (NKG7, PRF1, GZMA, GZMB, GZMM, FASLG) in M10-M13. The maximal expression of FCGR3A (CD16 marker), PRF1, and GZM genes defined the Temra cells,[15] in M14. The exhaustion transcription factor TOX[16] and immune checkpoint genes LAG3, TIGIT, and PDCD1 were weakly expressed in control γδ T cells (Figure 1b). When projected onto the trajectory, the scores for naïve (Tn) and effector (Teff) gene signatures evolved smoothly and inversely, reflecting the transcriptional progression along aging from naïve to Temra cells (Figure 1c). Projecting on this trajectory, several published T CD8 differentiation signatures (Table S1) confirmed the differentiation stage of each milestone (Supplemental Data Figure S2), whose ontogeny was validated by the dendrogram of milestones (Figure 1d).

Figure 1.

Circulating γδ T are more differentiated than T CD8 lymphocytes in most healthy individual’s blood. a, pseudotimed trajectory of PBMC-derived γδ T cells. Cell color key: pseudotime-ordered milestones. b, mean of gene expression per milestone. Genes shown blue and red: naïve (Tn) and effector (Teff) signatures, respectively. c, Single cell scores for Tn and Teff signatures across the trajectory. d, dendrogram of milestones (centroids, same color key as a) overlaid to the pseudotimed trajectory showing all γδ T cells (gray dots), marker genes, trajectory (dark line), and differentiation direction (lines). e-h, same as a-d for T CD8 cells extracted from the same PBMC. i, rates of γδ T and T CD8 cells at each differentiation stage. Abbreviations: Tn: naïve T, Tcm: central memory T, Teff: T effector, Tcm-mig: migration-prone central memory T, Tem: effector memory T, Tem-act: activated effector memory T, Tem-cyto: cytotoxic effector memory T, Temra: effector memory CD45RA T, Tscm: stem-like central memory T, Tpex: progenitor exhausted-like T

With 15 milestones, the pseudotimed trajectory of the (n = 3,680) PBMC-derived T CD8 lymphocytes approximated that of γδ T cells (Figure 1e), showing cells of stem central memory profile (IL7R, CD28)[17] in Tcm M4, and cells with progenitor exhausted-like (Tpex) profile (CD95, GZMK, TOX, TIGIT)[18] in M6. The activated Tem cells (IFNG) up-regulated PDCD1 and LAG3 in M11, while the T CD8 Temra (FCGR3A) cells peaked at final M15 instead of final M14 in PBMC (Figure 1f). These identifications were supported by projecting the above Tn, Teff, and other published signatures (Figure 1g, Supplemental Data Figure S3) and by the dendrogram of T CD8 milestones (Figure 1h).

Furthermore, (n = 137) γδ T and (n = 725) T CD8 lymphocytes were extracted from our CITEseq dataset of healthy donor’s PBMC immunostained for differentiation markers,[19] and were mapped onto the above trajectories. Their differentiation stages determined by the cell surface CD45RA and CD62L immunophenotypes were consistent with those assigned by gene signatures (Supplemental Data Figure S4). Reflecting distinct biologies and antigen specificities, however, the proportions of cells at each differentiation stage differed between the lineages and across individuals: γδ T were prominently Tem, while Tn predominated more often among T CD8 (Figures 1i, Supplemental Data Figure 4c). Therefore, in healthy adults PBMC, the γδ T are generally more advanced than T CD8 lymphocytes on their respective differentiation trajectories.

Convergent differentiations of γδ T and T CD8 lymphocytes infiltrating tumors

Then, to compare the differentiation hallmarks of γδ T and T CD8 lymphocytes infiltrating tumors (TILs), we performed 3ʹ chemistry scRNAseq of malignant lesions from (n = 3) splenic marginal zone lymphoma (SMZL) (GSE175785), and (n = 4) ovarian carcinoma (OVCAR) (GSE148162) cancer patients. Additional scRNAseq datasets from (n = 129) additional human tumors of various cancer types were also downloaded (summarized in Table S2), normalized, and integrated for comparison (Methods). Altogether, a total of (n = 4,648) γδ T and (n = 45,174) T CD8 TILs were extracted from a total of 136 tumors (Tables S3–S6), and analyzed as above (Supplemental data).

The tumors richest in γδ T TILs corresponded to the richest in T CD8 TILs (Fischer Exact p = .002), and the richest in TCRVγ9 were also the richest in TCRVγnon9 cells (Fischer Exact p = .0001). The trajectories of each lineage of TILs superimposed to those of PBMC (Figure 2a), but some milestones enriched in cancer revealed that both lineages of TILs (Figure 2b) were more differentiated than the PBMC (Figure 1i). The cell counts of γδ T and T CD8 TILs were strongly correlated (Pearson r = 0.7–0.99) in most cancer types except in lung carcinoma (LUAD) and melanoma (MEL). In Hodgkin’s lymphoma (HL) and head and neck carcinoma (HNSCC), the viral status of patients unveiled a higher correlation of the γδ T and T CD8 TILs counts in the EBV-positive and HPV-positive, respectively, than in their negative counterparts (Supplemental Data Figure 5).

Figure 2.

Memory γδ T and CD8 T cells are enriched in 136 human tumors from 11 cancer types. a, pseudotimed trajectories of TILs (red) of the (n = 4,648) γδ T and (n = 45,174) T CD8 cells featuring their dendrogram of milestones, overlaid to their PBMC-derived counterpart (green) shown on their respective trajectory. b, Rates of γδ T (blue) and T CD8 (orange) TILs with the specified differentiation stage in representative tumor samples. c, TIL cell counts and rates in each tumor sample (colored by cancer type). d, bubble plot of Tex and Ttrm cell counts per tumor sample, showing cancer type (color key as c), and infiltration pattern. e, representative example of trajectory maps of TILs from the inflamed tumor HNSCC#26, overlaid on the PBMC-derived cells (gray). Color key: TCRVγ9 cells (light blue), Tex TCRVγ9 cells (dark blue), TCRVγnon9 cells (orange), Tex TCRVγnon9 cells (dark red), T CD8 (yellow), Tex T CD8 (black). f, rates of Tex cells per differentiation stage and T cell lineage. Cell counts are indicated. g, violin plots and single cell scores of glycolytic gene enrichment in exhausted and non-exhausted cells from the specified group (Mann–Whitney P values). h, association of exhaustion with tissue residency in the Tcm and Tem stages of T γδ and T CD8 TILs (χ2 P values), absolute counts are specified. Abbreviations: Tn: naïve T, Tcm: central memory T, Tem: effector memory T, Temra: effector memory CD45RA T, Tex: exhausted T, Ttrm: tissue-resident memory

Classifying the tumor infiltration extents by TIL cell counts (Methods) indicated that the OVCAR and HNSCC were frequently inflamed (I) while the SMZL and hepatosplenic carcinoma (HCC) had classic (C) infiltrates, contrasting with the immune deserts (D) more frequent in many solid tumors (e.g. LUAD, MEL). The γδ T/T CD8 TIL ratio decreased from desert to inflamed infiltrates (mean D: 24%, C:10%, I: 5%) (Figure 2c, Table S2b). Identifying the Tex and Ttrm cells in all tumors (Methods) showed their respective rates were correlated between the two lineages (Pearson r = 0.8 in both cases) (Figure 2d). Confirming recent reports,[7-9] the majority of γδ TILs cancer-wide was of TCRVγnon-9 subtype, reflecting the extra-lymphoid tropism of this subset (Figure 2e, Table S2b). Combining trajectory and signature analyses revealed that, like Tex T CD8 TILs, all the Tex γδ TILs displayed Tem or Tcm differentiation stages (Figure 2e, f). Comparatively to their non-exhausted counterparts, these γδ T and T CD8 Tex cells shared a signature of 14 up-regulated genes involved in T cell exhaustion (Supplemental Data Figure 6a, b). These lymphocytes were also more glycolytic and more abundant among, but not restricted to, the Ttrm cells (Figure 2g, h, Supplemental Data Fig.6 c). In the non-exhausted compartment, the few Tn and early Tcm γδ TILs presented gene signatures of resting unexpanded T CD8 TIL clones (tumor singletons),[20] and the most mature Tem γδ T TILs had signatures of peripherally expanded large T CD8 clones[21] associated with response to immune checkpoint blockade (ICB) therapy[20,22] (Supplemental Data Figure 6c, d). As a result, the single-cell composition of all cancer samples showed that the differentiated states of γδ T and T CD8 were far more coherent in TILs than in PBMC.

Ex vivo immunophenotypes of γδ T and T CD8 TILs are correlated

To validate these conclusions, we obtained tumor biopsies from five patients with untreated ovarian cancer (OVCAR, n = 2), head and neck squamous cell carcinoma (HNSCC, n = 1) and cervical cancer (CC, n = 2), isolated their TILs, and immunostained the cell surface markers of γδ T and T CD8 lymphocytes for flow cytometry analysis.

Consistent with the scRNAseq data, both lineages of TILs from these lesions were prominently CCR7CD45RA Tem cells (Figure 3a). In these epithelial cancers, the tumor antigen-specific Tex CD8 TILs associated with response to ICB are highly cytotoxic PD-1hiTIM-3+ cells co-expressing CD39 and Ttrm markers.[23] So, we analyzed these features at the cell surface of the γδ T and T CD8 TILs from the above biopsies, and found sizable fractions of PD-1hiTIM-3+ cells and of PD-1hiCD39+ cells among both lineages (Figure 3b). Despite the previously depicted variability of TILs across tumors (Figure 2b), the samples with high Tex rates in T CD8 TILs had also high Tex rates in γδ T TILs (Methods) and the same applied for samples with low Tex rates (Figure 3b). In all biopsies, both the CD8 T and γδ T lymphocytes encompassed PD-1hiTIM-3+ and PD-1hiCD39+ cells, and across these samples, their respective rates of T CD8 and γδ T TILs were correlated (Figure 3c).

Figure 3.

The ex vivo phenotype of tumor-infiltrating γδ T and CD8 T lymphocytes consists of similarly differentiated cells a, right: Cell surface CCR7 and CD45RA phenotype of the γδ T and T CD8 TILs from a freshly biopsied HNSCC tumor, and differentiation stages in this tumor. Left: Rates of Tn, Tcm, Tem, and Temra cells in γδ T and T CD8 lineages from the five freshly biopsied tumors (means from two independent experiments). b, cell surface phenotype of PD-1 versus TIM-3 (left) and PD-1 versus CD39 (right) in γδ T and T CD8 TILs from the specified patients. Right: composite of results for rates of PD-1hiTIM-3+ and PD-1hiCD39+ cells among the γδ T and in T CD8 TILs (shown are results from three independent experiments from the five biopsies). c, Pearson correlation between the specified cells (same patients as above). d, e, single cell suspensions from TC1-engrafted murine HNSCC tumors (n = 12) were stained ex vivo for murine CD3, CD4, CD8, pan-TCRγδ, CD45.2, PD-1, TIGIT, CTLA-4 and TIM-3 and analyzed by flow cytometry. d, PD-1 versus TIM-3 expression in γδ T (top) and T CD8 (bottom) TILs from one mouse (representative result). e, proportions of subpopulations expressing the 16 possible combinations of PD-1, TIGIT, CTLA-4 and TIM-3 immune checkpoints in γδ T (top) and T CD8 (bottom) TIL. The five higher median proportions of TIL subpopulations are represented for γδ T (blue) and T CD8 TIL (orange). Composite of results means and sd. Abbreviations: Tn: naïve T, Tcm: central memory T, Tem: effector memory T, Temra: effector memory CD45RA T, Tex: exhausted T, Ttrm: tissue-resident memory

We further analyzed (n = 12) mice from a syngeneic murine model of HNSCC obtained by orthotopically implanting the lung carcinoma TC1 cells retrovirally transduced to express HPV16 E6/E7 antigens.[23] The TILs were isolated from established tumors for immunostaining and flow cytometry analyses. These assays showed the emergence of Tex γδ T and T CD8 TILs at the same time points of growth development, and the rates of PD-1hiTIM-3+ Tex cells in these γδ T and in T CD8 TILs were correlated (Figure 3d). In these mice, of the 16 possible co-expression patterns of the PD-1, TIGIT, CTLA-4, and TIM-3 immune checkpoints,[23] the same 5 phenotypic patterns predominated among both γδ T and T CD8 TILs (Figure 3e). Collectively, these ex vivo data confirmed the correlated differentiation patterns of γδ T and T CD8 TILs.

Both T CD8 and γδ T TILs are impacted by ICB therapy

The convergent differentiation of intratumoral γδ T and T CD8 lymphocytes reported above prompted us to investigate whether γδ T TILs are as affected as T CD8 TILs to immune checkpoint blocking (ICB) therapy. We downloaded two published scRNAseq datasets of 21 tumors from (n = 11) BCC[22] and (n = 10) MEL[24] patients-derived tumor biopsies sampled before and after ICB treatment. Their γδ T and T CD8 TILs were extracted as above, mapped on their respective differentiation trajectory, and analyzed with regard to patient’s clinical response, tumor infiltration pattern, and infiltrate’s reactivity to treatment (post-treatment TIL counts increase) (Methods).

In all samples, the T CD8 outnumbered the γδ T cells (Figure 4a). The pre-treated BCC lesions were more infiltrated than those of MEL (mean: 295 versus 86 TILs, respectively), and were more reactive to treatment (BCC: 8/11 reactive, mean counts: 306 (pre)/785 (post), MEL: 4/10 reactive, mean counts: 73 (pre)/149 (post)) (Figure 4a). In most BCC, the treatment increased the cell counts for both lineages and differentiation stages, including Tex, without affecting their relative rates of Tn, Tcm, Tem, and Temra cells (Figure 4b, c). Although the prominently desert infiltrates of MEL hampered further conclusions, the BCC results showed that ICB therapy strengthens the association of tumor infiltrating T CD8 and γδ T lymphocytes (Pearson r = 0.77 pre-Tt, 0.95 post-Tt). Among these latter, however, although similar counts of TCRVγ9 and TCRVγnon9 γδ cells infiltrated the pretreated BCC and MEL (on average 2.5 cells/pre-Tt sample), the TCRVγnon9 subset increased more significantly after therapy (on average: 5.3 TCRVγnon9 and 3.6 TCRVγ9 cells/post-Tt sample, Student p = .04).

Figure 4.

Both γδ T and CD8 T TILs respond to anti-PD1 therapy in patients with cutaneous cancers. a, summary of characteristics and TIL counts extracted from the specified BCC and MEL samples. Inflamed (dark red), classic (pink), desert infiltrates (white). b, representative examples of inflamed (BCC#009), classic (BCC#001), and desert (MEL#008) pre-treatment infiltrate patterns. The T TILs extracted from each specified tumor are overlaid on their respective pseudotimed trajectory (gray: control blood cells of the corresponding lineage; TCRVγ9 cells (light blue), Tex TCRVγ9 cells (dark blue), TCRVγnon9 cells (orange), Tex TCRVγnon9 cells (dark red), T CD8 (yellow), Tex T CD8 (black). Comparison of BCC samples before and after ICB treatment for: c, correlation of T CD8 and γδ T cell counts, d, distribution of differentiation stages in T CD8 and γδ T TILs. Abbreviations: Tn: naïve T, Tcm: central memory T, Teff: T effector, Tcm-mig: migration-prone central memory T, Tem: effector memory T, Tem-act: activated effector memory T, Tem-cyto: cytotoxic effector memory T, Temra: effector memory CD45RA T, Tscm: stem-like central memory T, Tpex: progenitor exhausted-like T

Finally, the patient’s clinical responses did not correlate to either of i) γδ T and/or T CD8 TIL counts alone or together, ii) TIL counts per milestone, iii) pre-treatment infiltration pattern, or iv) infiltrate reactivity to treatment, taken alone. Thus, beyond confirming that ICB therapy increases the tumor-infiltrating T CD8 and Tex regardless of clinical outcome,[22] our results demonstrated that such treatments affected the γδ T as well as T CD8 TILs.

Discussion

Understanding whether the functional differentiation of tumor-infiltrated γδ T lymphocytes is similar or unrelated to that of the T CD8 compartment is critical to implement γδ T cell-targeted cancer therapies. Here, we provided the most comprehensive cancer-wide atlas of these lymphocytes visualized in their continuous dynamics of differentiation. Thanks to the high resolution of pseudotimed trajectories, our study refined more and subtler transitions than the former four Tn-to-Temra stages of differentiation[15,25] in γδ T lymphocyte and T CD8 lymphocytes. By quantifying each cell’s progression along a continuous series of differentiation states, pseudotime represents a transform of true chronological time albeit through non-constant scales.[26-28] How explicit time is quantitatively converted into T cell pseudotime remains to be determined, but the related sizes of the γδ T and T CD8 trajectories suggests both lineages display the same time transform. Nonetheless, the ‘public’ trajectory paved here by γδ T lymphocytes from ~150 healthy and diseased individuals (Tables S3, S5) constitutes a reference roadmap of γδ T cell differentiation. Upon direct integration[29] of newer γδ T cells from any other scRNAseq dataset, each single cell coordinates on this trajectory will define its differentiation hallmarks, allowing automated characterizations based on this resource.

Although well characterized in the mouse,[30] the differentiation of peripheral human γδ T cell has remained less investigated over the last decade. A linear sequence of four successive and non-overlapping Tn to Temra stages initially established for TCRVγ9 cells[1] was assumed to also apply for TCRVγnon9 cells. This model was validated recently at the single cell level,[10] but neither its pseudotime scale nor milestones were known. Here, these features were defined for both subsets, and enlighted several subtle molecular transitions among Tcm and Tem toward stem-like, migration-prone, and exhausted cells, which so far were only characterized in the T CD8 lineage. Hence, most if not all T CD8 differentiation stages also exist in the TCRVγ9 and TCRVγnon9 subsets of γδ T lymphocytes. These subsets indeed have different tropisms at the steady state, with the predominance of TCRVγ9 cells in the blood and TCRVγnon9 cells in extra-lymphoid tissues.[31] This repartition disappears in tumors, however, where the TCRVγnon9 γδ TILs predominate in all (non-)hematopoietic cancer types.[7-9] Viral infections are major drivers of peripheral differentiation of T CD8 and γδ T lymphocytes. Thus, the locally phased differentiation of both γδ T and T CD8 lymphocytes strongly suggests a functional convergence toward effective antitumor immunity.[4,7,32] Substantiating this view, the highest correlation of γδ T and T CD8 TILs cell counts and differentiation observed in virus-positive patients is consistent with the repertoire of T CD8 TILs specificities prominently directed to viral antigens.[33] It is tempting to speculate that this functional convergence scenario underlies the better outcome of HPV over HPV HNSCC patients,[34] the Sars-CoV2-induced remission in EBV HL,[35] and the generally lower risk of cancer in CMV than CMV kidney transplant recipients.[36]

Since previous transcriptomic deconvolutions of human cancers showed a dissociation of γδ T and αβ T TILs,[4-6] the tight association of γδ T and T CD8 TILS unveiled here was unexpected. This discrepancy reflects not only the better resolution of scRNAseq transcriptomics than artificial intelligence-driven inferences from bulk microarray deconvolutions, but also the difference between T CD8 alone and the total αβ T compartment. Nevertheless, the coordination of γδ T and T CD8 TILs reported here might also reflect a broader orchestration, possibly involving a BTN3A1-mediated mechanism.[37] Furthermore, since T cell differentiation ensues antigenic activation, the advanced differentiation of TCRVγ9 and TCRVγnon9 γδ TILs questions the nature of their driving stimuli in tumor. Whether these are metabolites like bacterial phosphoantigens,[38] viral ligands,[39] stress-induced antigens,[40] or even HLA-restricted mutational peptide neoantigens remains to be determined more systematically.

The correlated exhaustions of intratumoral γδ T and T CD8 lymphocytes reflect a sustained activation and microenvironmental pressure in both cases. As a consequence of this coherence, the tumor infiltration by both cell types was affected by checkpoint blockade therapies. Although these are critical determinants of outcome, neither T CD8 TILs pre- or post-treatment[41] nor their exhaustion status[42-44] suffice to robustly predict the clinical responses of cancer patients. Hence, it is not surprising that the treatment-induced changes of γδ T and T CD8 TIL cell counts, alone or together, did not correlate either. In addition to immunotherapy-induced changes in the functional state of TILs,[45] the anatomical location of metastases, and tumor mutation of the IFNγ signaling or antigen-presentation pathways[46] are among the other determinants of therapeutic responses. Nevertheless, our recurrent detection of post-treatment TCRVγ9 γδ TILs in tumors without such cells before treatment suggests that ICB therapy harnesses a circulating γδ T cell contingent to the cancer immunity cycle.[47] Since therapeutic antibodies refuel tumor lesions with fresh, functional, and expanded lymphocytes from peripheral blood,[20] it is conceivable that the recruitment of any circulating subset of expanded γδ T lymphocytes allows intratumoral collaborations of synchronous γδ T and T CD8 effectors.[32]

Overall, our study builds a fundamental framework to characterize γδ T lymphocyte differentiation in human cancers, and represents an invaluable resource for the emerging cancer immunotherapies targeting γδ T lymphocytes.

Materials and Methods

Clinical patient sample collection

The SMZL samples were collected at the Institut Universitaire du Cancer de Toulouse-Oncopole (IUCT-O) and processed at the CRB Cancer from Hôpitaux de Toulouse, following the standard ethical procedures (Declaration of Helsinki), after obtaining written informed consent from all patients and upon approval by the Institutional Review Board (nDC-2009-989). The carcinoma samples were collected from patients with head and neck squamous cell carcinoma, ovarian cancer, and cervical cancer at the time of surgery at the Institut Universitaire du Cancer de Toulouse-Oncopole (IUCT-O) in accordance with the Declaration of Helsinki, upon approval by the Institutional Review Board (nDC-2016-2656 and DECIdE protocol, NCT03958240) and signed informed consent.

Mouse models

Female C57BL/6 mice, 7 weeks old, were purchased from Janvier Laboratory (Le Genest-Saint-Isle, France) and then maintained in pathogen-free facilities for experiments and for the specified analyses. Experimental protocols were approved by the Regional Ethic Committee of Toulouse Biological Research Federation (ref C2EA-01, FRBT) and by the French Ministry for Higher Education and Research. The European directive 2010/63/EU was followed for guidelines on animal welfare.

Cells and flow cytometry

Healthy donors PBMC were obtained after Ficoll-Hypaque (GE Healthcare) density centrifugation and cultured in RPMI 1640 medium (Dutscher) supplemented with 10% fetal bovine serum (ThermoFisher Scientific). Human tumor biopsy specimens were weighed, washed with PBS and mechanically disrupted. These preparations were centrifuged on a Ficoll gradient and frozen. When specified in text, T cells were isolated by negative selection using the Human T Lymphocytes Enrichment Set-DM (BD Biosciences) according to the manufacturer’s instructions. TC1 cells expressing the HPV 16 E6 and E7 proteins were developed in the laboratory of T.C. Wu (Johns Hopkins University). Resting PBMC control was processed for immunostaining with the specified antibodies or ADT (Biolegend) for 20 min at room temperature. For intracellular stainings, cells were fixed and permeabilized using the transcription factor staining buffer set (eBioscience) before staining with the proper antibodies or reagent. Used antibodies are listed on the Key Resource Table. Flow cytometry acquisitions were performed on a BD LSR II cytometer (BD Biosciences) and data were analyzed with Cytobank (http://www.cytobank.org).

Single cell RNA sequencing

The Chromium Controller Single-Cell Instrument and Chromium Single Cell 3′ Library & Gel Bead Kit (V3) were used to prepare individually barcoded single-cell RNA-Seq libraries, following manufacturer’s instruction (10x Genomics). Briefly, cells were thawed and viability was improved by using Dead Cell Removal Kit (Miltenyi Biotec). Then, cells were washed twice with calcium and magnesium-free phosphate buffered saline containing 0.04% bovine serum albumin, and re-suspended at a concentration comprised between 700 and 1200 cells/µL (viability>90%) in the same buffer. In the experiment, 10 000 cells of each sample were loaded to obtain a yield comprised between 2000 and 7000 sequenced cells. Reverse transcription was performed on the C1000 Touch Thermal Cycler, and cDNA was amplified for 12 cycles. cDNA quality control and dosage were performed on Fragment Analyzer, and then by PCR. Double-sided size selection of PCR products was performed with SPRIselect Reagent Kit (Beckman Coulter). The sequencing-ready library was cleaned up with SPRIselect beads (Beckman Coulter) and the sequencing was performed on a HiSeq 3000 instrument (Illumina) with a sequencing depth of at least 50,000 reads per cell.

scRNA-seq data pre-processing and QC

The sequencing output bcl2 files were converted to FASTQ format by using cellranger-mkfastq™ algorithm (10x Genomics), and cellranger-count was applied to align to the GRCh38 reference transcriptome and build the (cell, UMI) expression matrix for each sample. Preprocessing, normalization of UMI, and QC were assessed using the R package Seurat.[29,48] First, genes that correspond to human leukocyte antigens (HLA) and genes located in chromosome X and Y were removed, and the quality of cells was then determined by the number of detected genes per cell, and the proportion of mitochondrial genes. Additional datasets downloaded from 10x Genomics website, Gene Expression Omnibus (GEO), European Genome-phenome Archive (EGA), and ArrayExpress followed the same preprocessing and QC analyses.

Gene annotations

All gene symbols and protein annotations were from public databases (Uniprot, Genecards, and The Human Protein Atlas for humans) and transcription factors genes were defined as such from the GO annotations.

Integration of PBMC scRNAseq datasets

Here, seven different scRNAseq datasets were assembled to digitally extract control γδ T lymphocytes from PBMC of several healthy individuals. This comprised (n = 8,202) pure γδ T lymphocytes cell-sorted by cell surface TCRVδ1 or TCRVδ2 from the PBMC of n = 3 healthy donors[10] and (n = 25,736) PBMC from four other healthy donors (4k +8k from 10XGenomics-3ʹplatform V2, and 10k +1k from 10XGenomics-3ʹplatform V3 downloaded from 10xGenomics resource). After pre-processing and QC discarding of cell doublets and dying cells, the seven datasets were integrated using Seurat 3.0.[29] Briefly, the datasets were normalized together and variable features from each dataset were identified by FinVariable-Data before identification of anchors correspondences across all the datasets, to yield an integrated PBMC dataset containing 26,313 features across n = 25,252 control cells (n = 17,803 PBMC plus n = 7,449 purified γδ T lymphocytes). Principal component analysis (PCA) was performed on this integrated dataset and used to calculate uniform manifold approximation and projection (UMAP), as well as t-distributed Stochastic Neighbor Embedding (t-SNE)[49] coordinates. When visualized in UMAP space,[50] its main cell types were annotated by expression of the CD14, CD19, CD3D, and NCAM marker genes (Supplemental Data Figure 1a).

Single cell signatures and scores

Integrated and normalized scRNAseq datasets were then scored for the multigene signatures specified in text (Table S1) which were downloaded from the MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp) unless specified otherwise. Each single cell score for each geneset was computed as described.[51] Briefly, Single-Cell Signature Scorer computes the score of cell C for gene set GS as the sum of all UMI for all the GS genes expressed by C divided by the sum of all UMI expressed by Cj:

Score of cell C for geneset GS = (∑(UMI)GS)/(∑(UMI)C)[1]

For further visualization, all cell signatures scores are merged with each cell’s UMAP coordinates using Single-Cell Signature Merger.[51]

Score and Gate of γδ T and T CD8 cells from PBMC

Based on using single-cell signature scores to identify cells, this procedure allows the selective and exhaustive digital isolation from any scRNAseq dataset of any subset of cells specified by characteristic single-cell scores. The γδ T cell detection, which previously relied on a single-gene signature,[10] was improved here by implementing a sequence of positive and negative cell selections from several gene signatures, called ‘Score and Gate.’ The above (n = 25,252) PBMC-extracted cells were scored for the specified gene signatures (Table S1) using Single-Cell Signature Explorer[51] and scores were sequentially visualized on scatterplots by Single-Cell Virtual Cytometer.[19] The Score and Gate procedure comprises five sequential gating steps, here applied starting from the (n = 17,803) PBMC sub-dataset alone:

1. Positive selection of double negative cells [= T and NK lineages (n = 11,768)] from the scatterplot of ‘B cell’ and ‘Myeloid cell’ signatures.

2. Scatterplot of ‘B cell’ signature against the addition of ‘CD8AB’ and ‘TCRγδ’ signatures to positively select the double positive cells [= non-CD4 T cells (n = 6,162)].

3. Scatterplot of ‘CD8ab’ signature against ‘TCRα constant gene TRAC’ signature to digitally extract double-positive cells [= αβ T CD8 lymphocytes (n = 3,680)], on the one hand, and negative cells [= γδ T (n = 2,161) plus CD4CD8-double negative αβ T cells], on the other hand.

4. Scatterplot of the addition of ‘CD3 complex’ and ‘TCRα constant gene TRAC’ signatures against the addition of ‘CD3 complex’ and ‘TCRγδ’ signatures to extract the ‘CD3 complex’-‘TCRγδ’ double-positive cells [=γδ T lymphocytes (n = 580, 3% of initial PBMC)].

TCRVγ9 and TCRVγnon9 cell subtyping by compensated scores

Since γδ T lymphocytes express either the TRGC1-encoded TCRVγ9 or the TRGC2-encoded TCRVγnon9 in a mutually exclusive fashion, the extracted γδ T lymphocytes were classified on ‘TRDC, TRGC1ʹ positive cells [=TCRVγ9 cells (n = 312)] or ‘TRDC, TRGC2ʹ positive cells [=TCRVγ9 cells (n = 268)] using a compensated version of these two signatures. Compensated scores are obtained by multiplying the score of the geneset of interest, here GSSTCRV, by its difference to its complementary geneset (GSSTCRV minus GSSTCRV):

GSSTCRVcomp = GSSTCRV * (GSSTCRV − GSSTCRV)[1]

Scatterplots of compensated scores for the two complementary genesets GSS and GSS unambiguously identify the TCRVγ9 or TCRVγnon9 subset of each γδ T lymphocyte (Supplemental Data Figure 1b).

Validation of the Score and Gate method of digital extraction

The method was validated on two different scRNAseq datasets: our formely depicted dataset of purified γδ T lymphocytes from three healthy donors[10] and a CITEseq dataset of PBMC from another healthy donor.[19] When applied to the (n = 7,449) purified γδ T cell sub-dataset, this same score and gate strategy extracted (n = 5,542) γδ T cells (74% sensitivity). Given their already known TCRVδ1/TCRVδ2 subtype, the TCRVγ9/non9 subdivision delineated (n = 2,328) TCRVγ9Vδ2, (n = 1,993) TCRVγnon9Vδ1, (n = 817) TCRVγnon9Vδ2, and (n = 404) TCRVγ9Vδ1 cells. Notwithstanding the equal rates of annotated TCRVδ1 and TCRVδ2 cells deliberately composing this sub-dataset,[10] such results confirmed the prominence of TCRVγ9Vδ2 cells over all the other TCR subsets of blood γδ T, and the preferential pairings of TCRVγ9 (85%) with TCRVδ2, and of TCRVγnon9 (71%) with TCRVγ1 (Supplemental Data Figure 1c). In the UMAP, (n = 5,537 out of the 5,542) extracted γδ T cells are embedded in the T cell clusters (Supplemental Data Figure 1d) showing the γδ T cell detection specificity by score and gate had improved (99.9 vs 95%) over earlier method,[10] and successfully extracted those γδ T cells embedded in CD4 and CD8 T cells clusters. Since CITEseq datasets provide both transcriptomes and immune phenotypes of single cells,[52] the same score and gate sequence was also applied to our CITEseq dataset of 3k PBMC from a healthy donor (GSE144434).[19] This extracted (n = 137) γδ T cells (4.7% of PBMC, encompassing n = 81 TCRVγ9 plus n = 56 TCRVγnon9 cells) of cell surface CD3+CD4−CD8− phenotype, and embedded with T CD8 and NK cells in UMAP (Supplemental Data Figure 7), as expected. These results validated the score and gate method, and the lymphocytes extracted from the CITEseq dataset were integrated into the control datasets. The integrated (cell x gene) matrices of PBMC-derived γδ T and T CD8 single cells are provided in Tables S3, S4, respectively, together with their annotations for pseudotime and differentiation hallmarks in Tables S5, S6.

Extraction of TILs from scRNAseq datasets of tumor samples

The scRNA-seq of tumor lesions from three patients with splenic marginal zone lymphoma (SMZL) were produced for this study, and additional datasets of tumors from other cancers were downloaded from public resources (GEO, EGA, and AE). These encompassed tumor samples from cutaneous T cell lymphoma (CTCL),[53] Hodgkin lymphoma (HL),[54] lung adenocarcinoma (LUAD)[55,56], breast carcinoma (BRCA),[57] triple negative breast cancers (TNBC),[32] head and neck squamous cell carcinoma (HNSCC),[58] hepatocellular carcinoma (HCC),[59] basal cell carcinoma (BCC),[22] and melanoma (MEL)[24] (Table S2). The same score and gate procedure was applied to the pre-processed and quality controlled datasets for digital extraction of (n = 4,648) γδ T and (n = 45,174) T CD8 cells, referred to as TILs. From the CTCL samples, the malignant T lymphocytes were filtered out from TILs by gating out those cells overexpressing the marker gene TOX and the CTCL-specific 17 gene signature.[53] The integrated (cell and gene) matrices of both PBMC- and TILs-derived γδ T and T CD8 single cells are, respectively, provided in Tables S3, S4 with their respective annotations in Tables S5, S6.

Trajectory Inference After the digital extraction of γδ T and T CD8 lymphocytes from normal and cancer samples, both control and TILs were integrated, independently for each lineage (γδ T and T CD8), by using Seurat (R package),[29] as described above. Next, the γδ T and T CD8 integrated matrices were independently used for maturation trajectory inference using dynverse (collection of R packages), which provides more than 45 different methods for trajectory inference (TI). The TI was calculated following the recommended workflow.[60] Briefly, the normalized (log2) of raw expression data, and UMAP coordinates were used as an input. The most suitable and rapid method for our datasets determined by applying ‘guidelines’ function (dynverse on R) was Minimum Spanning Tree (MST) to infer the maturation trajectories of the γδ T and T CD8, independently (‘ti_mst’ function, dynverse on R), with default settings. Additional signatures (‘Tn’ and ‘Teff’) were scored, as described above, and together with the obtained MST attributes: dimension reduction coordinates 1 and 2 (MST1 and MST2), milestones and pseudotime, were used to define the root of the trajectory (‘add_root’ function) of the γδ T and T CD8 lymphocytes. The MST attributes allow multiple visualizations, so the data were visualized through pseudotimed trajectories.

Pseudotimed trajectories

The pseudotimed trajectory is a bidimensional representation of the trajectory inferred for any single cell (p) in which (x,y,z) are the MST coordinates for pseudotime, MST1, and MST2, respectively. The dim_proj method takes each cell’s (x,y,z) from this 3-dimension space (X: pseudotime, Y: MST1, Z: MST2) and projects them in a 2-dimensional plane where pseudotime is one of these dimensions (XoY’). Between the three orthonormal dimensions (XYZ) and the projection plane (XoY’), there is a rotation angle (ZoY’) called β (Supplemental Data Figure 8a). The new projected coordinates (x,y’) of the cell p are calculated as follows:

o’p = y2+ z2 (1)

cos(α+β) = o’y’/o’p (2)

So the new Y coordinate y’ of p projection in the plane is:

o’y’ = * cos(α+β) (3)

o’y’ = * cos(arctan()+β) (4)

To avoid negative values, the whole trajectory was translated by adding the minimal coordinates (y) of the trajectory cells to the initial (y, z) coordinates:

o’y’ = * cos(arctan(()+α) (5)

The resulting trajectory using only cell’s coordinates (x, y’) is called pseudotimed trajectory. Of note, although the trajectory itself is independent of the rotation angle β, selecting this latter allows us to optimize visualization of the same trajectory and its eventual branches (Supplemental Data Figure 8b).

Heatmaps

Heatmaps of gene expression level (UMI mean) per milestone were generated with color scale based on the z-score distribution per row, and shown from row min to row max. Heatmaps of multigene signature score per milestone (milestone mean) were generated by computing scores with Single-Cell Signature Explorer[51] and shown by color scale from the z-score distribution per row. Heatmaps of cell counts per milestone were generated as fractions of the row total.

Classifications

The cell classifications Ttrm and Tex were performed by ‘at least one binary’ as follows. The specified single cells were scored for several reference gene signatures (Table S1). For Tex classification, the reference Tex signatures were five published and partially overlapping Tex genesets.[16,23,61-63] For each signature, a cutoff was defined as the maximal score of the (n = 3,680) control γδ T CD8 lymphocytes extracted from healthy adults PBMC. This threshold defined the cell’s binary (1: cell score > threshold, 0: otherwise). The Tex cutoffs were: 3.9 for ‘Chihara_IL27_Coinhib_module’;[61] 0.52 for ‘Alfei_d20_tox’;[62] 0.22 for ‘Khan_Tox_OverExpressed_genes’;[16] 0.5 for ‘Tosolini_ NHL_IEGS33ʹ;[63] and 0.16 for ‘Balanca_QP_genes’.[23] For each TIL, the five Tex binaries were summed, and any TIL with a non-zero sum of binaries was classified Tex. The Ttrm classification was applied likewise using six reference and partially overlapping Ttrm signatures from the literature.[20,64] The binarizing cutoffs established with the control T CD8 lymphocytes were: 0.46 for ‘Kumar_13g_Ttrm’; 0.3 for ‘Kumar_3g_Ttrm’;[64] 3.72 for ‘Wu_Tcellcluster4.1_trm’; 2.6 for ‘Wu_Tcellcluster8.3_trm’; 2.8 for ‘Wu_Tcellcluster8.3b_trm’; 4.2 for ‘Wu_Tcellcluster8.3c_trm’.[20] These Ttrm binaries were summed, and any TIL with a non-zero sum of Ttrm binaries was classified Ttrm. The cytolytic T cell infiltrates of tumors were arbitrarily classified by total cytolytic (T CD8 + γδ T) TIL counts as deserts (D, <100 TILs), classic (C, 100–1000 TILs), or inflamed (I, >1000 TILs). The Tex classification from the flow cytometry data was based on both median (PD-1hiTIM-3+)>8% and median (PD-1hiCD39+)>12%. The reactivity of tumor infiltrates (TI) from ICB-treated MEL and BCC patients was classified as TI-R (reactive)/TI-NR (non-reactive) based on the difference of total TIL counts (post–pre-ICB therapy) by TI-R: ∆>0, else TI-NR.

Data visualization

Single-cell Multilayer Viewer, was developed as a serverless software for the visualization of up to 5 layers of information. It is available on GitHub repository: https://github.com/MarionPerrier/scMLV. Single-cell Multilayer Viewer allows us to merge and visualize up to 5 layers of quantitative and qualitative variables, alone and in combination.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.