Assessment of tumor-infiltrating TCRVγ9Vδ2 γδ lymphocyte abundance by deconvolution of human cancers microarrays

ABSTRACT

Most human blood γδ cells are cytolytic TCRVγ9Vδ2⁺ lymphocytes with antitumor activity. They are currently investigated in several clinical trials of cancer immunotherapy but so far, their tumor infiltration has not been systematically explored across human cancers. Novel algorithms allowing the deconvolution of bulk tumor transcriptomes to find the relative proportions of infiltrating leucocytes, such as CIBERSORT, should be appropriate for this aim but in practice they fail to accurately recognize γδ T lymphocytes. Here, by implementing machine learning from microarray data, we first improved the computational identification of blood-derived TCRVγ9Vδ2⁺ γδ lymphocytes and then applied this strategy to assess their abundance as tumor infiltrating lymphocytes (γδ TIL) in ∼10,000 cancer biopsies from 50 types of hematological and solid malignancies. We observed considerable inter-individual variation of TCRVγ9Vδ2⁺γδ TIL abundance both within each type and across the spectrum of cancers tested. We report their prominence in B cell-acute lymphoblastic leukemia (B-ALL), acute promyelocytic leukemia (M3-AML) and chronic myeloid leukemia (CML) as well as in inflammatory breast, prostate, esophagus, pancreas and lung carcinoma. Across all cancers, the abundance of αβ TILs and TCRVγ9Vδ2⁺ γδ TILs did not correlate. αβ TIL abundance paralleled the mutational load of tumors and positively correlated with inflammation, infiltration of monocytes, macrophages and dendritic cells (DC), antigen processing and presentation, and cytolytic activity, in line with an association with a favorable outcome. In contrast, the abundance of TCRVγ9Vδ2⁺ γδ TILs did not correlate with these hallmarks and was variably associated with outcome, suggesting that distinct contexts underlie TCRVγ9Vδ2⁺ γδ TIL and αβ TIL mobilizations in cancer.

KEYWORDS:

• Artificial intelligence
• cancer
• data mining
• deconvolution
• gamma delta lymphocyte
• machine learning
• microarray
• transcriptome

Introduction

In healthy human adults, circulating T lymphocytes comprise ∼1–3% of γδ T cells, most of which express TCRVγ9Vδ2, and mediate Th1 and cytolytic activities in response to stimulation by soluble phosphoantigens (PAg) or cancer cells.^1-3 Hence, this subset of T lymphocytes represent promising effectors for cancer immunotherapy, but the few γδ cell-based trials performed so far have proven difficult to interpret due to a lack of biomarkers for the stratification of patients.^4-7 Notably, the detection of tumor infiltrating γδ T lymphocytes (γδ TIL) has been hampered by their controversial detection by immuno-histochemical (IHC) methods due to their unreliable recognition by antibodies.

Alternative options would be to enumerate γδ TIL from gene expression profiles (GEP), which is now possible due to recently developed computational methodologies that can perform cell type-specific quantifications from the transcriptomic data of bulk samples (for review, see⁸). For example, the “deep deconvolution” CIBERSORT algorithm uses the transcriptome profile of reference leucocyte cell types to define a reference gene signature matrix which it then uses to deduce the leucocyte composition from the transcriptome of a tissue.⁹ Using CIBERSORT, the Alizadeh group were able to determine the relative proportions of 22 subsets of leucocytes that had infiltrated tumors from 30,000 transcriptomes representing 39 distinct malignancies, and identified γδ TILs as the most significant favorable cancer prognostic cell population.¹⁰ However, we show that applying CIBERSORT to pure γδ T cell transcriptomes failed to correctly differentiate these cells from αβ T CD4⁺, CD8⁺ and natural killer (NK) cells. Here we optimized this deconvolution and applied it to enumerate γδ TIL in human cancers.

Results

Optimization of CIBERSORT deconvolution for inferring γδ T cell abundance from a transcriptome

