Article Text

Original research
Intrahepatic CD69+Vδ1 T cells re-circulate in the blood of patients with metastatic colorectal cancer and limit tumor progression
  1. Elena Bruni1,2,
  2. Matteo Maria Cimino3,
  3. Matteo Donadon3,4,
  4. Roberta Carriero5,
  5. Sara Terzoli1,6,
  6. Rocco Piazza7,
  7. Sarina Ravens8,
  8. Immo Prinz8,9,
  9. Valentina Cazzetta1,2,
  10. Paolo Marzano1,2,
  11. Paolo Kunderfranco5,
  12. Clelia Peano10,
  13. Cristiana Soldani11,
  14. Barbara Franceschini11,
  15. Federico Simone Colombo12,
  16. Cecilia Garlanda6,13,
  17. Alberto Mantovani6,13,14,
  18. Guido Torzilli3,6,
  19. Joanna Mikulak1,2 and
  20. Domenico Mavilio1,2
  1. 1Laboratory of Clinical and Experimental Immunology, IRCCS Humanitas Research Hospital, Rozzano, Milan, Italy
  2. 2Department of Medical Biotechnologies and Translational Medicine, University of Milan, Milan, Italy
  3. 3Department of Hepatobiliary and General Surgery, IRCCS Humanitas Research Hospital, Rozzano, Milan, Italy
  4. 4Department of Health Science, Università del Piemonte Orientale, Novara, Italy
  5. 5Bioinformatics Unit, IRCCS Humanitas Research Hospital, Rozzano, Milan, Italy
  6. 6Department of Biomedical Sciences, Humanitas University, Pieve Emanuele, Milan, Italy
  7. 7Department of Medicine and Surgery, University of Milan-Bicocca, Monza, Italy
  8. 8Institute of Immunology, Hannover Medical School (MHH), Hannover, Germany
  9. 9Institute of Systems Immunology, Hamburg Center for Translational Immunology (HCTI), University Medical Center Hamburg-Eppendorf, Hamburg, Germany
  10. 10Institute of Biomedical Technologie, CNR Milan, Human Technopole, Milan, Italy
  11. 11Hepatobiliary Immunopathology Laboratory, IRCCS Humanitas Research Hospital, Rozzano, Milan, Italy
  12. 12Humanitas Flow Cytometry Core, IRCCS Humanitas Research Hospital, Rozzano, Milan, Italy
  13. 13IRCCS Humanitas Research Hospital, Rozzano, Milan, Italy
  14. 14The William Harvey Research Institute, Queen Mary University of London, London, UK
  1. Correspondence to Dr Domenico Mavilio; domenico.mavilio{at}humanitas.it
  • JM and DM are joint last authors.

Abstract

Background More than 50% of all patients with colorectal cancer (CRC) develop liver metastases (CLM), a clinical condition characterized by poor prognosis and lack of reliable prognostic markers. Vδ1 cells are a subset of tissue-resident gamma delta (γδ) T lymphocytes endowed with a broad array of antitumor functions and showing a natural high tropism for the liver. However, little is known about their impact in the clinical outcomes of CLM.

Methods We isolated human γδ T cells from peripheral blood (PB) and peritumoral (PT) tissue of 93 patients undergone surgical procedures to remove CLM. The phenotype of freshly purified γδ T cells was assessed by multiparametric flow cytometry, the transcriptional profiles by single cell RNA-sequencing, the functional annotations by Gene Ontology enrichment analyses and the clonotype by γδ T cell receptor (TCR)-sequencing.

Results The microenvironment of CLM is characterized by a heterogeneous immune infiltrate comprising different subsets of γδ tumor-infiltrating lymphocytes (TILs) able to egress the liver and re-circulate in PB. Vδ1 T cells represent the largest population of γδ TILs within the PT compartment of CLM that is greatly enriched in Vδ1 T effector (TEF) cells expressing constitutive high levels of CD69. These Vδ1 CD69+ TILs express a distinct phenotype and transcriptional signature, show high antitumor potential and correlate with better patient clinical outcomes in terms of lower numbers of liver metastatic lesions and longer overall survival (OS). Moreover, intrahepatic CD69+ Vδ1 TILs can egress CLM tissue to re-circulate in PB, where they retain a phenotype, transcriptional signature and TCR clonal repertoires resembling their liver origin. Importantly, even the increased frequencies of the CD69+ terminally differentiated (TEMRA) Vδ1 cells in PB of patients with CLM significantly correlate with longer OS. The positive prognostic score of high frequencies of CD69+ TEMRA Vδ1 cells in PB is independent from the neoadjuvant chemotherapy and immunotherapy regimens administered to patients with CLM prior surgery.

Conclusions The enrichment of tissue-resident CD69+ Vδ1 TEMRA cells re-circulating at high frequencies in PB of patients with CLM limits tumor progression and represents a new important clinical tool to either predict the natural history of CLM or develop alternative therapeutic protocols of cellular therapies.

  • immunity, cellular
  • lymphocytes, tumor-infiltrating
  • T-lymphocytes
  • liver neoplasms
  • immunologic surveillance

Data availability statement

Data are available on reasonable request.

http://creativecommons.org/licenses/by-nc/4.0/

This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See http://creativecommons.org/licenses/by-nc/4.0/.

Statistics from Altmetric.com

Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.

WHAT IS ALREADY KNOWN ON THIS TOPIC

  • The metastatic progression of colorectal cancer (CRC) to livers occurs in approximately 50% of patients and still represents an unmet clinical need lacking a cure.

  • The microenvironment of human liver is characterized by a heterogeneous distribution of distinct subsets of gamma delta (γδ) T lymphocytes with antitumor effector potential.

WHAT THIS STUDY ADDS

  • Intrahepatic CD69+ Vδ1 T cells in liver metastatic CRC (CLM) represent the predominant TIL subset with high antitumor effector functions that is also able to egress tumor and re-circulate in peripheral blood (PB).

  • Higher frequencies of both TI and PB CD69+ Vδ1 T cells are associated with a slower liver metastatic progression and a longer overall survival of patients with CLM.

HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY

  • The enrichment of tissue-resident CD69+ Vδ1 T cells re-circulating at high frequencies in PB represents a new important clinical tool to help predicting either the natural history of CLM in ‘liquid biopsies’ or develop alternative therapeutic protocols of cellular therapies.

Introduction

