# Single-Cell RNAseq Profiling of Human γδ T Lymphocytes in Virus-Related Cancers and COVID-19 Disease

Juan Pablo Cerapio, Marion Perrier, Frédéric Pont, Marie Tosolini, Camille

Laurent, Stephane Bertani, Jean-Jacques Fournie

► To cite this version:

Juan Pablo Cerapio, Marion Perrier, Frédéric Pont, Marie Tosolini, Camille Laurent, et al., Single-Cell RNAseq Profiling of Human γδ T Lymphocytes in Virus-Related Cancers and COVID-19 Disease. Viruses, 2021, 13 (11), pp.2212. 10.3390/v13112212. [hal-03413561](https://hal.science/hal-03413561)

HAL Id: hal-03413561

https://hal.science/hal-03413561

Submitted on 3 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Single-Cell RNAseq Profiling of Human γδ T Lymphocytes in Virus-Related Cancers and COVID-19 Disease

Juan Pablo Cerapio<sup>1,2,3,4,5,*</sup>, Marion Perrier<sup>1,2,3,4,5</sup>, Frédéric Pont<sup>1,2,3,4,5 D</sup>, Marie Tosolini<sup>1,2,3,4,5 D</sup>, Camille Laurent<sup>1,2,3,4,5,6,7</sup>, Stéphane Bertani<sup>2,8 D</sup> and Jean-Jacques Fournie<sup>1,2,3,4,5,6,*,†</sup>

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Cerapio, J.P.; Perrier, M.; Pont, F.; Tosolini, M.; Laurent, C.; Bertani, S.; Fournie, J.-J. Single-Cell RNAseq Profiling of Human γδ T Lymphocytes in Virus-Related Cancers and COVID-19 Disease. *Viruses* **2021**, *13*, 2212. https://doi.org/10.3390/v13112212

Academic Editor: Eric Champagne

Received: 14 September 2021

Accepted: 20 October 2021

Published: 3 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** The detailed characterization of human γδ T lymphocyte differentiation at the single-cell transcriptomic (scRNAseq) level in tumors and patients with coronavirus disease 2019 (COVID-19) requires both a reference differentiation trajectory of γδ T cells and a robust mapping method for additional γδ T lymphocytes. Here, we incepted such a method to characterize thousands of γδ T lymphocytes from (n = 95) patients with cancer or adult and pediatric COVID-19 disease. We found that cancer patients with human papillomavirus-positive head and neck squamous cell carcinoma and Epstein-Barr virus-positive Hodgkin's lymphoma have γδ tumor-infiltrating T lymphocytes that are more prone to recirculate from the tumor and avoid exhaustion. In COVID-19, both TCRVγ9 and TCRVγnon9 subsets of γδ T lymphocytes relocalize from peripheral blood mononuclear cells (PBMC) to the infected lung tissue, where their advanced differentiation, tissue residency, and exhaustion reflect T cell activation. Although severe COVID-19 disease increases both recruitment and exhaustion of γδ T lymphocytes in infected lung lesions but not blood, the anti-IL6R therapy with Tocilizumab promotes γδ T lymphocyte differentiation in patients with COVID-19. PBMC from pediatric patients with acute COVID-19 disease display similar γδ T cell lymphopenia to that seen in adult patients. However, blood γδ T cells from children with the COVID-19-related multisystem inflammatory syndrome are not lymphodepleted, but they are differentiated as in healthy PBMC. These findings suggest that some virus-induced memory γδ T lymphocytes durably persist in the blood of adults and could subsequently infiltrate and recirculate in tumors.

**Keywords:** human; gammadelta; lymphocyte; tumor; COVID-19; transcriptome; single cell; differentiation; trajectory

## 1. Introduction

Efficient protection against viruses and cancers requires innate and adaptive immunity. Due to their ontogeny, antigen recognition mode, and effector functions, γδ T lymphocytes are involved in both types of immune defense [1]. As a result, γδ T cells are activated by a wide array of viral infections as well as by many different types of cancer. Human CD4⁻CD8⁻ γδ T lymphocytes express either TCRVδ2 or TCRVδ1/3 chains prominently paired with TCRVγ9 or TCRVγ2/3/4/5/8 (TCRVγnon9), respectively. The TCRVγ9Vδ2 γδ T cells are prominently circulating lymphocytes specifically reactive to phosphoantigen-sensing butyrophilins while TCRVγnon9 γδ T cells are mostly tissular and recognize

<footnote n="1">¹ Centre de Recherches en Cancérologie de Toulouse, Inserm UMR 1037, 31037 Toulouse, France; marion.perrier@inserm.fr (M.P.); frederic.pont@inserm.fr (F.P.); marie.tosolini@inserm.fr (M.T.); camille.laurent@inserm.fr (C.L.)</footnote>

<footnote n="2">² Toulouse University, 31062 Toulouse, France; stephane.bertani@ird.fr</footnote>

<footnote n="3">³ CNRS UMR 5071, 31024 Toulouse, France</footnote>

<footnote n="4">⁴ Institut Universitaire du Cancer-Oncopole de Toulouse, 31100 Toulouse, France</footnote>

<footnote n="5">⁵ Laboratoire d'Excellence 'TOUCAN-2', 31037 Toulouse, France</footnote>

<footnote n="6">⁶ Institut Carnot Lymphome CALYM, 75014 Paris, France</footnote>

<footnote n="7">⁷ Centre Hospitalier Universitaire, 31300 Toulouse, France</footnote>

<footnote n="8">⁸ PHARMADEV, IRD UMR 152, 31062 Toulouse, France</footnote>

<footnote>* Correspondence: juan-pablo.cerapio-arroyo@inserm.fr (J.P.C.); jean-jacques.fournie@inserm.fr (J.-J.F.)</footnote>

<footnote>† Lead contact.</footnote>
more diversified antigens. Upon antigen activation, both subtypes differentiate from naive (Tn) to central memory (Tcm), effector memory (Tem), and terminally differentiated effector memory CD45Ra+ (Temra) cells [2]. Regardless of the γδ TCR subset, the strong cytotoxic properties of γδ Tem cells play a crucial role in the killing of virus-infected and cancer cells. For decades, TCRVγ9Vδ2 γδ T cells have been known to react to cells infected by viruses such as the Epstein–Barr virus (EBV) [3], herpes simplex viruses [4], human immunodeficiency virus (HIV) [5], [6], human papillomavirus (HPV) [7], hepatitis B virus (HBV) [8], Influenza viruses [9,10], dengue virus [11], and severe acute respiratory syndrome-related coronavirus (SARS-CoV) [12]; whereas the TCRVγnon9 lymphocytes [13] respond, inter alia, to cytomegalovirus (CMV), hepatitis E virus (HEV) [14], and Kaposi's sarcoma-associated herpesvirus (HHV-8) [15].

The overall T and γδ T cell response remain hitherto poorly understood in coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection. Like other respiratory viral infections, COVID-19 can exhibit distinct patterns of clinical features and disease severity, ranging from a virus-induced respiratory pathology due to a weak immune response to an overaggressive immunopathology caused by excessive immune activation. The heterogeneous immune profiles recently characterized in peripheral blood mononuclear cells (PBMC) of COVID-19 patients reflect this wide range of symptoms. COVID-19 patients have a quantitatively decreased compartment of γδ T cells that are more differentiated than control samples [16], accompanied by an increased number of immature neutrophils [17]. Similar alterations were reported with the other adaptive and unconventional T lymphocytes [18,19], supporting a scenario of broad T cell immunosuppression linked to fatal outcomes [20].

Likewise, the extent and quality of tumor-infiltrating T lymphocytes (TILs) in cancer are important determinants of the outcome of immunotherapy or untreated cancer patients. On the one hand, a weak immune response to cold tumors, peripherally excluding tumors, or tumors infiltrated by exhausted T lymphocytes allows for uncontrolled cancer and leads to fatal outcomes. On the other hand, too strong immune reactions are sometimes observed in untreated cancers (for review [21]) and systematically in patients receiving immune checkpoint blockade. Therefore, viral infections and cancer can present large and heterogeneous spectrums of immune response patterns. In human cancer, our research and other studies have shown that the abundance of γδ TILs, whether TCRVγ9 or TCRVγnon9, varies substantially across individuals and among cancer types, and that it is associated with patient outcomes [22,23]. As opposed to peripheral blood and hematopoietic malignancies dominated by TCRVγ9 γδ T cells, the predominant subset of γδ TILs in most solid tumors are TCRVγnon9 cells, whose number, cytotoxic differentiation, and functional status vary between individuals [24–26].

We recently depicted the pan-cancer landscape of human γδ TILs using single-cell RNA sequencing (scRNAseq). This strategy unveiled that the differentiation of cytotoxic γδ T cells from the TCRVγ9 and TCRVγnon9 subsets was related to the CMV status of the healthy donors [27]. Furthermore, despite the considerable heterogeneity of cancer patients, their infiltration and differentiation appeared phased with that of the adaptive T CD8 TIL compartment, yielding correlated rates of functional and/or exhausted γδ T and T CD8 TILs in human cancers [28].

Overall, these findings raised questions as to whether the viral status of cancer patients affected their γδ TILs, and whether the transcriptomic profiles correlated with those of T lymphocytes from COVID-19 patients' PBMC. Here, we addressed these questions through a comprehensive characterization of the differentiation and functional hallmarks of γδ T lymphocytes by analyzing scRNAseq datasets of tumors from cancer patients with known viral status, and PBMC from COVID-19 patients.
## 2. Materials and Methods

### 2.1. scRNAseq Datasets Pre-Processing

The (cells, gene) matrix data for the control γδ T isolated from healthy donor's PBMC has been previously published [28]. The scRNAseq of tumor lesions from Hodgkin's lymphoma (HL) (EGAS00001004085) [29], head and neck squamous cell carcinoma (HNSCC) (GSE139324) [30], PBMC from COVID-19 patients (GSE145926, GSE155224, GSE155249, GSE166489, and GSE167029) were downloaded and assembled to digitally extract their γδ T lymphocytes. After pre-processing and discarding cell doublets and dying cells, all datasets were integrated using the R package Seurat version 4.0 [31]. Principal components analysis (PCA) was performed on this integrated dataset before uniform manifold approximation and projection (UMAP) [32] to check the integrated dataset quality.

### 2.2. Single-Cell Signatures and Scores

After normalization and integration, the integrated scRNAseq dataset was scored for a collection of multigene signatures [28] using Single-Cell Signature Explorer software [33]. Briefly, the score of each single-cell Cj for gene set GSx was computed as the sum of all unique molecular identifiers (UMI) for all the GSx genes expressed by Cj divided by the sum of all UMI expressed by Cj:

$$ \text{Score of cell Cj for geneset GSx} = \frac{(\sum(\text{UMI})\text{GSx})Cj}{(\sum(\text{UMI})Cj)} \quad (1) $$

All cell signature scores were merged with each cell's trajectory coordinates using Single-Cell Signature Merger software [33]. For visualization, Single-Cell Multilayer Viewer was used, a serverless software allowing to merge up to five layers of quantitative and qualitative variables, alone or in combination [28].

### 2.3. Score and Gate for Digital Extraction of γδ T Lymphocytes

Starting from any scRNAseq datasets pre-processed as above, this procedure comprised five gating steps sequentially applied using Single-Cell Virtual Cytometer software [34]:

1. Positive selection of double negative cells (= T and NK lineages) from the scatterplot of 'B cell' and 'Myeloid cell' signature scores.

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).