We first applied the CIBERSORT algorithm and its LM22 gene signature to determine the fraction of γδ T cells from the GEP of 12 well-characterized samples of TCRVγ9Vδ2⁺ γδ lymphocytes purified from healthy donor peripheral blood mononuclear cells (PBMC) (purity >95%, assessed by staining for the TCRVγ9 and TCRVδ2 chains followed by FACS analysis).¹¹ CIBERSORT-LM22 detected 7–40% of γδ cells and overestimated the CD4⁺, CD8⁺ and NK cell content of these samples (Fig. 1A). In addition, CIBERSORT-LM22 enumerated up to 55% of γδ T cells in transcriptomes from purified CD4⁺, CD8⁺ and NK cells (Table S1). The GEP of these related cell subsets partly overlap¹¹ and the multicollinearity of this overlap is presumably what drives CIBERSORT to wrongly identify γδ T cells as T CD4⁺, T CD8⁺ or NK cells and vice versa. Thus, these data demonstrated the limits of CIBERSORT's LM22 reference matrix for the rating of TCRVγ9Vδ2⁺ γδ cells. The LM22 reference gene signature matrix was developed by machine learning from Affymetrix HG-U133A microarrays (∼8500 genes) from 113 reference samples downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/) and which included only two samples annotated as “blood γδ cell” whose purity and TCR phenotype, although most likely of the TCRVγ9Vδ2⁺ subset, were not specified.¹² Furthermore, of the 58 genes in the LM22 gene signature which are aimed at identifying γδ cells, the only TCRγδ-specific gene TRDC was actually not expressed by any of the 12 purified TCRVγ9Vδ2⁺ γδ cell samples described above. These observations prompted us to seek to improve the deconvolution of tissue expression profiles for assessing TCRVγ9Vδ2⁺ lymphocytes. We undertook this by increasing the CIBERSORT learning phase and reducing the number of cell subsets to discriminate between. The group of reference data sets was increased to 265 profiles, taken only from the Affymetrix HGU133 Plus 2.0 platform (which analyzes twice more genes than HGU133A microarrays), and included the above-depicted 12 samples of purified TCRVγ9Vδ2⁺ γδ T lymphocytes (Table S2). To optimize the CIBERSORT output, the number of subsets to determine was reduced from 22 to 7. These included B cells (pooling naïve, memory B cells and plasma cells), CD4⁺ cells (pooling CD4⁺ resting, CD4⁺ activated, regulatory T cells (Treg), and follicular helper T cells (TFH), CD8⁺ cells (resting and activated CD8⁺ cells), γδ T cells (resting and activated TCRVγ9Vδ2⁺ γδ cells), NK cells (resting and activated NK cells), “MoMaDC” cells (pooling monocytes, macrophages and DC), and granulocytes (arbitrarily pooling mastocytes, eosinophils and neutrophils). These 265 reference data sets were curated, assembled, normalized using the robust multiarray average (RMA) method and then collapsed to construct a novel leucocyte gene signature matrix using the CIBERSORT learning algorithm.⁹ After the removal of genes expressed by non-hematopoietic tissues the resulting matrix, called LM7 (Table S3), comprised 375 genes, 134 of which were present in LM22. This learning identified 23 genes specifically upregulated in human γδ T cells only: BCOR, BTBD5, CCR5, CREM, CXCR3, CXCR6, FAM46C, FLJ25967, GPR132, GPR171, LEPREL1, LOC653506, LOC653725, METRNL, NR4A3, PDE4B, PDE4D, PRKX, PRKY, RORA, SYTL2, SYTL3 and UPP1, while it found that the TCRγδ−encoding genes TRD@, TRDV2, TRGC2, TRGV2 and TRGV9 are strongly expressed in both γδ T, αβ T and NK cells (Table S3).

Figure 1. Deconvolution of microarrays from pure TCRVγ9Vδ2⁺ γδ T lymphocytes. (A) CIBERSORT detected only marginal rates of γδ cells when the original reference matrix LM22 was applied to deconvolution of the transcriptomes from 12 samples of purified (> 90% TCRVγ9Vδ2⁺)γδ T lymphocytes. (B) Using the reference matrix LM7, CIBERSORT detected respectively 92%, 76%, 92%, 93%, 87%, 71%, 94%, 84%, 95%, 85%, 81% and 85% of TCRVγ9Vδ2⁺ γδ cells in the same samples as above. (C) Deviation of CIBERSORT results from the rates specified in each sample annotation file. With LM7, the deconvolution results were closer to most rates specified in the annotation file of control samples of known leucocyte composition, compared with LM22. (D) Pearson correlation of results from CIBERSORT-LM22 or CIBERSORT-LM7 with FACS-based reference rates of CD4⁺, CD8⁺ and B cells in a series of follicular lymphoma biopsies (GSE 65135). (E) Results from CIBERSORT-LM22 and CIBERSORT-LM7 yielded consistent rates of αβ TILs but not of γδ TILs upon the deconvolution of ∼10,000 cancer data sets (see below). The rates of αβ TILs were taken as the sum of CD4⁺ and CD8⁺ cells. Pearson correlation values are specified; the smaller scale of the γδ TILs scatterplot reflects their lower abundance than αβ TIL in the biopsies. (F) Mean leucocyte abundance in normal tissues from a validation series of (n = 114) transcriptomes.