The metastatic progression of colorectal cancer (CRC) to livers occurs in approximately 50% of patients and still represents an unmet clinical need lacking a cure. The surgical resection of colon liver metastases (CLM) is the only available treatment increasing patients’ overall survival (OS), while definitive and shared clinical criteria that stratifies patients’ eligibility to surgery and/or neoadjuvant conventional/biological chemotherapies (na-CHT) are still missing.1

Targeting different aspects of antitumor immune responses in primary or metastatic tumors (MTs) is key for the development of new clinical tools to predict cancer progression and develop alternative therapeutic strategies. Indeed, since the diagnostic potential of resected/biopsied tumors is hampered by multiple methodological and ethical constraints, there is a growing interest in finding immune prognostic variables within the so-called ‘liquid biopsies’ from peripheral blood (PB). In this scenario, tumor-infiltrating (TI) lymphocytes (TILs) have been also found in the PB of patients affected by different cancers.2 3 The development of new single-cell technologies with an unprecedented degree of immunological resolution can now identify in ‘liquid biopsies’ T cells showing tumor specificities similar to those of their homologs in neoplastic lesions.4 5

Gamma delta (γδ) T cells are a unique subset of T lymphocytes endowed with innate-like features due to their major histocompatibility complex (MHC)-unrestricted ability to activate and respond to stimuli.6 They represent the first cluster of cells appearing during T cell differentiation and play a crucial role in cancer immune-surveillance given their potent cytotoxicity and production of several antitumor cytokines.7 Indeed, the presence of γδ T cells at high frequencies emerged as one of the most favorable and widest prognostic factors in cancer progressions (including CRC).8 Transcriptomic data and molecular analysis of human γδ TCR repertoire revealed the existence of two major subsets of γδ T cells, namely Vδ1 and Vδ2. The first one is highly enriched at tissue sites and owns features of local resident cells endowed with remarkable capacities in the clearance of malignant cells from different hematopoietic and solid cancers.9–13 Vδ2 T cells preferentially circulate in the bloodstream thus being targeted for developing protocols of in vitro expansion for adoptive cell transfer therapies.14 15 Moreover, a third population expressing Vδ3-TCR chain had been also described. These latter Vδ3 T cells exhibit a high liver tropism and are present at very low frequency in PB under homeostatic conditions, while expanding in bloodstream during viral infections and inflammatory diseases.16 17 However, a comparative cellular and molecular characterization of TI and PB γδ T cells in human cancers is still lacking.

Vδ1 T cells represent the main subset of γδ T lymphocytes in human liver. They include both tissue-resident and re-circulating cells playing a key role in hepatic immune homeostasis.18 19 Nothing is known about the clinical impact of hepatic Vδ1 T cells in limiting the progression of CLM. In this regard, we recently identified a specific subset of human intraepithelial gut NKp46+ Vδ1 T cells able to limit the metastatic diffusion of primary CRC to liver.20 Similar results have been also confirmed in mice.21 Moreover, increased frequencies of PB and intrahepatic γδ T cells have been reported during chronic hepatitis and liver inflammation.22 23

The present study provides an in-depth and integrated phenotypical and transcriptional characterization of TI and PB γδ T cell subset in patients with CLM. In particular, we identified a subset of CD69+ Vδ1 TILs whose higher frequencies in tumor specimens are associated with lower numbers of metastatic lesions and longer OS of patients with CLM. Importantly, these latter immune cell subsets can egress the tumor and re-circulate in PB as CD69+ terminally differentiated (TEMRA) Vδ1 cells. Taken together, these findings could remarkably improve our ability to predict the prognosis of patients with CLM and pave the ground for developing new strategies of adoptive cell transfer immunotherapies.

Materials and methods

Study design

The study characterized γδ T cells in a real-life cohort of total 93 patients with CRC underwent hepatic resection of rectal metastasis mostly represented by older patients with mean age of 64 years, mainly synchronous clinical profile who received na-CHT treatment (table 1). We assessed phenotype, function, γδ TCR clonotype and transcriptional profiles by single cell RNA-sequencing (scRNA-seq) of γδ T cells isolated from PB, tumor-free liver parenchyma surrounding CLM and MT.

Table 1

Cohort characterization of patients with CLM

Cell preparation

Peripheral blood mononuclear cells obtained from the buffy coats of healthy donors (HDs) or from blood of patients with CLM were isolated by Lympholyte Cell Separation density gradient solution (Cederlane Laboratories, Burlington, Canada) according to the manufacturer’s instructions. Liver specimens were dissociated by enzymatic digestion with 2 mg/mL of collagenase D (Roche Diagnostic; Indianapolis, Indiana, USA) in gentleMACS Dissociator (Miltenyi; Bergisch Gladbach, Germany) as previously described.20 Lymphocytes were separated by 70%/30% discontinuous Percoll gradient (GE Healthcare; Chicago, Illinois, USA) and frozen in fetal bovine serum (Lonza, Verviers, Belgium) with 10% of dimethylsulfoxide (Lonza) before being stored in liquid nitrogen for further analysis.

Additional materials and methods are included in online supplemental materials.

Supplemental material

Results

Vδ1 T cells in CLM predict longer patients’ overall survival

We first characterized by immunohistochemistry (IHC) the frequency and distribution of TILs within CLM. Among CD3+ T cells, γδ TILs are significantly higher in numbers within peritumoral tissue (PT) compared with the core of highly fibrotic MT (figure 1A,B and online supplemental figure 1A). Additionally, multiparametric flow cytometry data showed that Vδ1 T cells represent the main subset of γδ TILs in CLM, as their frequencies are significantly higher compared with Vδ2 T cells (figure 1C and online supplemental figure 1B). To assess the phenotype of γδ TILs, we evaluated the expression of CD27 and CD45RA that are recognized markers of αβ T cell differentiation. Moreover, the surface levels of both these markers are also indicative of the homeostatic status and effector functions of human γδ T cells including proliferative capacity, resistance to drug treatment, cytokine secretion and cytotoxic response.24–28 Vδ1 TILs with the CD27 phenotype characterize two effector (TEF) cell subpopulations that include TEMRA (CD45RA+CD27) and T effector memory (TEM; CD45RACD27) cells. These cells represent the larger fractions of Vδ1 in CLM, while phenotypes of naïve T (TNAIVE) (CD45RA+CD27+) and T central memory (TCM; CD45RACD27+) Vδ1 T cells are present at very low frequencies. This is not the case of Vδ2 TILs, whose distribution in CLM is characterized by the predominant expansion of cells with the phenotype of TEM together with lower frequencies of TCM and TEMRA (figure 1D). As expected,18 29 both Vδ1 and Vδ2 TILs express high levels of the tissue-residency and activating marker CD69. In accordance with more mature phenotype, Vδ1 TILs display low levels of CD27 and CD28, both markers of immature T cells (figure 1E and online supplemental figure 1C). Compared with the Vδ2 counterpart, Vδ1 TILs present higher frequencies of CD57, a T cell marker of proliferative-senescence found to be increased in Vδ2 T cells of patients with CLM.24 Importantly, expression of CD57 in Vδ1 TILs inversely correlates with the expression of CD69 and several effector molecules such as NKG2D and CD56 (online supplemental figure 1C).