3. Scatterplot of 'CD8ab' signature against 'TCRα constant gene TRAC' signature to digitally extract double-positive cells (= αβ T CD8 lymphocytes) on the one hand, and negative cells (= γδ T 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' and 'TCRγδ' double-positive cells (=γδ T lymphocytes).

The datasets of the finally identified γδ T cells were digitally extracted from the original scRNAseq dataset and further subtyped and injected on the trajectory (see below).

### 2.4. Subtyping of the TCRVγ9 and TCRVγnon9 γδ T Cells

This was performed as previously described [28]. Briefly, γδ T lymphocytes express either the TRGC1-encoded TCRVγ9 or the TRGC2-encoded TCRVγnon9 in a mutually exclusive fashion, so the extracted γδ T lymphocytes were categorized as TCRVγ9 cells based on ('TRDC, TRGC1' positive cells) or TCRVγnon9 cells based on ('TRDC, TRGC2' positive cells) classifications using a compensated score of these two signatures. Compensated scores were obtained by multiplying the score of the gene set of interest, here GSSTCRVγ9, by its difference to its complementary gene set (GSSTCRVγ9 minus GSSTCRVγnon9):

$$ \text{GSS}_{\text{TCRV}\gamma9\text{comp}} = \text{GSS}_{\text{TCRV}\gamma9} \times (\text{GSS}_{\text{TCRV}\gamma9} - \text{GSS}_{\text{TCRV}\gamma\text{non9}}) \quad (2) $$
Scatterplots of GSS<sub>TCRVγ9</sub> and GSS<sub>TCRVγnon9</sub> finally identified the TCRVγ9 or TCRVγnon9 subset of each γδ T lymphocyte [28].

## 2.5. Cells Injection on Public Pseudotimed Differentiation Trajectory

The maturation trajectory of γδ T cells was computed by minimum spanning tree and visualized as pseudotimed trajectories. Pseudotimed trajectories are differentiation trajectories shown with cell pseudotime on the x-axis and a projection of both trajectory dimensions (here MST1 and MST2) on the y-axis [28]. We will refer below to them as 'public γδ T lymphocytes trajectory' the pseudotimed differentiation trajectory of all γδ T lymphocytes from ~150 tissues samples, including healthy donors' PBMC and cancer patients' tumors (the (cell, gene) matrix data is available from [28]). New γδ T lymphocytes digitally extracted from additional datasets, such as those from COVID-19 patients, were injected onto this 'public γδ T cell trajectory'. Briefly, a total of (n = 9) datasets for COVID-19-derived tissues were downloaded. These COVID-19 datasets were first filtered for quality control, integrated as described above, and their γδ T lymphocytes were digitally extracted by score and gate [28]. These new γδ T lymphocytes were then injected onto the public trajectory by using Seurat's Integration and Map-to-Query workflow to determine their respective (MST1, MST2) coordinates and then infer their pseudotime within the public trajectory. This procedure encompasses the following four steps:

Step 1:

* Query 1 = γδ T lymphocytes digitally extracted from COVID-19 samples;

* Reference = public trajectory;

* Integration of gene expression data from Query 1 and Reference (Seurat's Integration).

Step 2:

* Query 2 = Integrated (Query and Reference);

* Reference = public trajectory;

* Transferring MST coordinates from Reference metadata to Query 2 (Seurat's Map-to-Query).

Step 3:

* Query 2 = pseudotime calculation using MST coordinates and R package Dynverse.

Step 4:

* Visualization of new data (Query 1) on the public trajectory (Reference) using Single-Cell Multilayer Viewer [28].

## 2.6. Classifications

The TILs classifications as Ttrm/non-Ttrm and Tex/non-Tex were performed by 'at least one binary'. Briefly, the single cells were scored for several reference gene signatures as previously described [28]. For Tex classification, the reference Tex signatures were based on five published and partially overlapping Tex gene sets [35–39]. For each signature, a cut-off was defined as the maximal score of the (n = 3680) control γδ T lymphocytes extracted PBMC of healthy adults. This threshold defined the cell's binary (1: cell score > threshold, 0: otherwise). The Tex cut-offs were: 3.9 for 'Chihara_IL27_Coinhib_module' [35]; 0.52 for 'Alfei_d20_tox' [36]; 0.22 for 'Khan_Tox_OverExpressed_genes' [37]; 0.5 for 'Tosolini_NHL_IEGS33' [38]; and 0.16 for 'Balanca_QP_genes' [39]. For each TIL, the five Tex bina-ries were summed, and if the sum was non-zero, the TIL was classified as Tex, or otherwise as non-Tex. The Ttrm classification was applied likewise using six references and par-tially overlapping Ttrm signatures published previously [40,41]. The binarizing cutoffs established with the control γδ T lymphocytes were: 0.46 for 'Kumar_13g_Ttrm'; 0.3 for 'Kumar_3g_Ttrm' [40]; 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' [41]. These Ttrm bi-naries were summed, and if the sum was non-zero, the TIL was classified as Ttrm, or otherwise as non-Ttrm (recirculating).
## 3. Results

### 3.1. Tumor Infiltrating γδ T Lymphocytes from Virus-Positive and -Negative Cancer Patients

We downloaded the published scRNAseq dataset from a total of (n = 26) human tumors of head and neck carcinoma (HNSCC), including (n = 8) samples positive and (n = 18) negative for HPV, as well as (n = 9) Hodgkin's lymphoma (HL) of mixed cellularity subtype in which (n = 5) samples were positive and (n = 4) negative for EBV. Their respective γδ T lymphocytes were digitally purified, and their TCR subtype and differentiation stage were determined by injection onto the 'public γδ T cell trajectory' (see Methods) and cross-labeling with reference gene signatures from external single-cell datasets of human TILs [28].

A total of (n = 1055) γδTILs ((n = 401) TCRVγ9 plus (n = 654) TCRVγnon9 cells) were extracted from HNSCC, and (n = 113) γδ TILs ((n = 51) TCRVγ9 plus (n = 62) TCRVγnon9 cells) from HL. In contrast with PBMC from healthy donors, the cell counts extracted from tumors varied considerably across individuals and in both diseases (range: 1–33 γδ TILs in HL patients, and 0–111 in HNSCC patients), as previously reported in other cancer types [28]. On average, HNSCC comprised more γδ TILs (mean = 40 cells, 15 TCRVγ9 + 25 TCRVγnon9) than HL (mean = 13 cells, 6 TCRVγ9 + 7 TCRVγnon9; t-test P = 3 × 10⁻⁴). The rates of their Tn, Tcm, Tem, and Temra stages varied between patients but without reaching significant difference between cancers and viral status (mean in HNSCC: 15% Tn, 37% Tcm, 47% Tem, 0.9% Temra; mean in HL: 9% Tn, 35% Tcm, 55% Tem, 0.8% Temra) (Figure 1a). In both cancers, these recirculating (non-Ttrm TILs, Methods) γδ TILs comprised less exhausted cells than their tissue-resident (Ttrm) counterparts (15% and 37%, respectively). Furthermore, the rates of recirculating γδ TILs were significantly associated with patient's viral status (χ² P value = 7.8 × 10⁻¹⁸): tumors of viral-positive HL and HNSCC included far more recirculating γδ TILs than viral-negative tumors (80% vs. 45% of γδ TILs, respectively) (Figure 1b).

Altogether, these data demonstrate that patients with HPV-positive HNSCC and EBV-positive HL have γδ TILs more prone to recirculate from the tumor and avoid exhaustion.

### 3.2. γδ T Lymphocytes from COVID-19 Patients

The scRNAseq datasets from PBMC of (n = 43) COVID-19 patients and (n = 21) healthy controls, as well as from broncho-alveolar lavage fluids (BALF) from (n = 20) COVID-19 patients and (n = 4) healthy controls were downloaded. A total of (n = 7473) γδ T lymphocytes were digitally extracted from these datasets, and their TCR subtype and differentiation stage were characterized as described above (see Methods).

In spite of similar dataset sizes, each sample of COVID-19 PBMC yielded on average three times less γδ T lymphocytes than the healthy controls (100 vs. 330, respectively, t-test P = 4 × 10⁻⁵), and both TCR subsets had decreased cell counts (on average 63 vs. 142 TCRVγ9 and 37 vs. 187 TCRVγnon9 in COVID-19 vs. control samples, respectively). Parallel analyses showed on average n = 561 and n = 479 T CD8 cells per sample in the same groups indicating that the COVID-19-related lympho-reduction from PBMC is specific to the γδ T cell lineage (Chi-2 P value = 5 × 10⁻²⁷).

However, their differentiation stages remained unchanged between healthy and COVID-19 PBMC (7% Tn, 25% Tcm, 66% Tem, 3% Temra among TCRVγ9 cells; 19% Tn, 30% Tcm, 50% Tem, 2% Temra among TCRVγnon9 cell) (Figure 2a,b). In contrast, BALF of COVID-19 patients encompassed more γδ T lymphocytes than control BALF (on average 16 vs. 0 γδ T per BALF, respectively). Such COVID-19 BALF-derived γδ T encompassed both subsets (on average nine TCRVγ9 and seven TCRVγnon9 cells) and were more differentiated than in PBMC (on average 1% Tn, 18% Tcm, 81% Tem, 0% Temra among TCRVγ9 cells; 3% Tn, 36% Tcm, 61% Tem, 0% Temra among TCRVγnon9). Furthermore, the Ttrm gene signature was expressed by scarce (0.3%) γδ T lymphocytes in the control PBMC, but by 4% in COVID-19 PBMC and 17% in the COVID-19 BALF. Likewise, the Tex signature was up-regulated by 0%, 1%, and 24% γδ T lymphocytes, respectively, from these COVID-19 BALF samples (Figure 2a,b).
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Differentiation and functional hallmarks of γδ T lymphocytes infiltrating HNSCC and HL tumors according to viral status of cancer patients. (a) Examples of γδ TILs extracted from representative HNSCC and HL tumors with viral status, overlaid on the public trajectory of γδ T lymphocytes reference (grey). Key: TCRVγ9 cells (light blue), Tex TCRVγ9 cells (dark blue), TCRVγnon9 cells (orange), Tex TCRVγnon9 cells (dark red), recirculating (non-Ttrm) cells (open circles), Ttrm (cross). Differentiation pseudotime scale: 0–5: Tn, 5–20: Tcm, 20–50: Tem; > 50: Temra. Together all the EBV-negative HL tumors totalized 53 mostly functional (43/53, 81%) γδ TILs, including 13 Ttrm with 39% of exhausted cells (5/13 Ttrm), while all the EBV-positive HL tumors totalized 60 mostly functional (56/60, 93%) γδ TILs, including 2 non-exhausted Ttrm γδ TILs. The HPV-negative HNSCC tumors totalized 698 mostly functional (489/698, 70%) γδ TILs, including 326 Ttrm with 39% of exhausted cells (127/326 Ttrm), while the HPV-positive HNSCC tumors totalized 357 mostly functional (310/357, 87%) γδ TILs, including 81 Ttrm with 31% of exhausted cells (25/81 Ttrm). (b) Right: rates of recirculating and Ttrm among γδ TILs (both subsets) per viral status of patients (pooled); left: rates of exhausted cells among the recirculating and Ttrm γδ TILs.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Differentiation and functional hallmarks of γδ T lymphocytes from control donor's PBMC and BALF. (a) γδ T lymphocytes extracted from PBMC and BALF of healthy individuals and COVID-19 patients. No γδ T lymphocyte was found in any of the healthy control's BALF (representative examples, cells are shown overlaid on the public trajectory of γδ T lymphocyte differentiation, same legend key as in Figure 1). (b) Rates of cells at each differentiation stage among the TCRVγ9 and TCRVγnon9 γδ T lymphocytes (group means are shown).</figcaption></figure>

These data suggest that in adult COVID-19 patients, both subsets of γδ T lymphocytes relocalize from PBMC to the infected lung tissue, where their advanced differentiation, tissue residency, and exhaustion hallmarks witness of a recent T cell activation.
### 3.3. γδ T Lymphocytes in Clinical Groups of Adult Patients with COVID-19

These γδ T cell analyses were then extended to clinical subgroups of adult COVID-19 patients. The total number of γδ T lymphocytes in BALF of adult COVID-19 patients varied considerably between individuals without reaching a statistical difference between mild (n = 2) and severe (n = 6) disease (on average 39 vs. 22 γδ T cells per sample, respectively). The TCRVγ9 cell counts were similar in both groups (on average 13 vs. 14 cells per sample, respectively), whereas the TCRVγnon9 cell counts decreased in severe COVID-19 (on average 26 vs. 7) despite maintaining the same differentiation stage pattern (on average 1% Tn, 37% Tcm, 62% Tem, 0% Temra). Of note, with 0% Tn, 20% Tcm, 80% Tem, and 0% Temra in both groups, the TCRVγ9 cells were more differentiated than the TCRVγnon9 cells. The γδ T lymphocytes in BALF of patients with mild disease included 9% Ttrm and 12% Tex cells, contrasting to 40% Ttrm and 40% Tex in severe disease, and both subsets behaved similarly (Figure 3a). Hence, acute COVID-19 increases both tissue residency and exhaustion of γδ T lymphocytes in infected lung lesions.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Differentiation and functional hallmarks of γδ T lymphocytes from adult COVID-19 patients. (a) BALF-derived γδ T lymphocytes from adult patients with mild or severe COVID-19 disease. (b) PBMC-derived γδ T lymphocytes from adult COVID-19 patients untreated or treated with Tocilizumab. (Representative examples overlaid on the public trajectory of γδ T lymphocyte differentiation, same legend key as in Figure 1).</figcaption></figure>
The PBMC of adult COVID-19 patients treated (n = 12) or not (n = 4) by Tocilizumab, an anti-IL-6R mAb drug, had similar cell counts of total γδ T lymphocytes (on average 59 vs. 45 γδ T cells per sample, respectively), as for both TCRVγ9 and TCRVγnon9 subsets. Meanwhile, Tocilizumab clearly promoted both subsets' differentiation (TCRVγ9 untreated: 8% Tn, 29% Tcm, 61% Tem, 2% Temra; TCRVγ9 treated: 2% Tn, 11% Tcm, 85% Tem, 2% Temra; TCRVγnon9 untreated: 30% Tn, 31% Tcm, 39% Tem, 0% Temra; both subsets treated: 9% Tn, 16% Tcm, 74% Tem, 1% Temra). This treatment did not affect their low rates of Tex cells (Figure 3b). Therefore, anti-IL6R treatment favors γδ T lymphocyte differentiation in COVID-19 patients.

### 3.4. γδ T Lymphocytes in Clinical Groups of Pediatric Patients with COVID-19

Finally, we analyzed the γδ T cell landscape in PBMC of children with acute COVID-19 (n = 5), multisystem inflammatory syndrome (MIS-C) (n = 14), and healthy pediatric controls (n = 7). The PBMC of pediatric controls comprised more γδ T lymphocytes of both TCR subsets than the PBMC of adult controls (n = 4) (average γδ T cell counts: 393 vs. 220, respectively), although their differentiation was similar (TCRVγ9: 5% Tn, 43% Tcm, 51% Tem, 2% Temra; TCRVγnon9: 14% Tn, 45% Tcm, 40% Tem, 1% Temra) (Figure 4a).

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Differentiation and functional hallmarks of γδ T lymphocytes from pediatric COVID-19 patients. (a) PBMC-derived γδ T lymphocytes from adult and pediatric controls. (b) PBMC-derived γδ T lymphocytes from pediatric COVID-19 patients with either MIS-C or acute disease. (Representative examples overlaid on the public trajectory of γδ T lymphocyte differentiation, same legend key as in Figure 1).</figcaption></figure>
PBMC from acute and MIS-C COVID-19 children comprised less γδ T lymphocytes of both subsets than the pediatric controls (average γδ T cell counts: 71 (acute), 163 (MIS-C) vs. 393 (controls)). They were also less differentiated in acute disease than in MIS-C and pediatric controls (acute: TCRVγ9: 14% Tn, 52% Tcm, 34% Tem, 0% Temra; TCRVγnon9: 37% Tn, 44% Tcm, 19% Tem, 0% Temra; MIS-C: TCRVγ9: 10% Tn, 26% Tcm, 61% Tem, 3% Temra; TCRVγnon9: 22% Tn, 28% Tcm, 50% Tem, 1% Temra). Last but not least, the PBMC of pediatric controls contained no Tex cell among γδ T lymphocytes (in both subsets), as in the acute COVID-19 samples, but *n* = 2 Tex out of 2285 cells were detected in MIS-C samples (Figure 4b).

Together, these results indicated that acute COVID-19 disease induces a similar γδ T cell lymphopenia PBMC of pediatric and adult patients. In contrast, the PBMC-derived γδ T cells from MIS-C patients with SARS-CoV-2 infection are not lymphodepleted and resemble to those from healthy PBMC.

## 4. Discussion

In single-cell transcriptomics, pseudotime describes the progression of each cell along- side a continuous series of differentiation states, through a non-linear transformation of real chronological time [42–44]. How explicit time is quantitatively converted into γδ T cell pseudotime remains to be determined, and depends on the transcriptome diversity and the total number of single cells in the dataset. All the features of differentiation trajectories are strictly dependent upon the number of single cells selected. Hence, once a trajectory has been inferred from a given single-cell dataset, adding *ab extra* single cells to this dataset will change the trajectory, preventing their mapping on the former map. Here, we addressed this issue by incepting a method to inject without distortion any amount of *ab extra* cells on the reference trajectory.

In this way, we were able to map and thus characterize thousands of γδ T lymphocytes from newer datasets onto a formerly-built reference map called here 'public γδ T cell trajectory' [28]. Of note, this trajectory encompasses a remarkable set of γδ T lymphocytes, being currently paved with ~30,000 γδ T cells from ~220 healthy and diseased adult individuals. As a result, it achieves a much higher level of exhaustivity and resolution than is possible with the cells from a single individual alone. Together, the public γδ T cell trajectory and the method for mapping newer γδ T cells provide a unique resource for characterizing these lymphocytes from diverse sources at the highest level of resolution.

As a result of high-resolution mapping of thousands of γδ T lymphocytes from patients with cancer or COVID-19, we were able to gain a wealth of information about γδ T lymphocyte differentiation dynamics at the single-cell scale. As reviewed in this series, the role of human γδ T cells in viral infections is highly investigated and debated in the context of SARS-CoV-2 infection [45–48]. Here, we found in PBMC of COVID-19 patients that γδ T lymphocytes were decreased compared to control samples, a finding consistent with flow cytometry studies [16, 20, 49, 50]. Although a peripheral lymphopenia in COVID-19 was initially thought to indicate T cell immunosuppression, our present maps and other studies [18] show that γδ T cells increase concurrently in BALF of COVID-19 patients. T cell exhaustion is rarely observed among the γδ T cells in PBMC from adult and pediatric patients with COVID-19, even during its acute phase as reported with HEV infection [14]. Most Tex γδ T cells were instead found in the infected lung rather than in the circulating blood. Hence, COVID-19-induced γδ T cell exhaustion likely represents a scar of recent localized T cell activation in both TCRVγ9 and TCRVγnon9 subsets of cells rather than in a single subset alone [17]. Moreover, their advanced differentiation stages in blood from children with multisystem inflammatory syndrome associated with SARS-CoV-2 infection further support this view [51]. Tocilizumab, an IL-6 pathway inhibitor, effectively reduces mortality of severely ill and lymphopenic COVID-19 patients, notably by attenuating inflammation and normalizing their circulating T cell rates [52, 53]. Consistent with these reports, the increased effector memory cells we observed here indicates that Tocilizumab
also promotes γδ T cell differentiation, which constitutes a bioactivity of dual interest for both antiviral [45] and anticancer [54] immunity.

Single-cell studies have shown that the differentiation dynamics of human γδ T lym-
phocytes in both TCRVγ9 and TCRVγnon9 subsets primarily reflect the transcriptional
emergence of their cytotoxic function, which peaks at the Tem and Temra stages [27]. In
many solid and hematological human cancers, both lineages of γδ T and T CD8 TILs
present strikingly coherent differentiation profiles, and a sizeable fraction of γδ T TILs are
tissue-resident memory and exhausted cells [28]. Many T CD8 TILs in melanoma are by-
stander cytotoxic lymphocytes with a TCR specificity for viral antigens (V-spe), while only a
minority represent true tumor-antigen-specific (T-spe) lymphocytes [55]. Among such CD8
TILs, the V-spe clonotypes prominently display functional Tcm and Tem differentiation
stages whereas the majority of T-spe clonotypes are dysfunctional Tex cells [56]. In parallel,
we found that HPV-positive HNSCC and EBV-positive HL patients have γδ TILs that are
more prone to recirculate from the tumor and avoid exhaustion. Considering the above
findings, it is tempting to hypothesize that a contingent of peripheral γδ T lymphocytes
induced by viral infections does persist in the long term in the blood of adult individuals as
central/effector memory cells, and can readily infiltrate tumors to mediate a large spectrum
of HLA-unrestricted cytotoxic activities despite true antigen specificity.

For its clinical relevance to cancer immunotherapy, future studies from our team will seek to discern at the single-cell level the T-spe and V-spe clonotypes persisting among the human γδ T lymphocytes in the peripheral blood and infiltrating tumors. Further, our study provides a fundamental framework for characterizing γδ T lymphocytes in viral diseases and cancers, and provides the community with a unique resource for the monitoring of human γδ T lymphocytes at high-resolution.

**Author Contributions:** J.-J.F. and J.P.C. conceived and designed the study; J.P.C. conceived methods, processed, and performed bioinformatic data analyses; J.-J.F. and J.P.C. generated figures, and wrote the manuscript; M.P., F.P. and M.T. conceived algorithms; C.L., S.B. and J.-J.F. supervised the study; all authors revised the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by institutional grants from the Institut National de la Santé et de la Recherche Médicale (INSERM), the Université Toulouse 3-Paul Sabatier (UPS), the Centre National de la Recherche Scientifique (CNRS); Laboratoire d'Excellence Toulouse Cancer (TOUCAN-2) under grant contract ANR11-LABX; the Fondation ARC under grant contract PGA1-RF2019-0208691; the Institut Universitaire du Cancer de Toulouse-Oncopole under contract CIEL; the Institut Carnot Lymphome under contract CALYM. This work was granted access to the HPC resources of CALMIP supercomputing center under allocation 2020-T19001.

**Institutional Review Board:** Not applicable.

**Informed Consent Statement:** Not applicable, the presented study did not involve novel tissue samples from humans or animals.

**Data Availability Statement:** Further information and reasonable requests for resources and pro-
cessed data, such as integrated (cell, UMI) matrices of extracted γδ T and their corresponding cell
annotations, should be directed to J.-J.F. (jean-jacques.fournie@inserm.fr) or J.P.C. (juan-pablo.cerapio-
arroyo@inserm.fr). Source scRNAseq datasets used in this study are available at GEO with the ac-
cession numbers: GSE139324, GSE145926, GSE155224, GSE155249, GSE166489, and GSE167029; and
at the European Genome-Phenome Archive (EGA) with the accession number: EGAS00001004085.
Single-Cell Multilayer Viewer (scMLV) is available on GitHub repository: https://github.com/
MarionPerrier/scMLV (accessed on 31 March 2021). Single-Cell Signature Scorer and Single-Cell
Virtual Cytometer are available at: https://sites.google.com/site/fredsoftwares/products (accessed
on 12 October 2021). Seurat v4 is available from the Comprehensive R Archive Network (CRAN)
and further details in the installation can be found on https://satijalab.org/seurat/ (accessed on
12 October 2021).

**Acknowledgments:** The authors thank CRCT collaborators for their support and their critical discussions.

**Conflicts of Interest:** The authors declare no relevant conflict of interest.
## Abbreviations

The following abbreviations are used in this manuscript

<table><tr><td>CMV</td><td>cytomegalovirus</td></tr><tr><td>COVID-19</td><td>coronavirus disease 2019</td></tr><tr><td>EBV</td><td>Epstein-Barr virus</td></tr><tr><td>HEV</td><td>hepatitis E virus</td></tr><tr><td>HL</td><td>Hodgkin's lymphoma</td></tr><tr><td>HNSCC</td><td>head and neck squamous cell carcinoma</td></tr><tr><td>HPV</td><td>human papillomavirus</td></tr><tr><td>MIS-C</td><td>multisystem inflammatory syndrome in children</td></tr><tr><td>PBMC</td><td>peripheral blood mononuclear cells</td></tr><tr><td>SARS-CoV</td><td>severe acute respiratory syndrome-related coronavirus</td></tr><tr><td>scRNAseq</td><td>single-cell RNA sequencing</td></tr><tr><td>Tcm</td><td>central memory T</td></tr><tr><td>Tem</td><td>effector memory T</td></tr><tr><td>Temra</td><td>effector memory CD45RA T</td></tr><tr><td>Tex</td><td>exhausted T</td></tr><tr><td>TIL</td><td>tumor-infiltrated lymphocyte</td></tr><tr><td>Tn</td><td>naive T</td></tr><tr><td>Ttrm</td><td>tissue-resident memory T</td></tr></table>

## References

1. Hayday, A.C.; Vantourout, P. The Innate Biologies of Adaptive Antigen Receptors. *Annu. Rev. Immunol.* **2020**, *38*, 487–510. [CrossRef] [PubMed]

2. Dieli, F.; Poccia, F.; Lipp, M.; Sireci, G.; Caccamo, N.R.; DI Sano, C.; Salerno, A. Differentiation of Effector/Memory Vδ2 T Cells and Migratory Routes in Lymph Nodes or Inflammatory Sites. *J. Exp. Med.* **2003**, *198*, 391–397. [CrossRef] [PubMed]

3. Sturm, E.; Braakman, E.; Fisch, P.; Vreugdenhil, R.J.; Sondel, P.; Bolhuis, R.L. Human V Gamma 9-V Delta 2 T Cell Receptor-Gamma Delta Lymphocytes Show Specificity to Daudi Burkitt's Lymphoma Cells. *J. Immunol.* **1990**, *145*, 3202–3208. [PubMed]

4. Sciammas, R.; Johnson, R.M.; I Sperling, A.; Brady, W.; Linsley, P.S.; Spear, P.G.; Fitch, F.W.; Bluestone, J. Unique Antigen Recognition by a Herpesvirus-Specific Tcr-Gamma Delta Cell. *J. Immunol.* **1994**, *152*, 5392–5397.

5. Wallace, M.; Scharko, A.M.; Pauza, C.D.; Fisch, P.; Imaoka, K.; Kawabata, S.; Fujihashi, K.; Kiyono, H.; Tanaka, Y.; Bloom, B.R.; et al. Functional gamma delta T-lymphocyte defect associated with human immunodeficiency virus infections. *Molecular Med. Eng.* **1997**, *3*, 60–71. [CrossRef]

6. Poccia, F.; Gougeon, M.-; Agrati, C.; Montesano, C.; Martini, F.; Pauza, C.; Fisch, P.; Wallace, M.; Malkovsky, M. Innate T-Cell Immunity in HIV Infections: The Role of Vg9Vd2 T Lymphocytes. *Curr. Mol. Med.* **2002**, *2*, 769–781. [CrossRef]

7. Hong, K.; Greer, C.; Ketter, N.; Van Nest, G.; Paliard, X. Isolation and characterization of human papillomavirus type 6-specific T cells infiltrating genital warts. *J. Virol.* **1997**, *71*, 6427–6432. [CrossRef]

8. Ozaki, S.; Ogasahara, K.; Kosaka, M.; Inoshita, T.; Wakatsuki, S.; Uehara, H.; Matsumoto, T. Hepatosplenic gamma delta T-cell lymphoma associated with hepatitis B virus infection. *J. Med Investig.* **1998**, *44*, 215–217.

9. Qin, G.; Mao, H.; Zheng, J.; Sia, S.F.; Liu, Y.; Chan, P.L.; Lam, K.; Peiris, J.S.M.; Lau, Y.; Tu, W. Phosphoantigen-Expanded Human γδ T Cells Display Potent Cytotoxicity against Monocyte-Derived Macrophages Infected with Human and Avian Influenza Viruses. *J. Infect. Dis.* **2009**, *200*, 858–865. [CrossRef]

10. Sant, S.; Jenkins, M.R.; Dash, P.; Watson, K.; Wang, Z.; Pizzolla, A.; Koutsakos, M.; Nguyen, T.H.; Lappas, M.; Crowe, J.; et al. Human γδ T-cell receptor repertoire is shaped by influenza viruses, age and tissue compartmentalisation. *Clin. Transl. Immunol.* **2019**, *8*, e1079. [CrossRef]

11. Tsai, C.-Y.; Liong, K.H.; Gunalan, M.G.; Liang, L.D.S.; Lim, D.S.L.; Fisher, D.A.; Macary, P.A.; Leo, Y.S.; Wong, S.-C.; Puan, K.J.; et al. Type I IFNs and IL-18 Regulate the Antiviral Response of Primary Human γδ T Cells against Dendritic Cells Infected with Dengue Virus. *J. Immunol.* **2015**, *194*, 3890–3900. [CrossRef]

12. Poccia, F.; Agrati, C.; Castilletti, C.; Bordi, L.; Gioia, C.; Horejsh, D.; Ippolito, G.; Chan, P.; Hui, D.; Sung, J.J.Y.; et al. Anti-Severe Acute Respiratory Syndrome Coronavirus Immune Responses: The Role Played by Vγ9Vδ2 T Cells. *J. Infect. Dis.* **2006**, *193*, 1244–1249. [CrossRef]

13. Déchanet, J.; Merville, P.; Lim, A.; Retière, C.; Pitard, V.; Lafarge, X.; Michelson, S.; Méric, C.; Hallet, M.-M.; Kourilsky, P.; et al. Implication of γδ T cells in the human immune response to cytomegalovirus. *J. Clin. Investig.* **1999**, *103*, 1437–1449. [CrossRef]

14. Barragué, H.; Fontaine, J.; Abravanel, F.; Mauré, E.; Péron, J.-M.; Alric, L.; Dubois, M.; Izopet, J.; Champagne, E. Mobilization of γδ T Cells and IL-10 Production at the Acute Phase of Hepatitis E Virus Infection in Cytomegalovirus Carriers. *J. Immunol.* **2021**, *206*, 1027–1038. [CrossRef]

15. Barcy, S.; De Rosa, S.C.; Vieira, J.; Diem, K.; Ikoma, M.; Casper, C.; Corey, L. γδ+T Cells Involvement in Viral Immune Control of Chronic Human Herpesvirus 8 Infection. *J. Immunol.* **2008**, *180*, 3417–3425. [CrossRef]
16. Rijkers, G.; Vervenne, T.; Van Der Pol, P. More bricks in the wall against SARS-CoV-2 infection: Involvement of γ9δ2 T cells. *Cell. Mol. Immunol.* **2020**, *17*, 771–772. [CrossRef]

17. Carissimo, G.; Xu, W.; Kwok, I.; Abdad, M.Y.; Chan, Y.-H.; Fong, S.-W.; Puan, K.J.; Lee, C.Y.-P.; Yeo, N.K.-W.; Amrun, S.N.; et al. Whole blood immunophenotyping uncovers immature neutrophil-to-VD2 T-cell ratio as an early marker for severe COVID-19. *Nat. Commun.* **2020**, *11*, 1–12. [CrossRef]

18. Jouan, Y.; Guillon, A.; Gonzalez, L.; Perez, Y.; Boisseau, C.; Ehrmann, S.; Ferreira, M.; Daix, T.; Jeannet, R.; François, B.; et al. Phenotypical and functional alteration of unconventional T cells in severe COVID-19 patients. *J. Exp. Med.* **2020**, *217*. [CrossRef]

19. Li, C.-X.; Chen, J.; Lv, S.-K.; Li, J.-H.; Li, L.-L.; Hu, X. Whole-Transcriptome RNA Sequencing Reveals Significant Differentially Expressed mRNAs, miRNAs, and lncRNAs and Related Regulating Biological Pathways in the Peripheral Blood of COVID-19 Patients. *Mediat. Inflamm.* **2021**, *2021*, 1–22. [CrossRef]

20. Kalicińska, E.; Szymczak, D.; Zińczuk, A.; Adamik, B.; Smiechowicz, J.; Skalec, T.; Nowicka-Suszko, D.; Biernat, M.; Bogucka-Fedorczuk, A.; Rybka, J.; et al. Immunosuppression as a Hallmark of Critical COVID-19: Prospective Study. *Cells* **2021**, *10*, 1293. [CrossRef]

21. Shah, A.A.; Casciola-Rosen, L.; Rosen, A. Review: Cancer-Induced Autoimmunity in the Rheumatic Diseases. *Arthritis Rheumatol.* **2015**, *67*, 317–326. [CrossRef]

22. Meraviglia, S.; Presti, E.L.; Tosolini, M.; La Mendola, C.; Orlando, V.; Todaro, M.; Catalano, V.; Stassi, G.; Cicero, G.; Vieni, S.; et al. Distinctive features of tumor-infiltrating γδ T lymphocytes in human colorectal cancer. *OncoImmunology* **2017**, *6*, e1347742. [CrossRef]

23. Gentles, A.J.; Newman, A.; Liu, C.L.; Bratman, S.; Feng, W.; Kim, D.; Nair, V.S.; Xu, Y.; Khuong, A.; Hoang, C.D.; et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. *Nat. Med.* **2015**, *21*, 938–945. [CrossRef]

24. Wu, Y.; Kyle-Cezar, F.; Woolf, R.T.; Naceur-Lombardelli, C.; Owen, J.; Biswas, D.; Lorenc, A.; Vantourout, P.; Gazinska, P.; Grigoriadis, A.; et al. An innate-like Vδ1 + γδ T cell compartment in the human breast is associated with remission in triple-negative breast cancer. *Sci. Transl. Med.* **2019**, *11*, eaax9364. [CrossRef]

25. Mikulak, J.; Oriolo, F.; Bruni, E.; Roberto, A.; Colombo, F.S.; Villa, A.; Bosticardo, M.; Bortolomai, I.; Presti, E.L.; Meraviglia, S.; et al. NKp46-expressing human gut-resident intraepithelial Vδ1 T cell subpopulation exhibits high antitumor activity against colorectal cancer. *JCI Insight* **2019**, *4*, 125884. [CrossRef]

26. Foord, E.; Arruda, L.C.M.; Gaballa, A.; Klynning, C.; Uhlin, M. Characterization of ascites- and tumor-infiltrating γδ T cells reveals distinct repertoires and a beneficial role in ovarian cancer. *Sci. Transl. Med.* **2021**, *13*, 0192. [CrossRef]

27. Pizzolato, G.; Kaminski, H.; Tosolini, M.; Franchini, D.-M.; Pont, F.; Martins, F.; Valle, C.; Labourdette, D.; Cadot, S.; Quillet-Mary, A.; et al. Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVδ1 and TCRVδ2 γδ T lymphocytes. *Proc. Natl. Acad. Sci. USA* **2019**, *116*, 11906–11915. [CrossRef]

28. Cerapio, J.-P.; Perrier, M.; Balança, C.-C.; Gravelle, P.; Pont, F.; Devaud, C.; Franchini, D.-M.; Féliu, V.; Tosolini, M.; Valle, C.; et al. Phased differentiation of γδ T and T CD8 tumor-infiltrating lymphocytes revealed by single-cell transcriptomics of human cancers. *OncoImmunology* **2021**, *10*, 1939518. [CrossRef]

29. Aoki, T.; Chong, L.C.; Takata, K.; Milne, K.; Hav, M.; Colombo, A.; Chavez, E.A.; Nissen, M.; Wang, X.; Miyata-Takata, T.; et al. Single-Cell Transcriptome Analysis Reveals Disease-Defining T-cell Subsets in the Tumor Microenvironment of Classic Hodgkin Lymphoma. *Cancer Discov.* **2019**, *10*, 406–421. [CrossRef]

30. Cillo, A.R.; Kürten, C.H.; Tabib, T.; Qi, Z.; Onkar, S.; Wang, T.; Liu, A.; Duvvuri, U.; Kim, S.; Soose, R.J.; et al. Immune Landscape of Viral- and Carcinogen-Driven Head and Neck Cancer. *Immuni* **2020**, *52*, 183–199.e21. [CrossRef]

31. Stuart, T.; Butler, A.; Hoffman, P.; Hafemeister, C.; Papalexi, E.; Mauck, W.M., III; Hao, Y.; Stoeckius, M.; Smibert, P.; Satija, R. Comprehensive Integration of Single-Cell Data. *Cell* **2019**, *177*, 1888–1902. [CrossRef] [PubMed]

32. McInnes, L.; Healy, J. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. Available online: https://arxiv.org/abs/1802.03426 (accessed on 12 October 2021)

33. Pont, F.; Tosolini, M.; Fournié, J.J. Single-Cell Signature Explorer for comprehensive visualization of single cell signatures across scRNA-seq datasets. *Nucleic Acids Res.* **2019**, *47*, e133. [CrossRef] [PubMed]

34. Pont, F.; Tosolini, M.; Gao, Q.; Perrier, M.; Madrid-Mencía, M.; Huang, T.S.; Neuvial, P.; Ayyoub, M.; Nazor, K.; Fournié, J.-J. Single-Cell Virtual Cytometer allows user-friendly and versatile analysis and visualization of multimodal single cell RNAseq datasets. *NAR Genom. Bioinform.* **2020**, *2*, qaa025. [CrossRef] [PubMed]

35. Chihara, N.; Madi, A.; Kondo, T.; Zhang, H.; Acharya, N.; Singer, M.; Nyman, J.; Marjanovic, N.D.; Kowalczyk, M.S.; Wang, C.; et al. Induction and transcriptional regulation of the co-inhibitory gene module in T cells. *Nat. Cell Biol.* **2018**, *558*, 454–459. [CrossRef]

36. Alfei, F.; Kanev, K.; Hofmann, M.; Wu, M.; Ghoneim, H.E.; Roelli, P.; Utzschneider, D.T.; Von Hoesslin, M.; Cullen, J.G.; Fan, Y.; et al. TOX reinforces the phenotype and longevity of exhausted T cells in chronic viral infection. *Nature* **2019**, *571*, 265–269. [CrossRef]

37. Khan, O.; Giles, J.R.; McDonald, S.; Manne, S.; Ngiow, S.F.; Patel, K.P.; Werner, M.T.; Huang, A.C.; Alexander, K.A.; Wu, J.E.; et al. TOX transcriptionally and epigenetically programs CD8+ T cell exhaustion. *Nature* **2019**, *571*, 211–218. [CrossRef]

38. Tosolini, M.; Algans, C.; Pont, F.; Ycart, B.; Fournié, J.-J. Large-scale microarray profiling reveals four stages of immune escape in non-Hodgkin lymphomas. *OncoImmunology* **2016**, *5*, e1188246. [CrossRef]
39. Balança, C.-C.; Scarlata, C.-M.; Michelas, M.; Devaud, C.; Sarradin, V.; Franchet, C.; Gomez, C.M.; Gomez-Roca, C.; Tosolini, M.; Heaugwane, D.; et al. Dual Relief of T-lymphocyte Proliferation and Effector Function Underlies Response to PD-1 Blockade in Epithelial Malignancies. *Cancer Immunol. Res.* **2020**, *8*, 869–882. [CrossRef]

40. Kumar, B.V.; Ma, W.; Miron, M.; Granot, T.; Guyer, R.S.; Carpenter, D.J.; Senda, T.; Sun, X.; Ho, S.-H.; Lerner, H.; et al. Human Tissue-Resident Memory T Cells Are Defined by Core Transcriptional and Functional Signatures in Lymphoid and Mucosal Sites. *Cell Rep.* **2017**, *20*, 2921–2934. [CrossRef]

41. Wu, T.D.; Madireddi, S.; de Almeida, P.E.; Banchereau, R.; Chen, Y.-J.J.; Chitre, A.S.; Chiang, E.Y.; Iftikhar, H.; O'Gorman, W.E.; Au-Yeung, A.; et al. Peripheral T cell expansion predicts tumour infiltration and clinical response. *Nat. Cell Biol.* **2020**, *579*, 274–278. [CrossRef]

42. Alpert, A.; Moore, L.S.; Dubovik, T.; Shen-Orr, S.S. Alignment of single-cell trajectories to compare cellular expression dynamics. *Nat. Methods* **2018**, *15*, 267–270. [CrossRef]

43. Trapnell, C.; Cacchiarelli, D.; Grimsby, J.; Pokharel, P.; Li, S.; Morse, M.; Lennon, N.J.; Livak, K.J.; Mikkelsen, T.S.; Rinn, J.L. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. *Nat. Biotechnol.* **2014**, *32*, 381–386. [CrossRef]

44. Berge, K.V.D.; de Bézieux, H.R.; Street, K.; Saelens, W.; Cannoodt, R.; Saeys, Y.; Dudoit, S.; Clement, L. Trajectory-based differential expression analysis for single-cell sequencing data. *Nat. Commun.* **2020**, *11*, 1–13. [CrossRef]

45. Lo Presti, E.; Dieli, F. and Meraviglia, S. Lymphopenia in COVID-19: Gammadelta T Cells-Based Therapeutic Opportunities. *Vaccines (Basel)* **2021**, *9*, 562. [CrossRef]

46. Caron, J.; Ridgley, L.A.; Bodman-Smith, M. How to Train Your Dragon: Harnessing Gamma Delta T Cells Antiviral Functions and Trained Immunity in a Pandemic Era. *Front. Immunol.* **2021**, *12*, 666983. [CrossRef]

47. Yazdanifar, M.; Mashkour, N.; Bertaina, A. Making a case for using γδ T cells against SARS-CoV-2. *Crit. Rev. Microbiol.* **2020**, *46*, 689–702. [CrossRef]

48. Orumaa, K.; Dunne, M.R. The role of unconventional T cells in COVID-19. *Ir. J. Med Sci.* **2021**, 1–10. [CrossRef]

49. Lei, L.; Qian, H.; Yang, X.; Zhang, X.; Zhang, D.; Dai, T.; Guo, R.; Shi, L.; Cheng, Y.; Zhang, B.; et al. The phenotypic changes of γδ T cells in COVID-19 patients. *J. Cell. Mol. Med.* **2020**, *24*, 11603–11606. [CrossRef]

50. Moratto, D.; Chiarini, M.; Giustini, V.; Serana, F.; Magro, P.; Roccaro, A.M.; Imberti, L.; Castelli, F.; Notarangelo, L.D.; Quiros-Roldan, E. Flow Cytometry Identifies Risk Factors and Dynamic Changes in Patients with COVID-19. *J. Clin. Immunol.* **2020**, *40*, 970–973. [CrossRef]

51. Carter, M.J.; Fish, M.; Jennings, A.; Doores, K.J.; Wellman, P.; Seow, J.; Acors, S.; Graham, C.; Timms, E.; Kenny, J.; et al. Peripheral immunophenotypes in children with multisystem inflammatory syndrome associated with SARS-CoV-2 infection. *Nat. Med.* **2020**, *26*, 1701–1707. [CrossRef]

52. Guo, C.; Li, B.; Ma, H.; Wang, X.; Cai, P.; Yu, Q.; Zhu, L.; Jin, L.; Jiang, C.; Fang, J.; et al. Single-cell analysis of two severe COVID-19 patients reveals a monocyte-associated and tocilizumab-responding cytokine storm. *Nat. Commun.* **2020**, *11*, 1–11. [CrossRef]

53. Xu, X.; Han, M.; Li, T.; Sun, W.; Wang, D.; Fu, B.; Zhou, Y.; Zheng, X.; Yang, Y.; Li, X.; et al. Effective treatment of severe COVID-19 patients with tocilizumab. *Proc. Natl. Acad. Sci. USA* **2020**, *117*, 10970–10975. [CrossRef]

54. Presti, E.L.; Corsale, A.M.; Dieli, F.; Meraviglia, S. γδ cell-based immunotherapy for cancer. *Expert Opin. Biol. Ther.* **2019**, *19*, 887–895. [CrossRef]

55. Simoni, Y.; Becht, E.; Fehlings, M.G.; Loh, C.Y.; Koo, S.-L.; Teng, K.W.W.; Yeong, J.P.S.; Nahar, R.; Zhang, T.; Kared, H.; et al. Bystander CD8+ T cells are abundant and phenotypically distinct in human tumour infiltrates. *Nature* **2018**, *557*, 575–579. [CrossRef]

56. Oliveira, G.; Stromhaug, K.; Klaeger, S.; Kula, T.; Frederick, D.T.; Le, P.M.; Forman, J.; Huang, T.; Li, S.; Zhang, W.; et al. Phenotype, specificity and avidity of antitumour CD8+ T cells in melanoma. *Nat. Cell Biol.* **2021**, *596*, 119–125. [CrossRef]