With the above reference data set (n = 265 samples), the relative leucocyte proportions predicted by CIBERSORT-LM7 and CIBERSORT-LM22 were concordant with the compositions annotated in most GEO sample record (GSM), with the exception of the above-depicted TCRVγ9Vδ2⁺ γδ samples whose values were still perfectible though significantly improved with LM7 (which calculated >90% γδ cells on average, versus 31% with LM22) (Fig. 1B and Table S1). We then downloaded from GEO and processed a validation data set of tissues (n = 114) whose composition had previously been determined (Table S4). Deconvolution of this data set using CIBERSORT-LM7 and CIBERSORT-LM22 matched the specified fractions of B cells and CD4⁺ cells, while the predicted levels of cytolytic (CD8⁺, NK, TCRVγ9Vδ2⁺ γδ), myeloid (monocytes, macrophages and DC) and granulocyte subsets improved with LM7 (Fig. 1C, D). However, deconvolution with CIBERSORT-LM7 of (n = 6) transcriptomes from purified cord blood γδ T cells, which encompass both TCRVγ9Vδ2 and non-TCRVγ9Vδ2 γδ T cells dominated by TCRVγ8Vδ1 cells in cytomegalovirus-infected samples,¹³ gave quite low results (detection of 24% γδ cells on average; range 11–51%), which were better than with CIBERSORT-LM22 (detection of < 5% γδ cells on average). Further, deconvolution of the transcriptomes of ∼10,000 tumor biopsies of unknown composition (section below) by CIBERSORT-LM22 and CIBERSORT-LM7 was concordant for αβ tumor-infiltrating lymphocytes (TILs) (taken as sum of CD4⁺ and CD8⁺ T cells, Pearson r = 0.72) but not for γδ TILs (Pearson r = −0.12) (Fig. 1E), and was variably concordant for the other subsets: Mo-Ma-DC (Pearson r = 0.8), B lymphocytes (Pearson r = 0.76), NK (Pearson r = 0.3) and Granulocytes (Pearson r = 0.2) (data not shown).

Although CIBERSORT-LM7 was able to determine the leucocyte composition of tissues arrayed with Affymetrix HG-U133 Plus 2.0 microarrays, it was possible that the use of other platforms which lack some of the LM7 genes could alter its results. We estimated this bias by evaluating the LM7-generated leucocyte subset fractions from the training data set with and without the genes missing from the other platforms (Fig. S1) and found that the incomplete data introduced negligible (< 5%) deviations from the original deconvolution results (Fig. S2).

CIBERSORT determines the relative fractions of leucocytes within a sample, and not their abundance, whereas the sample enrichment score (SES) of leucocytes (SES-LM7) (see Materials and methods) reflects the abundance of leucocytes in each sample¹⁴ and parallels CD45 gene expression levels (Fig. S3). Thus, the abundance of each leucocyte subset in each sample was obtained by computing [(SES-LM7) × (CIBERSORT-LM7 fraction)] using the DeepTIL algorithm (see Materials and methods) which yields abundance results expressed as arbitrary units. The leucocyte abundance estimates determined by this method from HGU133 Plus 2.0 microarrays of normal tissues (n = 494) were consistent with the annotations of purified cell samples and the composition of normal tissues in that leucocytes were found to be most abundant in PBMC, lymph nodes and lung tissue while scarce in normal brain and kidney tissue (Fig. 1F). This validated the DeepTIL algorithm as a useful tool for obtaining more accurate quantification of leucocyte subset abundance, particularly the blood-derived TCRVγ9Vδ2⁺ γδ T cells.

Abundance of tumor infiltrating leucocytes in a spectrum of cancers