Supplemental material

Figure 1

Cellular characterization and clinical impact of intrahepatic gamma delta (γδ) tumor-infiltrating lymphocytes (TILs) in colon liver metastatic cancer (CLM). (A) Representative immunohistochemistry (IHC) images showing CD3+ (upper panels) and γδ T cell receptor (TCR)+ (lower panels) T cells within peritumoral tissue (PT) and metastatic tumor (MT) from one CLM liver section (out of 10) at 10× (left panel) and 20× (right panel) magnified views. Arrows indicate γδ TCR+ cells in lower panels. (B) Statistical graph from IHC data showing the mean frequency (%) of γδ TCR+ TILs among CD3+ T cells in PT and MT areas (n=10). (C) Statistical graph from multiparametric flow cytometry data showing the mean (±SEM) frequency (%) of matched Vδ1 and Vδ2 T cells among CD3+ T cells in PT (n=82). (D) Pie charts of multiparametric flow cytometry data showing the mean distribution (%) of naïve T (TNAIVE), T central memory (TCM), T effector memory (TEM) and terminally differentiated (TEMRA) Vδ1 (upper chart) and Vδ2 (lower chart) T cells within PT (n=62). (E) Heatmaps of multiparametric flow cytometry data showing the expression of several surface markers between matched Vδ1 (left panel) and Vδ2 (right panel) TILs from PT (n=62). Group of markers which statistically differ between Vδ1 and Vδ2 TILs are annotated under the heatmap as ‘≥****’ and indicates the range of significant p values; ‘ns’ indicates a group of markers not statistically significant. (F) Kaplan-Meier curve of postsurgical overall survival (OS; %) of patients with CLM based on the number of liver metastases (MT, cut-off ≥4; n=60). (G) Spearman’s rank correlation between the frequency (%) of liver PT CD69+ Vδ1 T cells and the number of surgically removed liver metastasis (n=60). (H) Kaplan-Meier curves showing postsurgical OS (%) of patients with CLM based on the median frequencies (%) of CD69+CD27 (cut-off ≤2%; left panel) and CD69+CD28 liver Vδ1 T cells in PT (cut-off ≤5%; right panel) (n=63) Statistically significant p values are represented with the following number of asterisks (*): *p ≤ 0.05; ***p ≤ 0.001; ****p ≤ 0.0001.

In addition, the surface levels of activating receptors NKG2D, CD56 and CD16 are similar between Vδ1 and Vδ2 TILs. On the other hand, Vδ2 TILs are characterized by the constitutive high expression of the inhibitory checkpoint NKG2A (figure 1E).25 30 31 We then assessed the prognostic values of Vδ1 TILs. In line with reported clinical data,32 we found that patients with CLM with less than four metastatic lesions had a significant longer OS compared with those individuals having four or more lesions (figure 1F). Importantly, higher frequencies of CD69+ Vδ1 TILs cells in CLM significantly correlate with lower numbers of metastatic lesions (figure 1G). We then combined flow cytometry data with the clinical outcome of patients with CLM stratified for high and low cut-offs of CD69+CD27 and CD69+CD28 Vδ1 effector T cells. Our results showed that high frequencies of both these CD69+ Vδ1 TIL subsets are associated with significantly longer OS, thus indicating their key role in limiting the metastatic progression of CRC to liver (figure 1H).

Transcriptional heterogeneity of γδ TILs in CLM

To dissect the nature and functional relevance of γδ TILs in the CLM microenvironment, we performed scRNA-seq analysis. To avoid sample selection bias, for scRNA-seq experiments we selected three representative patients (table 1) who underwent limited liver resection for synchronous CLM disease and were treated with standard combination na-CHT, prior surgery. PT and distal tumor-free (DT) γδ TILs were performed starting from total CD3+ cells as described in online supplemental materials and methods. The unsupervised clustering obtained by performing dimensionality reduction based on uniform manifold approximation and projection (UMAP) identified eight clusters (c0–c7) among total γδ T cells in PT and DT compartments (figure 2A). The analysis of TCR δ chain expression distinguishes three clusters of Vδ1 cells (c0, c3 and c4), while two cell clusters (c2 and c5) show a Vδ2-TCR repertoire. We also identified two additional clusters of cells (c1 and c6) expressing Vδ3 chain, thus confirming the presence of Vδ3 T cells in the human liver (figure 2B).18 Cells lacking expression of TRDV1-3 (c7) were excluded from further analysis. In line with our flow cytometry results (figure 1C,D), the relative enrichment analysis showed that Vδ1 TILs represent the dominant subset among all γδ T cells in both PT and DT compartments (figure 2C). In particular, the cluster distribution in DT is characterized by the enrichment of clusters c0 and c4 (Vδ1 TILs) and c2 (Vδ2 TILs). On the other side, PT was mainly enriched in clusters belonging to Vδ1 (c0 and c3) and Vδ3 (c1) TILs (figure 2D). Regarding Vγ chain recombination, Vδ1 TILs were mainly associated with Vγ3 and Vγ4, Vδ2 TILs with Vγ9 and Vδ3 TILs with Vγ5 chain (figure 2E).

Figure 2

Single cell RNA-sequencing of gamma delta (γδ) tumor-infiltrating lymphocytes (TILs) in colon liver metastatic cancer (CLM). (A) Uniform manifold approximation and projection (UMAP) projection of γδ T cells from peritumor (PT) (n=1333) and distal tumor-free (DT) (n=424) compartments of CLM from three patients underwent surgical liver resection. UMAP graph identifies eight specific clusters (left panel) and their distribution among DT (orange cells) and PT (light blue cells) compartments of CLM (right panel). (B) Heatmap showing the average of T cell receptor (TCR) δ chain expression within the eight γδ T cell clusters identified by UMAP analysis. The expression values are zero-centered and scaled for each gene. (C) Pie charts showing the relative enrichment of Vδ1, Vδ2 and Vδ3 T cells (%) among DT (left panel) and PT (right panel) compartments of CLMs. (D) Bar plot graph showing the distribution (%) of γδ T cell clusters among DT and PT compartments of CLM within the eight clusters identified from UMAP analysis. The proportion of γδ T cells for each cluster across the DT and PT compartments were calculated as ratio between number of cells in each cluster and total number of cells in DT and PT, respectively. (E) Violin plot graph showing the expression of TRGV genes among the seven γδ T cell clusters grouped according to their Vδ1, Vδ2 and Vδ3 T cell origin. (F) Violin plot graphs of selected genes among all γδ T cell clusters grouped according to their Vδ1, Vδ2 and Vδ3 T cell origin. Each graph shows genes associated with cell differentiation status, transcription factors, tissue-affinity, cytotoxicity, inhibitory and activating molecules.