We then investigated the abundance of TCRVγ9Vδ2⁺ γδ T cells and other LM7-defined leucocyte subsets determined by this computational method from a compendium of ∼10,000 bulk tumor biopsies from 16 hematological and 30 solid types of cancers (listed in Table S6). This showed that the hematological malignancies had a distinct leucocyte abundance and composition dominated by the cell-of-origin lineage (COO), except for Hodgkin's lymphoma (HL) and T-cell/histiocyte rich B cell lymphoma (THRBL) whose composition reflected the tumor microenvironment.⁹ In follicular lymphoma (FL), CIBERSORT-LM7 detected less TCRVγ9Vδ2⁺ γδ TILs (on average 0.22 a.u. corresponding to 0.5% of leucocytes) than previously reported (∼5% of leucocytes infiltrating FL and diffuse large B cell lymphoma (DLBCL) according to¹⁰). Our result (Fig. 2) are in agreement with earlier IHC and FACS studies showing significantly less γδ T cells in FL lymph nodes than in inflammatory adenitis lymph nodes.¹⁵ Of note, deconvolution of the same FL and DLBCL transcriptomes by CIBERSORT-LM-7 detected 10-fold less γδ cells than CIBERSORT-LM22, illustrating its more stringent identification of γδ cells (Fig. 1C, E). Besides a few cases with a γδ COO,^16,17 most peripheral T cell lymphoma (PTCL) and the other T cell malignancies did not present abundant TCRVγ9Vδ2⁺ γδ cell infiltrates, whereas the cohorts of B-ALL, CML and AML included in this study revealed an abundant TCRVγ9Vδ2⁺  γδ TIL compartment (corresponding to ∼10% of leucocytes). This result is higher than has previously been reported for AML samples,^9,18,19 and even for normal bone marrow (∼2.5% of leucocytes,²⁰). Among the eight subtypes of AML defined by the French–American–British (FAB) classification, acute megakaryoblastic leukemia (FAB-M7 subtype of AML) presented the greatest infiltration of total leucocytes and acute promyelocytic leukemia (FAB-M3 subtype of AML) presented the highest abundance of TCRVγ9Vδ2⁺ γδ TILs. Solid tumors generally presented less abundant leucocyte infiltration than blood malignancies, and half of their infiltrated content was consistently of myeloid lineage. Quantitatively, these infiltrates were found to be most important in lung, pancreas, prostate, breast, stomach, renal (papillary and clear cell subtypes), colorectal, bladder, liver, esophagus, uterus, liver, head and neck carcinoma, melanoma, thyroid, uterus, endometrium, cervix, and ovary carcinoma tissues (ranked by decreasing magnitude). The least abundant infiltrates were found in cancers carrying the lowest number of somatic mutations, such as kidney chromophobe carcinoma, Ewing sarcoma, kidney oncocytoma and indeed in central nervous system (CNS) cancers. The infiltrates of brain cancers are mostly of myeloid lineage, probably the monocyte-related microglia.²¹ The abundance of  TCRVγ9Vδ2⁺ γδ TILs varied across the different cancers, being quite high in prostate, esophagus, cervix, head and neck, lung, pancreas and breast carcinoma. Both CD8⁺ and TCRVγ9Vδ2⁺ γδ TILs were abundant in inflammatory breast carcinoma, whose samples generally consisted of triple negative breast cancers. In melanoma, both αβ and TCRVγ9Vδ2⁺ γδ TILs were also found to be more abundant in primary lesions than in metastases (Fig. 2). Overall, the extent of leucocyte infiltration was quite heterogeneous across cancers and the average abundance of αβ and TCRVγ9Vδ2⁺ γδ TILs varied likewise between cancer types (Table S5). This large inter-individual fluctuation (ranging across several orders of magnitude) was also seen in FACS analysis of TILs from small cohorts of AML and prostate cancer patients (not shown). Plotting the individual abundance of αβ TILs at the level of patients and ranking cancer types on this criterion confirmed that the abundance of immune infiltrates fluctuated strongly among patients both inside and across cancer types. The abundance of TCRVγ9Vδ2⁺ αβ TILs among patients was also found to be highly variable, though they were consistently less abundant than αβ TILs. Strikingly however, the individual abundance of TCRVγ9Vδ2⁺ γδ TILs did not appear to correlate with the abundance of the αβ TILs neither within nor across the different types of cancer (Fig. 3).

Figure 2. Abundance of TILs across 55 types of hematological and solid cancers. Shown are the mean leucocyte composition of each specified cancer type. Blood malignancies were found to be as “leucocyte-rich” as normal PBMC or lymph nodes, whereas solid tumors included less leucocytes. Lung carcinomas displayed the most abundant leucocyte infiltrates, whereas brain tumors had a particularly low leucocyte composition.

Figure 3. Individual abundances of tumor-infiltrating αβ T and TCRVγ9Vδ2+ γδ T lymphocytes across 55 types of hematological and solid cancers. Individual abundances of αβ and TCRVγ9Vδ2+ γδ TILs across cancer types. Blood and solid cancers are successively presented in order of increasing αβ TIL abundance. Each dot is a sample, red bars show the median of each cancer type, blue dotted lines show the mean abundance of the specified TILs across hematological (top) and solid (bottom) cancers. For each type of cancer, the Wilcoxon test was used to calculate the statistical significance of deviation from the mean abundance of TIL in hematological or solid cancers. The Benjamini–Hochberg (BH) correction was applied to determine the FDR for multiple comparisons, only the significant BH-corrected p-values are specified (right).

No correlation between αβ TIL and TCRVγ9Vδ2+ γδ TIL abundance across the cancer spectrum

The above results suggested that tumor infiltration by the αβ and TCRVγ9Vδ2⁺ γδ TILs occur in distinct immunological contexts. To characterize these contexts, we computed SES¹⁴ of all the gene sets from the MSigDB²² for the 10,000 cancer samples, and determined whether these correlated with the abundance of αβ and TCRVγ9Vδ2⁺ γδ TILs. This approach revealed significant differences between the TIL subsets. Across all cancers, the abundance of αβ TILs positively correlated with the enrichment of “TCR signaling” (Spearman r = 0.78), “Toll-like receptor signaling pathway” (Spearman r = 0.45) and related gene sets, whereas the abundance of TCRVγ9Vδ2⁺ γδ TILs did not correlate with these (Spearman r = −0.02 and −0.03, respectively). These distinct features were best exemplified within specific cancer types. For example, in colorectal carcinoma, the abundance of αβ TILs correlated with the enrichment of the “Response to IFNα”, “Inflammatory response”, “Antigen processing and presentation”, “T cell activation” and “Cytolytic activity” gene sets, as well as with the abundance of Mo-Ma-DC cells (Fig. 4A). We found that all of these metrics also cross-correlated, and patients with “αβ TIL-high” tumors presented significantly higher overall survival (OS) rates than those with “αβ TIL-low” tumors (Fig. 4B). This result is in accordance with the well-known association of αβ TIL density with the outcome of colorectal cancer patients.^23,24 In contrast, the abundance of TCRVγ9Vδ2⁺ γδ TILs did not correlate with any of these transcriptomic readouts (Spearman r < 0.12), since correlation with these reference MSIgDB gene sets was the same as with a collection of random gene sets (RGS, data not shown). Nevertheless, colorectal cancer patients with “TCRVγ9Vδ2⁺ γδ TIL-high” tumors presented better OS rates than those with “TCRVγ9Vδ2⁺ γδ TIL-low” tumors, irrespective of their abundance of αβ TIL. Hence, the presence of TCRVγ9Vδ2⁺ γδ TIL was like αβ TIL, associated with a favorable outcome in colorectal cancer, even though these lymphocytes are recruited through a distinct immunological context. The same conclusions were drawn from the analyses of additional cohorts of CLL, AML and prostate cancer patients (Fig. 4B). These and the above results (Fig. 3) indicate that the recruitment of αβ and TCRVγ9Vδ2⁺ γδ TILs occurs in different immunological contexts.

Figure 4. Different correlates of αβ and TCRVγ9Vδ2⁺ γδ TIL abundance. (A) In the specified groups of cancer patients, heatmaps of spearman correlation between abundance of αβ TILs (top) or TCRVγ9Vδ2⁺ γδ TILs (bottom) with SES for the gene sets indicated. The abundance of αβ TILs never shows the same correlation profile as the abundance of TCRVγ9Vδ2⁺ γδ TILs. (B) Kaplan–Meier plots for the overall survival of colorectal carcinoma, CLL and AML patients or for the biochemical relapse (BCR) of prostate cancer patients, according to the abundance of the αβ TILs only (left), the TCRVγ9Vδ2⁺ γδ TILs superimposed on the αβ TILs-low only (middle) and the TCRVγ9Vδ2⁺ γδ TILs on αβ TILs-high only (right). Abundance of both subsets of TILs is associated to favorable outcome in these cohorts of cancer patients.