We then analyzed the different transcriptional profiles of the seven γδ TIL clusters identified in CLM (figure 2F, figure 3). The tissue tropism of all these clusters was confirmed by expression of tissue and liver residency markers CD69 and CXCR6, respectively. However, lower transcript levels of CXCR6 and CD69 were observed in c4 (ie, Vδ1 T cells), thus likely indicating their re-circulation from PB. Consistently with our flow cytometry data, Vδ1 TILs (c0, c3, c4) and Vδ3 TILs (c1, c6) express a CD27low/CD28low phenotype. On the other hand, Vδ2 TILs (c2, c5) express higher transcript levels of CD27, CD28 and IL7R, thus showing an immature phenotype. Of note, CD27low Vδ2 TILs in c5 are characterized by the specific signature of previously reported type 3 Vδ2 T cells (CCR6, RORC, IL23R, DPP4).33

Figure 3

Differentially expressed genes among gamma delta (γδ) tumor-infiltrating lymphocyte (TIL) clusters in colon liver metastatic cancer (CLM). Heatmap showing the values of 238 differentially expressed genes (DEGs) (adjusted p<0.05) coming from pairwise comparison between cluster 0 (c0) against each identified γδ T cluster. Expression values are zero-centered and scaled for each gene.

Besides some transcriptional factors (TFs) such as ZNF331 and JUNB that are widely expressed in all γδ TIL clusters, we detected specific TF signatures among different γδ T cell types. Indeed, Vδ2 (c2 and c5) and Vδ3 (c1 and c6) T cells are characterized by high transcript levels of ZBTB16 (PLZF) and EOMES, respectively. Conversely, the expression of IKZF2 together with the concomitant lack of EOMES defines Vδ1 T cells (c0, c3, c4). The expression of genes associated with inhibitory signaling and tumor-promoting tolerance also revealed specific signatures among γδ TILs in CLM. Indeed, high transcriptional levels of inhibitory Killer-cell immunoglobulin-like receptors (iKIRs) are present on Vδ1 (c0, c3, c4) and Vδ3 (c1, c6) TILs, while the expressions of NKG2A (KLRC1) and KLRG1 distinguish Vδ2 TILs. Regarding immune-checkpoints expression, we found low/negative expression of PDCD1 (programmed cell death protein 1) and CTLA4, across all γδ TIL clusters with the only exception of TIGIT that resulted highly expressed mainly on Vδ3 T cells (c1, c6).

In general, all seven identified γδ TIL clusters express high levels of several cytotoxic molecules (GZMA, GZMB, GZMH, GZMK, GZMM, PRF1 and CTSW) and the activating receptor NKG7. Specifically, Vδ1 TILs (c0, c3, c4) can be distinguished for their high transcript levels of KLRC3 (NKG2E) and NCR3 (NKp30), while the expression of KLRB1 (CD161) characterizes Vδ2 (c2, c5) and Vδ3 (c1, c6) TILs in CLM. Taken together, these transcriptional profiles demonstrate the overall high cytotoxic potentials of intrahepatic TILs in CLM and make it possible to define specific pathways and functional commitments among the three subsets of γδ T cells infiltrating human liver during metastatic progression of CRC.

Functional commitments of γδ TILs in CLM

The functional annotations of γδ TILs were performed by Gene Ontology (GO) enrichment analysis on sc-RNA seq data from CLM specimens. We first observed that all γδ TILs share functional pathways highly associated with T cell activation, immune effector responses, leukocytes proliferation and regulation of cell-to-cell adhesion (figure 4A). The different ranges of functional annotations across specific γδ TIL subsets are due to the diverse degrees of gene transcripts and their response/adaptation to tumor microenvironment. In this context, Vδ2 TILs can be distinguished for the enrichment in pathways associated with T cell differentiation and B cell activation. As expected,33 type 3 Vδ2 TILs in c5 are specifically linked to the functional pathway of T helper 17 (Th-17)-IL-17-producing cells. However, we did not observe any IL-17A production. Instead, we found higher transcriptional level of tumor necrosis factor (TNF) among type 3 Vδ2 TILs (figures 3 and 4C).

Figure 4

Functional annotations of gamma delta (γδ) tumor-infiltrating lymphocytes (TILs) in colon liver metastatic cancer (CLM). (A) Statistical dot plot graph showing the biological processes (BO) obtained by Gene Ontology (GO) enrichment analysis calculated for total Vδ1, Vδ2, Vδ3 TILs and for the clusters c3 of Vδ1, c5 of Vδ2 and c6 of Vδ3 TIL subsets. GO enrichment analyses are performed by using differentially expressed genes (DEGs) with adjusted p≤0.05 and log-foldchange >0. Only enriched GO terms with adjusted p≤0.05 and more than five genes are reported. (B) Statistical bar plot graph showing interferon-gamma (IFN-γ), tumor necrosis factor (TNF) and X-C motif chemokine ligand 2 (XCL2) cytokine modular score calculated for each γδ T cell clusters. (C) Uniform manifold approximation and projection (UMAP) density plot graphs showing IFN-γ, TNF and XCL2 gene density distribution among γδ TIL cell clusters. (D) Scatter dot plot distribution with Pearson’s correlation between TRDV1 and Vδ1-Cy-SI expression in The Cancer Genome Atlas (TCGA) and The Genotype-Tissue Expression (GTEx) cohorts of patients with CRC calculated by Gene Expression Profiling Interactive Analysis 2 (GEPIA2). (E) Kaplan-Meier curve showing postoperative disease-free survival (%) ranking in TCGA cohorts of low-frequency microsatellite instability (MSI-L) and microsatellite stability (MSS) patients with CLM stratified by the low (blue curve) or high (red curve) Vδ1-Cy-SI expression (HR 0.5; median high cut-off 50%–50%; n=110). (F) Kaplan-Meier curve (left panel) showing the OS (%) in TCGA cohorts of MSI-L and MSS patients with CLM stratified by the low (blue curve) or high (red curve) Vδ1-Cy-SI ranking (HR 0.5; median 50%–50% high cut-off; n=110). Forest plots (right panel) showing HR with 95% CI, p value (p) and number (n) of patients obtained for different high Vδ1-Cy-SI cut-off values in TCGA cohorts of MSI-L and MSS patients with CRC (n=110). Dashed line at HR=1 indicates the numerical distance from no survival benefit. The cox proportional HR of the Vδ1-Cy-SI high and low-expression cohort was calculated by GEPIA2, while 95% CI was calculated as ‘exp [ln(HR)±z×SE], with SE’.51 Statistically significant p values are represented with the following number of asterisks (*): *p ≤ 0.05; ***p ≤ 0.001; ****p ≤ 0.0001; ns, not statistically significant.