Discussion

We report here the application of an emerging method of machine learning for microarray deconvolution to assess the abundance of leucocyte infiltration in tumor biopsies. Since its first application in the deconvolution of blood transcriptomes from patients with systemic lupus erythematosus through iterative linear least-squares fitting,²⁵ this field has rapidly expanded and has progressively validated the use of computational deconvolution to identify the individual component cells from RNA mixtures. The current in silico approaches for characterizing cellular heterogeneity from GEP encompass either the analysis of cell type-specific genes and gene set enrichments, or the direct algorithmic deconvolution of transcriptomes to deduce cell-type compositions (reviewed in⁸). These both now yield results for the rapid profiling of TILs that match those achieved by conventional flow cytometry (FACS) and IHC techniques, but without the biases due to the mechanical dissociation of samples or the antibodies used for reliable detection of γδ TILs by IHC of formalin-fixed paraffin-embedded tumor biopsies. Optimizing the discrimination of ontologically related cell types is not trivial, however, as it requires the multicollinearity of their highly similar transcriptomes be addressed. Several techniques of automated deconvolution have recently been developed toward this aim (reviewed in²⁶), among which the most robust and performant is currently the machine learning algorithm CIBERSORT.⁹ CIBERSORT applies nu-support vector regression to enumerate 22 leucocyte subsets from bulk tumor biopsies, and has allowed the identification of genes and leucocyte predictors of survival in a meta-analysis of thousands of cancers.¹⁰ More recently, CIBERSORT also enabled the successful leucocyte enumeration by deconvolution of microarray data sets from PBMC,²⁷ glioma,²⁸ lung carcinoma²⁹ and triple negative breast cancers.³⁰ Here, however, we found that CIBERSORT-LM22 does not discriminate between γδ T cells and the closely related T CD4⁺, T CD8⁺ and NK lymphoid subsets due to insufficient training from only two γδ Τ cell transcriptomes assessed on first-generation microarrays. The improved learning from more than 260 second-generation microarrays (encompassing ∼21,000 genes) and the reduced output of only 7 instead of 22 subtypes of leucocytes was expected to reduce the multicollinearity of their individual gene profiles, and thereby improve output reliability. Our report shows that CIBERSORT-LM7 detects TCRVγ9Vδ2⁺ γδ T lymphocytes more accurately than CIBERSORT-LM22, but it has learnt this from 12 microarrays of purified, resting and activated polyclonal γδ T cells expressing only TCRVγ9Vδ2, raising the question of whether it detects as efficiently the other (non-TCRVγ9Vδ2) subsets of γδ T lymphocytes. So far, the scarce publicly available microarrays for such purified γδ T cells do not allow to address this point. Nevertheless, the less performant deconvolution of transcriptomes for γδ cell bulks highly purified from fetal blood,¹³ although far better with LM7 than LM22, indicates that the current learning process enabled CIBERSORT-LM7 to reliably identify as γδ cells those cells sharing the gene signature of blood TCRVγ9Vδ2 γδ T cells. Future integration of other γδ T cell subsets purified from normal tissues or tumor samples and profiled on microarray or RNA Seq will enlarge the learning of CIBERSORT-LM7 to empower its identification capacities. However, since learning from more than one hundred CD4⁺ and CD8⁺ microarrays does not empower CIBERSORT to discriminate between the various TCR-defined subsets of αβ T cells, the learning from smaller sampling data sets might set limits for performing TCR-defined γδ T cell identification.

The overall output of CIBERSORT-LM7 matches that of CIBERSORT-LM22 but with a better reliability for the detection of γδ T cells, which were previously mistakenly identified as CD4⁺ (31%) or NK cells (23%), while 8% of CD8⁺ T cells were mistakenly classified as γδ T lymphocytes (Table S4). As a consequence, the association of γδ cells with a favorable outcome across a large number of cancer patients should be taken with caution.¹⁰ Our study reports TCRVγ9Vδ2⁺ γδ TILs levels that are fully consistent with the body of literature about various types of cancers (^31-34 and see below). For B-ALL, however, this result could reflect either the actual recruitment of γδ TILs or the bi-phenotypic profile of such leukemia, which frequently co-express rearranged Ig and γδ TCR genes.^35,36 In CML, bone marrow-infiltrating TCRVγ9Vδ2⁺ γδ T lymphocytes have not been reported, but the reported average of 8% CD3, 0% CD4⁺ and <4% CD8⁺ T cells in such tissues³⁷ is compatible with our results. In AML, the few studies which have quantified TCRVγ9Vδ2⁺ γδ TILs in bone marrow biopsies have found them to be at similar levels as in the blood, but these figures might differ across patients and FAB subtypes.^19,38 In solid tumors, independent studies have mentioned an abundance of γδ TILs in cancers of the pancreas,³⁴ breast,³¹ head and neck,³³ colorectal³⁹ and lung,⁴⁰ but there have been no reports for prostate, cervix and esophagus cancers. Moreover, beside of the pioneering CIBERSORT studies which sub-optimally detect γδ cells,^9,10 there are currently no published meta-analyses of TCRVγ9Vδ2⁺ γδ TIL across different human cancers.

The results presented here detected highly heterogeneous amounts of TCRVγ9Vδ2⁺ γδ TILs in cancers, and their relative abundance in the above solid cancers as well as in prostate and esophagus carcinoma. The high level of the TCRγ gene expression by prostatic epithelial cells⁴¹ might set pitfalls to identify γδ TILs in prostate cancer by conventional methods, but not to the deconvolution depicted here which takes in account all the LM7 genes for this aim. So far however, the abundance of TCRVγ9Vδ2γδ TILs in prostate cancer was not reported, so this study brings further rationale to harness these cells for immunotherapeutic purposes.^42-45 This work also reports abundant TCRVγ9Vδ2 γδ TILs in breast carcinoma, a finding in line with earlier studies demonstrating that these γδ cells are abundant but associated with immunosuppressive activity, tumor progression, lymph node metastases and correlate with poor outcome.³¹ Likewise, the stroma of ductal adenocarcinoma was reported to comprise abundant γδ TIL with unfavorable role,³⁴ supporting our result for this disease. Actually, a detrimental role is reportedly associated to γδ TILs in breast, colorectal and pancreas cancer,^31,34,39 raising the question of whether their abundance generally predicts a favorable or unfavorable outcome. The results presented here illustrate an association of TCRVγ9Vδ2⁺ γδ TILswith favorable clinical evolution in four cohorts of patients with different types of cancer (Fig. 4). This seemingly contrasts with an earlier study depicting promotion of colorectal cancer in patients with abundant Th17-producing γδ TILs.³⁹ Nevertheless, the clinical significance of γδ TILs abundance in cancer most probably depends on the TCR subsets of γδ cells,⁷ their maturation and cytokine profile⁴⁶ and the tumor contexts as demonstrated in AML,⁴⁷ breast,³¹ colorectal³⁹ and pancreatic³⁴ cancers.

The abundance of αβ and TCRVγ9Vδ2⁺ γδ TIL across the spectrum of cancers tested varied considerably. The ordering of cancers according to their abundance of αβ TILs presented a similar pattern as that seen when they are ordered by mutational load.^48,49 This was illustrated by children's cancers, which typically have few mutations and here presented the least αβ TILs whereas cancers associated with durable exposure to carcinogens had the greatest number of these cells (Fig. 3). This is in agreement with the finding that the diversity and amount of αβ TIL clonotypes increases with the amount of non-synonymous mutations within tumors.⁵⁰ Nevertheless, we detected a different distribution of TCRVγ9Vδ2⁺ γδ TILs, suggesting that their incentive for infiltrating tumors is not mutational as for αβ TILs. The lack of identification of transcriptomic correlates for the abundance of  TCRVγ9Vδ2⁺ γδ TILs across the cancer spectrum could be due to several diverse reasons: the screen for correlates involved only MSIgDB-predefined human gene sets not unknown ones, this human transcriptome-based strategy cannot detect associations with most viral⁵¹ or microbiotal contexts,^52,53 this strategy might be intrinsically unable to detect the non-transcriptional hallmarks of cancers such as specific somatic mutations or the accumulation of endogenous (onco) metabolites. Due to the metabolic and microbial origin of the PAg that activate most human blood TCRVγ9Vδ2⁺ γδ lymphocytes,² such hypotheses deserve further investigation. Future research exploiting the increasing body of RNASeq and metabolomics data sets from human cancers will certainly help in addressing these issues.