On the other side, Vδ3 TILs are characterized by the high transcript levels of HLA-DR coded molecules associated with antigen processing and presentation. Moreover, Vδ3 TILs show high responsiveness to interferon-gamma (IFN-γ) in association with increased transcriptions of heat shock proteins (HSPs) that are, in turn, linked to biological pathways of protein stability and folding. These latter functional features mainly characterize c6 of Vδ3 TILs.

The GO pathways enriched at the highest levels among Vδ1 TILs encounter innate immune response, cytokine production and cell killing. In particular, Vδ1 TILs show great pro-inflammatory cytokine-secretion signature (Cy-SI) characterized by high expression of IFN-γ, TNF and X-C motif chemokine ligand 2 (XCL2) (figure 4A–C and online supplemental figure 2). Considering the positive prognostic impact of high frequencies of CD69+ Vδ1 TILs in the clinical outcomes of CRC and CLM (figure 1F–H),20 we evaluated whether the functional annotations of this specific Cy-SI among Vδ1 TILs played any role in limiting the clinical progression of metastatic CRC. To do so, we searched for patients with CRC in The Cancer Genome Atlas (https://cancergenome.nih.gov) datasets by using the Gene Expression Profiling Interactive Analysis 2 (GEPIA2) web server. In the absence of available metadata on a specific cohort of patients with CLM, we selected individuals affected by CRC at highest risk of liver metastases. These patients are characterized by low-frequency microsatellite instability or microsatellite stability and represent approximately 90% of all patients progressing toward CLM.34 35 Our results showed that the transcript levels of Vδ1 TILs in CRC tumor environment significantly correlate with the pro-inflammatory and antitumor Cy-SI (figure 4D). Importantly, higher Vδ1-Cy-SI predicts better disease-free survival in these patients. Moreover, we also found that patients with CRC ranked with a higher Vδ1 (TRDV1)-Cy-SI are associated with significantly longer OS at different cut-off values (figure 4E,F).

Distribution and transcriptional signatures of circulating and tumor-infiltrating γδ T cells from patients with CLM

We then performed an integrated scRNA-seq analysis to compare the distribution and transcriptional signatures of γδ T cells from matched PB and tumor specimens collected from patients undergone surgical resection of CLM. The differential expression of TRDV transcripts identified six clusters (c0–c5) of γδ T cells showing different repertoires of Vγ chains. We then excluded c4 from our further analyses, due to its low expression of δ and γ chains, respectively (figure 5A,B, online supplemental figure 3A).

Figure 5

Single cell RNA-sequencing (scRNA-seq) integrated analysis comparing the distribution of circulating and tumor-infiltrating gamma delta (γδ) T cells and transcriptional profiles of Vδ3 T cells from patients with colon liver metastatic cancer (CLM). (A) Uniform manifold approximation and projection (UMAP) graph showing the γδ T clusters from cells purified from matched peritumoral (PT) and peripheral blood (PB) samples of three patients with CLM undergoing surgical resection of tumors. (B) Heatmap showing the average of T cell receptor (TCR) δ chain expression along the six identified γδ T cell clusters showed in panel A. (C) UMAP graph (left panel) and bar plot graph (right panel) showing the distribution (%) of the γδ T cell clusters within PT (gray) and PB (red) anatomic compartments. (D) UMAP density plot graphs showing CD27 and CX3CR1 gene density distribution among γδ T cell clusters. (E) Heatmap showing the top 30 differently expressed genes (DEGs) (adjusted p<0.05) between PB and PT Vδ3 T cells in patients with CLM (left panel). Representative UMAP graph showing the distinct molecular signatures of PB and PT Vδ3 T cells (right panel).

The distribution of the different γδ T cell clusters within PB and CLM tissue specimens shows that c1, c2 and c3 are present in both anatomic compartments, while c0 and c5 resemble respectively two Vδ1 and Vδ3 T cell subsets infiltrating the tumor mass and retained in situ (figure 5C). We observed that c1 comprises high frequencies of PB Vδ3 T cells, a subset of γδ T cells that had been previously reported to be preferentially enriched in liver tissues under homeostatic conditions.18 Interestingly, the comparison of DEGs between PT-Vδ3 and PB-Vδ3 T cells in c1 revealed that the highest cytotoxic potential is held by the circulating ones displaying increased levels of GNLY, GZMB, PRF1 and of the activating receptors NKG7, NCR3 and FCGR3A. For more accurate immunophenotyping in CLM, we compared different γδ T cell subsets in relation to their CD27high/low transcriptional pattern. All Vδ3 T cells share CD27low transcriptional profile of mature T cells (figure 5D). In addition, as a peculiarity of Vδ3 T cells, their CD27low profile in the PB correlates with high expression of CX3CR1. The specific role of this marker in γδ T cells is still unknown. However, it has been recently proposed as a marker of differentiated and long-lived tissue resident T cells (figure 5D,E).36

Previous studies reported that the release of Vδ3 T cells in the bloodstream occurs preferentially in the presence of pathological conditions, including liver cancers.37 38 In this context, our data confirmed this phenomenon by showing that terminally differentiated Vδ3 T cells can egress from the liver to PB also in patients with CLM. Moreover, Vδ3 T in c5 can be easily discriminated due to their high transcript level of HSPs, thus confirming their liver-identity in the context of an immune response within the CLM microenvironment.

We then analyzed the transcriptional signatures of Vδ1 and Vδ2 TILs that are characterized by higher degrees of activation (eg, FOS, FOSB, JUN, JUNB-D, KLF6) and increased transcription of chemokine/cytokine transcripts compared with their PB counterparts. Moreover, as already shown in several cancer settings,39 Vδ1 and Vδ2 TILs express decreased levels of mitochondrial (mt)DNA-encoded genes (MT-CYB, MT-ATP6-8, MT-ND1-6, MT-CO1-3) compared with their circulating homologs (online supplemental figure 3B,C). A more in-depth analysis revealed that among the Vδ2 compartment in c2, there is a population of approximately 40% of cells able to infiltrate CLM (figure 5C). c2 is composed of both less mature CD27high and more differentiated CD27low Vδ2 T cells. These latter two populations show distinct transcriptional signatures, are mainly enriched in the PB compartment and share a similar capacity to infiltrate the tumor. The third subset comprised in c2 is composed by type 3 Vδ2 cells, a population that had been previously reported of being able to polarize into Th-17.33 These latter lymphocytes are mainly enriched within the PT compartments of CLM and express a distinct transcriptional signature, thus being tightly associated with tumor-associated activation within the CLM microenvironment (figure 5D, figure 6A,B).36

Figure 6

Single cell RNA-sequencing (scRNA-seq) integrated analysis comparing the transcriptional profiles of circulating and tumor infiltrating Vδ1 and Vδ2 T cells from patients with colon liver metastatic cancer (CLM). (A) Heatmap showing the top 30 differently expressed genes (DEGs) (adjusted p<0.05) between peripheral blood (PB) and peritumoral (PT) Vδ2 T cells from patients with CLM (left panel). Representative uniform manifold approximation and projection (UMAP) graph showing the distinct molecular signatures of PT-enriched type 3, PB CD27high and CD27low Vδ2 T cells (right panel). (B) Bar plot graph showing the distribution (%) of type 3, CD27high and CD27low Vδ2 T cells in PB (red) and PT (gray) anatomic compartments from patients with CLM. (C) Heatmap showing the top 30 DEGs (adjusted p<0.05) between PT and PB Vδ1 T cells from patients with CLM (left panel). Representative UMAP graph showing the distinct molecular signatures of the PT and PB CD69high and CD69low Vδ1 T cells subsets in patients with CLM (right panel). (D) Bar plot graph showing the distribution (%) of the CD69high and CD69low Vδ1 T cells in PB (red) and PT (gray) anatomic compartments from patients with CLM.

Vδ1 T cells are contained within c0 and c3 and all have in common an effector CD27low signature. While Vδ1 T cells in c0 show a high tissue-affinity resembling their liver residency, Vδ1 T cells in c3 are mainly enriched within the PB compartment (figure 5C,D). Our in-depth integrated analyses identified two subsets of circulating Vδ1 T cells showing different transcription levels of CD69. In particular, PB CD69high Vδ1 T cells express a ‘liver transcriptional signature’ characterized by higher expressions of CXCR6, CREM, RGS1, ZNF331, SRRT and GZMK compared with PB CD69low Vδ1 T cells. In addition, PB CD69high Vδ1 T cells show an activated profile (eg, increased levels of XCL2, JUNB, ZFP36, DUSP4 and FYN) and lower expression of (mt)DNA-encoded genes (MT-ND4-6, MT-CO2, MT-ATP6) compared with PB CD69low Vδ1 T cells (figure 6C,D), thus further confirming their liver tropism. Taken together, our integrated analyses strongly suggest that the subset of CD69+ Vδ1 T cells in CLM has the capacity to egress the liver and migrate to PB.

Circulating CD69+ TEMRA Vδ1 cells share similar phenotypic and TCR repertoires with their intrahepatic counterparts

Considering that Vδ1 TILs represent the main γδ T cell subset enriched in CLM specimens (figure 1C) and given the similarities in the transcriptomic signatures of CD69+ Vδ1 T cells in PB and tumor mass (figure 6C,D), we proceeded to assess the frequencies, phenotypes and TCR repertoires of circulating Vδ1 T cells in patients with CLM. Multiparametric flow cytometry results showed that the frequencies of PB Vδ1 T cells increase in patients with CLM compared with HDs. In particular, we found significantly higher percentages of TEMRA Vδ1 cells in the PB of patients with CLM compared with HDs (figure 7A,B). Interestingly, circulating Vδ1 T cells from patients with CLM express significantly higher surface levels of CD69 compared with their counterparts in HDs. This phenomenon is mainly due to the expansion of PB CD69+ TEMRA Vδ1 T cells that account up to 70% of total Vδ1 TEMRA cells in CLM (figure 7C,D). Moreover, the high frequencies of PB CD69+ TEMRA Vδ1 T cells significantly correlate with the enrichment of their cellular homologs in the tumor mass of the same patients with CLM (figure 7E).

Figure 7

Identification of intrahepatic Vδ1 tumor-infiltrating lymphocytes (TILs) egressing tumor and re-circulating in peripheral blood (PB) of patients with colon liver metastatic cancer (CLM). (A) Statistical bar graph showing the mean (±SEM) frequency (%) of PB Vδ1 T cells among CD3+ T lymphocytes in health donors (HDs) (n=66) and patients with CLM (n=75). (B) Pie charts showing the mean frequency distribution (%) of PB naïve T (TNAIVE), T central memory (TCM), T effector memory (TEM) and terminally differentiated (TEMRA) Vδ1 T cells in HDs (n=33) and in patients with CLM (n=50) (upper panel) either in the absence (n=13) (lower left) or in the presence (n=38) (lower right) of neoadjuvant conventional/biological chemotherapies (na-CHT) (lower panel). (C) Heatmaps of multiparametric flow cytometry data showing the expression of several surface markers between PB Vδ1 T cells isolated from patients with CLM (left) and HDs (n=49) (right). Group of markers which statistically differ between two groups are annotated under the heatmap as ‘‘≥** and indicates the range of significant p values; ‘ns’ indicates a group of markers not statistically significant. (D) Statistical bar graph showing the mean (±SEM) frequency (%) of PB CD69+ TEMRA Vδ1 cell subset in HDs (n=22) and patients with CLM (n=45). (E) Pearson’s correlation of flow cytometry percentages (%) between matched PB CD69+ TEMRA and PT CD69+ TEMRA Vδ1 TILs (n=45). (F) Statistical bar graph showing the mean (±SEM) frequency (%) of PB CD69+ TEMRA Vδ1 cells from patients with CLM either in the absence (n=11) or in the presence (n=45) of na-CHT. (G) Heatmap of flow cytometry data showing Spearman’s rank between the mean of different cell markers expression (%) in PB Vδ1 T cells from patients with CLM (n=49). (H) Overlap analysis of T cell receptor (TCR) repertoires from FACS-sorted PB CD69+ and CD69 and liver PT Vδ1 T cells in patients with CLM (n=5). Shared top 20 TRG clones of one representative patient with CLM are represented as colored bands between columns (left panel) and included in the summary statistic (right plot). (I) Kaplan-Meier curve showing postoperative OS (%) of patients with CLM stratified on the basis of the median frequency (%) of the PB CD69+CD28 Vδ1 T cells in patients with CLM (cut-off ≤10%; n=53). Statistically significant p values are represented with the following number of asterisks (*): *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001; ns, not statistically significant.

We previously reported that toxicity of na-CHT administered prior surgery to elder patients with CLM increases circulating Vδ2 T cells with highly differentiated TEMRA phenotype (online supplemental figure 4A).24 To assess whether a similar phenomenon occurs also for PB Vδ1 T cells, we divided our cohort in two subgroups either undergone or not to conventional na-CHT and immunotherapies (table 1). Surprisingly, we did not observe any significant differences, in terms of frequencies and phenotype of circulating Vδ1 T cells between the two subgroups of patients with CLM (figure 7B and online supplemental figure 4B). Likewise, different biological agents used in combination with na-CHT (table 1) showed no effect on the frequencies of PB CD69+ TEMRA Vδ1 cells (figure 7F). Importantly, our data also showed that the frequencies of CD69+ TEMRA Vδ1 cells are independent from age, sex and human cytomegalovirus status of the enrolled patients with CLM (online supplemental figure 4D,E).24 37 40 Even more interestingly, similar to CD69+ Vδ1 TILs in the liver, the expression of CD69 in PB Vδ1 T cells of patients with CLM displays phenotype which correlates with the expression of CD56 and CD161 and inversely associates with CD28 and CD57 in the entire patients with CLM (figure 7G). Indeed, t-distributed Stochastic Neighbor Embedding (t-SNE) analysis shows striking phenotyping overlap of PB CD69+ Vδ1 TEMRA phenotype with the expression of CD56 and CD161 and low expression of CD28 (online supplemental figure 4F,G).

We then compared the TCR repertoires of FACS-sorted CD69 and CD69+ Vδ1 TILs with the ones of their PB homologs from the same patients with CLM. To this end, we performed high-throughput sequencing of V(D)J regions that paired the Vδ1 chain with different γ chains (online supplemental figure 4H,I). Although analyses of cell clonality revealed high interchanges of Vδ1 T cell clones between PB and CLM, PB CD69+ Vδ1 T cells share a significantly higher number of Vδ1 clones with their CLM counterparts when compared with PB CD69 Vδ1 T cells (figure 7H). These data demonstrate that PB CD69+ Vδ1 T cells present in the blood of patients with CLM include Vδ1 TIL clones egressing from tumor and keeping their high expression of the tissue-residency marker CD69 within the circulating TEMRA compartment.

In a more clinical context, we then observed that higher frequencies of PB CD69+CD28 Vδ1 T cells are associated with longer OS of patients with CLM (figure 7I). These latter findings mirror the positive prognostic impact on OS of higher frequencies of CD69+CD28 Vδ1 TILs in CLM specimens detected at the time of surgery (figure 1H).

Discussion

The present study reveals the heterogeneous distribution, phenotypes, functional relevance and prognostic values of γδ T cells either infiltrating CLM or re-circulating in PB. Our results showed that CLM tissue specimens comprise all main subsets of γδ TILs including Vδ3 T cells, a subset naturally resident in healthy human liver without egressing the organ under physiological conditions.17 18 scRNA-seq results identified seven γδ TIL clusters enriched in CLM specimens and expressing distinct transcriptional signatures and functional annotations. In particular, we found that CLM microenvironment is enriched in type 3 Vδ2 T cells, a subset known for their high plasticity and potential to polarize into pro-inflammatory γδ Th-17 cells.33 40 Although γδ Th-17 cells have been reported to promote cancer development,41 secretion of IL-17A by human Vδ2 T cells is still a matter of discussion. In our study, we did not observe any IL-17A production at transcriptional level that, instead, revealed a preferential production of TNF-α by type 3 Vδ2 TILs. Hence, additional investigations are required to disclose the impact of type 3 Vδ2 T cells in the natural history of CLM. Our findings also disclosed the molecular signature of mature cytotoxic Vδ3 T cells, whose activation and response to CLM is associated with a strong stress-HSP-induced transcriptional profile. The integrated analyses of scRNA-seq data from matched PB and CLM specimens also revealed the capacity of intrahepatic Vδ1 and Vδ3 TIL clusters to egress tumor and re-circulate in PB.

Vδ1 T cells represent the largest fraction of γδ TILs enriched in CLM and highly impact the clinical outcomes of these patients. Indeed, we demonstrate here that liver CD69+ Vδ1 T cells highly infiltrate CLM and limit the metastatic progression of CRC to liver. In specific, our data show an inverse correlation between the percentage of TEF CD69+ Vδ1 TILs and the number of metastatic lesions in the liver that, in turn, predicts longer patients’ OS after surgery. These findings indicate the existence of local CD69+ Vδ1 TIL-mediated antitumor response at tissue sites. In line with recently reported data during liver inflammatory diseases,18 these anticancer activities in patients with CLM are exerted by two subsets of the liver CD69+ TEF Vδ1 cell: the TEM Vδ1 cells that reside permanently in the liver tissue and the TEMRA Vδ1 cells that can re-circulate toward PB.

The ability of Vδ1 TILs to control CRC progression is further confirmed by their strong pro-inflammatory Cy-SI characterized by high expression of IFN-γ, TNF and XCL2 compared with other γδ TIL subsets. Certainly, future studies are necessary to evaluate whether these properties of Vδ1 TILs can be translated to the CLM clinic, nevertheless, our data support a model in which Vδ1 TILs create large cytokine field able to penetrate tumor and whose concentrations depend on the collective Vδ1 T cell activation. In this context, recently it was observed that while outset of IFN-γ secretion requires local T cell response, prolonged activation diffuses IFN-γ to distant areas and extensively alters the tumor phenotype.42

CD69 is a C-type lectin receptor widely used to track immune cell activation and tissue-resident T lymphocytes that persist throughout adult life. Indeed, after entering tissues, all T lymphocytes upregulate expression of CD69 which, in turn, promotes cell retention into peripheral organs.29 The expression of CD69 on tissue T lymphocytes is constitutive and not associated with degrees of activation as its surface levels do not correlate with the ones of activation markers such as CD25 or CD38.43 However, the functional relevance of CD69 on immune cells at tissue sites is still being discussed.29 We characterize here a distinct subset of CD69+ TEMRA Vδ1 TILs that can egress from CLM to re-circulate in the bloodstream where it retains a phenotypic/transcriptional profiles resembling its liver tropism and origin. In this regard, all γδ T cells residing in tissues constitutively express CD69 together with specific repertoires of homing receptors different from those of their circulating counterparts.18 44 45 The re-circulation of Vδ1 TILs in PB of the same patients with CLM is further confirmed by TCR sequencing experiments showing that these two anatomic compartments share high frequencies of identical CD69+ Vδ1 T clones. In addition, scRNA-seq analysis identified in the PB of patients with CLM the existence of two CD27low Vδ1 T cells with distinct transcriptional profiles that include high expression levels of CD69. Hence, the ability of intrahepatic CD69+ Vδ1 TILs to migrate toward PB defines de facto a novel liver-blood barrier in CLM. The mechanisms regulating the re-circulation of intrahepatic CD69+ TEMRA TILs in patients with CLM are unclear. Moreover, the extent of clonal exchange between PB-PT compartments should be further investigated for Vδ3 T cells that were detected at high frequencies in the PB of patients with CLM. The inducible expression of CX3CR1 might contribute to explain the egressing of Vδ3 TILs in CLM considering that CX3CR1+CD8+ T cells were recently identified as tissue peripheral memory cells with unique migratory and self-renewal properties.36 In addition, expression of CX3CR1 was recently proposed as a blood-based T-cell biomarker in response to cancer treatment.46 Interestingly, it has been reported that skin-resident CD4+ T cells with a distinct memory tissue-signature can egress to PB even under homeostatic conditions.47 Therefore, further investigations are required to better understand the homing of long-term persistent liver Vδ1 and Vδ3 T cells.

Higher frequencies of both PB and TI CD69+ TEMRA Vδ1 cells predict a longer OS of patients with CLM. This phenomenon represents a potent and reliable prognostic marker immediately transferable to everyday clinical practice through the implementation of liquid biopsies without the need to always analyze liver tissue. Our gained knowledge can remarkably improve our capacity to predict tumor progression or response to therapies. Indeed, the identification of novel prognostic biomarkers in CLM is hampered by the fact that the majority of patients receive na-CHT. In this regard, the toxic side effects of these therapeutic regimes induce cell death and induce proliferative-senescence, as we observed in Vδ2 T cells.24 This is not the case of Vδ1 T cells. Indeed, despite the many na-CHT options available in CLM, the frequency and the phenotype of mature PB Vδ1 cells are not altered or affected by drug toxicity. This implies that circulating CD28CD69+ Vδ1 cells are valuable prognostic markers predicting CLM clinical outcome regardless of the administration of na-CHT.

Considering the adaptive nature of human Vδ1 T cells,48 the identification of specific clones of circulating antitumor Vδ1 T cells is also key for developing new protocols of adoptive cell transfer therapies. Up to date, the vast majority of currently active clinical trials relies on the transfer of in vitro expanded Vδ2 T cells or polyclonal γδ T cells.10 11 CD69+ TEF Vδ1 cells represent the perfect target for setting novel protocols of cellular immunotherapy in CLM, as this subset of lymphocytes egresses from tumor sites, is endowed with high antitumor potential, is present at high frequencies both in PB of patients with CLM and is expandable in vitro.49 50

Conclusion

The tumor microenvironment of CLM is highly enriched in a heterogenous infiltrate of γδ TIL subsets with potent antitumor potentials. The present study identifies a specific subset of liver CD69+ Vδ1 TEMRA cell able to egress and re-circulate in the PB. The higher are its frequencies both at tissue site and in PB, the longer is patients’ OS regardless of na-CHT administered to patients with CLM. Taken together, these findings demonstrate that CD69+ Vδ1 TEF cells in CLM specimens limit the metastatic progression of CRC to liver and can also serve as a reliable prognostic marker in ‘liquid biopsies’ to identify those patients with better clinical outcomes.

Supplemental material

Data availability statement

Data are available on reasonable request.

Ethics statements

Patient consent for publication

Ethics approval

This study was approved by Institutional Review Board (IRB) of Istituto Clinico Humanitas (approval 168/18). Participants gave informed consent to participate in the study before taking part.

Acknowledgments

We thank all the patients who participated in this study.

References

Supplementary materials

Footnotes

  • Contributors EB performed experiments, analyzed data and wrote the manuscript; MMC recruited patients, provided reagents/specimen and analyzed clinical data; MD and GT recruited patients and provided reagents/specimen; RC, ST and RP performed scRNA-seq analysis; SR and IP designed and performed γδ TCR-seq analysis; VC and PM performed experiments; PK, ST and CG provided assistance with the scRNA-seq experiments; CS and BF performed immunohistochemistry study; FSC provided assistance with flow cytometry and cell sorting; CG and AM provided reagents and assistance in the manuscript’s preparation; JM and DM directed and designed the study, analyzed data and wrote the manuscript. J.M and D.M are the guarantors of this study.

  • Funding This work was supported by Associazione Italiana per la Ricerca sul Cancro (IG 14687 to DM and 5 X 1000-9962 to AM), Italian Ministry of Health (Bando Ricerca Finalizzata PE-2016-02363915 to DM and RF-2018-12367150 to MD), intramural research and clinical funding programs including Humanitas Research Hospital (5 X 1000 to DM and MD), University of Milan (to DM) and FOR2799 and RESIST (to SR and IP). EB and VC are recipients of competitive fellowships awarded from the PhD program of Experimental Medicine from University of Milan. ST is a recipient of competitive fellowship awarded from the Data Science in Medicine and Nutrition (DASMEN) PhD program from Humanitas University. The purchase of a FACS Symphony A5 was defrayed in part by a grant from the Italian Ministry of Health (Agreement 82/2015).

  • Competing interests None declared.

  • Provenance and peer review Not commissioned; externally peer reviewed.

  • Data ara available on reasonable request from the corresponding author D. Mavilio [domenico.mavilio@humanitas.it] All data relevant for the study are included in the article or uploaded as supplementary information. The raw data from single cell RNA sequencing and TCR sequencing are available on reasonable request from the corresponding author D. Mavilio [domenico.mavilio@humanitas.it]

  • Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.