Materials and methods

Signature matrix

Public raw data of (n = 265) transcriptomes from purified cell types obtained using Affymetrix HGU133 Plus 2.0 microarrays were downloaded from https://www.ncbi.nlm.nih.gov/gds, normalized together and collapsed to HUGO gene symbols using chipset definition files available from the NCBI GEO. These are referred to as reference data (Table S2). As a filtration step, genes with enriched expression (mean log₂ (expression level) >7) in all the normal non-hematopoietic tissues from the Human Body Index (study GSE7307) were removed. We then defined seven leucocyte subsets (B cells, CD4⁺ T cells, CD8⁺ T cells, γδ T cells, Mo-Ma-DC cells and granulocytes) which were imported into the CIBERSORT software (https://cibersort.stanford.edu/)⁹ to build the LM7 matrix (Table S3). For validation, CIBERSORT-LM7 was applied to a variety of external normalized public data sets (Table S4) containing purified leucocyte cells subsets including non-TCRVγ9Vδ2 γδ T cells isolated from human cord blood⁵⁴ or mixed populations of known composition.

Abundance calculations

Abundances were calculated from the CIBERSORT results and the SES, as depicted in the text. The SES¹⁴ were computed by applying the open source software AutoCompare-SES (https://sites.google.com/site/fredsoftwares/products/autocompare_ses) with normalized settings. Then, the open source software DeepTIL (https://sites.google.com/site/fredsoftwares/products/deeptil) was used to automatically compute the abundance of the seven leucocyte subsets in the samples. DeepTIL was written using the cross-platform Julia programming language⁵⁵; its typical computation time is 3–5 sec per sample. A parallel version of DeepTIL is also available for computing data tables simultaneously. From CIBERSORT and SES calculations, DeepTIL computes the following steps: (1) selection of LM7 genes with a mean log₂ (expression level) >11.5 in each leucocyte type to constitute its defining gene set. In each sample, the SES of the seven leucocyte defining gene sets are computed against the testing alternative “greater than”.¹⁴ The SES of each type of leucocyte is also computed for each of the seven pure types of leucocyte to determine the seven corresponding correcting factors. (2) a linear combination of the seven leucocyte type SESs weighted by their respective correcting factors is computed for each sample. (3) for each sample, CIBERSORT fractions are finally multiplied by the linear combination of SES to produce the abundance of each leucocyte subset. The mean abundance of each leucocyte subset in ∼10,000 samples of human cancers is summarized (Table S6).

Statistical analysis

Assessments of leucocyte fractions from the specified transcriptomes were performed using CIBERSORT (https://cibersort.stanford.edu/)⁹ with 500 Monte Carlo iterations and the matrix indicated in the text. Comparisons of TIL (αβ and γδ) abundances across cancers (Fig. 3) were performed by Wilcoxon tests, assessing the abundance of TIL in the specified cancer type comparatively to the same variable across the whole set of hematological or solid cancers. The Benjamini–Hochberg (BH) correction⁵⁶ was applied for multiple comparisons. The functionally defined collection of gene sets was downloaded from MSIgDB (http://software.broadinstitute.org/gsea/msigdb),²² while for control experiments, we generated a collection of 500 random gene sets of the same size as MSIgDB, referred to as RGS. For the transcriptome samples specified in the text, the SES of each gene set was computed as depicted¹⁴ and their respective correlations (Pearson and Spearman) with the abundance of the specified cell types were automatically computed using the “Correlations by lines in list” software (https://sites.google.com/site/fredsoftwares/products/correlation-by-lines-in-list). For Kaplan–Meier plots, optimal cutoffs were determined with the survival R package⁵⁷ and Log-Rank p-values were corrected using the BH method.⁵⁶

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Acknowledgments

We acknowledge JJF's team members for their stimulating comments and Kelly Thornber from Scientific Scripts for English editing of the text.

Funding

This work was funded by the Institut National de la Santé et de la Recherche Médicale (INSERM), the Université Toulouse III: Paul Sabatier, the Center National de la Recherche Scientifique (CNRS), the Laboratoire d'Excellence Toulouse Cancer (TOUCAN) (contract ANR11-LABX) and the Program Hospitalo-Universitaire en Cancérologie CAPTOR (contract ANR11-PHUC0001).