Multiomics profiling reveals the benefits of gamma-delta (γδ) T lymphocytes for improving the tumor microenvironment, immunotherapy efficacy and prognosis in cervical cancer

Abstract

Background

As an unconventional subpopulation of T lymphocytes, γδ T cells can recognize antigens independently of major histocompatibility complex restrictions. Recent studies have indicated that γδ T cells play contrasting roles in tumor microenvironments—promoting tumor progression in some cancers (eg, gallbladder and leukemia) while suppressing it in others (eg, lung and gastric). γδ T cells are mainly enriched in peripheral mucosal tissues. As the cervix is a mucosa-rich tissue, the role of γδ T cells in cervical cancer warrants further investigation.

Methods

We employed a multiomics strategy that integrated abundant data from single-cell and bulk transcriptome sequencing, whole exome sequencing, genotyping array, immunohistochemistry, and MRI.

Results

Heterogeneity was observed in the level of γδ T-cell infiltration in cervical cancer tissues, mainly associated with the tumor somatic mutational landscape. Definitely, γδ T cells play a beneficial role in the prognosis of patients with cervical cancer. First, γδ T cells exert direct cytotoxic effects in the tumor microenvironment of cervical cancer through the dynamic evolution of cellular states at both poles. Second, higher levels of γδ T-cell infiltration also shape the microenvironment of immune activation with cancer-suppressive properties. We found that these intricate features can be observed by MRI-based radiomics models to non-invasively assess γδ T-cell proportions in tumor tissues in patients. Importantly, patients with high infiltration levels of γδ T cells may be more amenable to immunotherapies including immune checkpoint inhibitors and autologous tumor-infiltrating lymphocyte therapies, than to chemoradiotherapy.

Conclusions

γδ T cells play a beneficial role in antitumor immunity in cervical cancer. The abundance of γδ T cells in cervical cancerous tissue is associated with higher response rates to immunotherapy.

Correspondence to Professor Zhiyong Yuan; [email protected] ; Dr Jie Chen; [email protected]

WHAT IS ALREADY KNOWN ON THIS TOPIC

• γδ T cells play a tumor-promoting or tumor-suppressing role in a variety of cancers. The results of studies based on cell lines or animal models are divergent in their understanding of the role of γδ T cells in the microenvironment of cervical cancer.

WHAT THIS STUDY ADDS

• In patients with cervical cancer, γδ T cells play a role in suppressing tumor progression mainly through direct cytotoxicity and shaping the immune-activated microenvironment. Compared with chemoradiotherapy, patients with high levels of γδ T-cell infiltration may benefit more from immunotherapy, including immune checkpoint inhibitors and autologous tumor-infiltrating lymphocyte therapy.

HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY

• Upon further validation, assessment based on γδ T-cell status can stratify cervical cancer patients for individualized treatment strategy development, especially immunotherapy. Cellular immunotherapy based on γδ T cells demonstrates promising prospects.

Background

γδ T cells represent an unconventional and infrequent subset of T cells, accounting for approximately 1%–5% of T cells in the peripheral blood of humans.¹ They mainly reside in peripheral mucosal barriers. These cells play crucial roles in maintaining physiological homeostasis and orchestrating immune responses during pathological conditions.² γδ T cells demonstrate independence from the major histocompatibility complex (MHC) and do not necessitate MHC-mediated antigen presentation. Notably, they have been identified as indispensable components that play a dual immune role within the tumor immune microenvironment (TIME). For instance, γδ T cells exhibit potent cytotoxicity and produce interferon (IFN)-γ, thereby serving as effectors in protective immune responses. However, a subset of γδ T cells that produces interleukin-17 can exert tumor-promoting functions, thereby highlighting their potential deleterious role.^{3 4} The association of γδ T cells with a favorable prognosis has been demonstrated in non-small cell lung cancer and gastric cancers, while conversely, it has been linked to an unfavorable prognosis in gallbladder cancer and acute myeloid leukemia.²

The successful isolation of γδ T cells from the human cervix was accomplished three decades ago.⁵ Cervical cancer ranks second among gynecologic malignancies worldwide in terms of both incidence and mortality, with approximately 604,000 new cases and 342,000 deaths globally in 2020.⁶ Among these, approximately 90% of new cases and fatalities related to cervical cancer are concentrated in less developed regions.⁷ Despite the remarkable advancements in vaccination, early screening and treatments, recurrence is observed in approximately one-third of patients with invasive cervical cancer subsequent to initial treatment, typically occurring within a span of 3 years.⁸ Currently, the National Comprehensive Cancer Network Clinical Practice Guidelines in Oncology (NCCN Guidelines) recommends combination therapies involving immune checkpoint inhibitors (ICIs), such as pembrolizumab and nivolumab, as a treatment option for recurrent or metastatic cervical cancer.

As our comprehension of the role played by common cell subtypes in tumor immunity deepens, it becomes increasingly important to elucidate the immune mechanisms of rare immune cell types that may also exert significant influence, particularly their association with immunotherapy. However, the precise role of γδ T cells in cervical cancer development, treatment response, and prognostic impact remains elusive due to technical challenges arising from their limited abundance and the lack of conservation of the T-cell receptor (TCR) repertoire and function between humans and mice.⁹ Here, we leveraged plentiful data from single-cell and bulk transcriptome sequencing, mutational spectra from whole exome sequencing and genotyping array, and radiomics based on MRI to explore the role of γδ T cells in the tumor microenvironment (TME) of cervical cancer and explore their clinical significance through a multiomics analysis strategy.

Methods

The overall flow chart of this study is shown in figure 1.

Figure 1. Schematic design of the study. CNAs, copy-number alterations; FACS, fluorescence-activated cell sorting; GSE, Gene Expression Omnibus series; ICIs, immune checkpoint inhibitors; PBMCs, peripheral blood mononuclear cells; RNA-seq, RNA sequencing; scRNA-seq, single-cell RNA sequencing; SNVs, single nucleotide variants; TCGA, The Cancer Genome Atlas; TCIA, The Cancer Imaging Archive; TILs, tumor-infiltrating lymphocytes; T2WI, T2 weighted image; WXS, whole exome sequencing.

Subjects and data acquirement

In accordance with the principles stated in the Declaration of Helsinki, approval from the ethics committee and written informed consents from the participants were obtained. First, formalin-fixed paraffin-embedded (FFPE) tissue specimens and the survival status of 47 patients with cervical cancer who underwent surgery from March 15, 2013 to December 10, 2018 were retrieved retrospectively. Second, biopsies were performed for three patients with cervical cancer prospectively, including two primary lesions and one osseous metastasis lesion, between July 22, 2022 and October 24, 2022. Third, bulk transcriptome profiling, variations and clinical data were obtained from the CESC (cervical squamous cell carcinoma and endocervical adenocarcinoma) cohort of The Cancer Genome Atlas (TCGA) Research Network (release: July 27, 2022) through the National Cancer Institute Genomic Data Commons.^{10 11} Fourth, the single-cell RNA sequencing (scRNA-seq) data of peripheral blood γδ T cells and tumor-infiltrating lymphocytes (TILs) of cervical cancer were obtained from Gene Expression Omnibus Series GSE128223¹² and GSE190075¹³ respectively. Fifth, the MRI image data were obtained from The Cancer Imaging Archive (TCIA).^{14 15}

Immunohistochemistry

FFPE tissue sections were subjected to immunohistochemical staining using the HRP-polymer method. After dewaxing, hydration, heat induced antigen retrieval with 1X Tris-EDTA buffer (pH 9.0), endogenous peroxidase inactivation and non-specific site blocking, the samples were incubated with TCR γ/δ monoclonal antibody TCR1061 (clone 5A6.E9, Thermo Fisher Scientific, Waltham, Massachusetts, USA) at 4°C overnight and were rewarmed to room temperature. Next, the samples were incubated with horseradish peroxidase-conjugated secondary antibodies against mouse IgG at room temperature for 30 min. Then, a 3,3’-diaminobenzidine peroxidase substrate kit was applied, and cell nuclei were counterstained with hematoxylin. Following dehydration, transparency and encapsulation, five high-power fields (HPF; 200×) per slide were captured randomly on an Olympus BX53 microscope. Finally, the numbers of γδ TCR⁺ cells were quantified by pathologists blinded to the sample information.

Single-cell transcriptome sequencing and data analysis

Single-cell sequencing

scRNA-seq was performed on prospectively collected samples. The processes of single-cell preparation, barcoding, library construction and sequencing are described in the online supplemental material. The scRNA-seq data generated in this study were deposited in the Genome Sequence Archive¹⁶ of the National Genomics Data Center¹⁷ with reference number HRA005178.¹⁸

Preprocessing and quality control of scRNA-seq data

Reads alignment, basic filtering, barcode counting and unique molecular identifier counting of scRNA-seq raw data were performed using Cell Ranger to generate analysis-ready feature-barcode matrices. The basic downstream analysis was performed in Seurat (V.4.3)¹⁹ based on R (V.4.2.3). The routine quality control (QC) requirements were implemented for each scRNA-seq data set (online supplemental figure 1), and doublets were identified by DoubletFinder (V.2.03).²⁰ Linear and graph-based dimensionality reduction was performed, and cells were clustered using the Louvain algorithm as described in the online supplemental material.

Annotation of cell clusters

Cell clusters were annotated based on a previous study²¹ and databases,^{22 23} and γδ T cells were distinguished (expressing PTPRC, CD3D, CD3E, TRDC, TRGC1 and TRGC2 highly, expressing CD8A and CD8B very lowly, and not expressing CD4, FGFBP2 or KLRC1). Annotations for other cell clusters and the marker genes can be found in the online supplemental material.

Trajectory analysis and cell–cell communication analysis

Between-sample batch effects were removed before trajectory construction using Harmony²⁴ (V.0.1.1). Trajectory analysis was conducted using Slingshot²⁵ (V.2.8). Monocle²⁶ (V.2.28) was also used for trajectory exploration and visualization. CytoTRACE²⁷ (V.0.3.3) was used to assess the differentiation potential of cell states thereby determining the temporal order of trajectory. Cell–cell communication analysis was conducted using CellChat²⁸ (V.1.6.1). Cell–cell contact and secreted signaling were observed separately.

Estimation of cell infiltration from bulk transcriptome data

The CIBERSORTx²⁹ deconvolution algorithm was employed to quantitatively estimate γδ T-cell infiltration levels from bulk sequencing data of cervical cancer tumor tissue in absolute mode. The scRNA-seq data set GSE128223,¹² consisting of human TCR Vδ1 and TCR Vδ2 γδ T cells purified from peripheral blood samples of healthy adults, was used as the reference for generating the signature matrix. The detailed parameters for the deconvolution can be found in the online supplemental material. Ultimately, we quantified the absolute levels of γδ T-cell infiltration in bulk sequencing data by assigning a numerical score. The estimation of infiltration levels of other immune cells was based on the LM22 reference set³⁰ using a similar approach.

Differential expression and gene set analyses

Differential expression analyses were performed using the limma³¹ package (V.3.56.2) or a Wilcoxon rank-sum test. Gene set analyses were performed based on the AUCell package,³² gene set variation analysis (GSVA)³³ or gene set enrichment analysis (GSEA).³⁴ The clusterProfiler³⁵ package (V.4.8.2) was used for visualization. A more detailed description is attached in the online supplemental material.

Analysis of focal somatic variants

Masked somatic mutation data in mutation annotation format were obtained for the analysis of single nucleotide variants (SNVs). GISTIC³⁶ (V.2.0.23) was used for localization of focal somatic copy-number alterations (CNAs) based on masked copy number segments. The somatic SNVs and CNAs were comprehensively characterized using Maftools³⁷ (V.2.16).

In silico estimation of response to ICIs

The EaSIeR³⁸ package (V.1.6.3) was used to compute response scores, enabling the prediction of response to ICIs based on interpretable systems biomarkers.

Radiomics feature extraction, selection and model establishment

The T2-weighted MRI images of 53 patients with cervical cancer in the TCIA database, were subjected to radiomics analysis. Technical details of the preprocessing and feature extraction of the image data can be found in the online supplemental material. In brief, radiomic features were extracted using the ontology-guided radiomics analysis workflow,³⁹ which operates based on PyRadiomics⁴⁰ (V.3.7). The majority of the extracted features were in accordance with the guidelines set forth by the Image Biomarker Standardization Initiative.⁴¹ After Z-score normalization, radiomic features significantly correlated with γδ T-cell abundance (r>0.3, p<0.05) were initially screened based on Pearson’s correlation. After determining the optimal λ value through leave-one-out cross validation, feature selection was performed using least absolute shrinkage and selection operator (LASSO) regression. The resulting features were used to construct a linear regression model for calculating γδ T-cell scores, and the effectiveness of the model was assessed using the glmnet package⁴² (V.4.1–7).

Statistics

Kaplan-Meier survival analysis and corresponding statistical tests were performed using MedCalc Statistical Software (V.19.6.1). The grouping cut-off value for continuous variables was determined based on the median. According to the applicability of various statistical test methods for different types or distributions of data, Welch’s t-test, Wilcoxon rank-sum test, adjusted χ² test or Fisher’s exact test were performed in R (V.4.3.1) or MedCalc Statistical Software (V.19.6.1). The ggstatsplot⁴³ package (V.0.12) was also used, and the effect size, Hedge’s g, was computed.

Results

γδ T cells infiltrate human cervical carcinoma tissue

To investigate whether γδ T cells infiltrate human cervical carcinoma tissue, FFPE sections were stained for γδ TCR by immunohistochemistry. We observed that γδ TCR⁺ cells were mainly present in the tumor stroma (figure 2A), and the results showed heterogeneity in their abundance with a median count of 13 cells/HPF (200×) (figure 2B,C). Moreover, there was a tendency for the number of γδ T cells to vary in patients with different survival status (median 15 vs 8 cells/HPF, Wilcoxon test, p=0.1716) (figure 2D). The overall survival (OS) was significantly improved (HR=0.303 (95% CI: 0.097 to 0.947); log-rank test, p=0.040) by the infiltration of γδ T cells (figure 2E). Therefore, further verification and interpretation of the effect of γδ T cells is necessary.

Figure 2. The infiltration level of [gamma][delta]T cells in tumor tissue was associated with an overall favorable prognosis for patients with cervical cancer. (A) Example immunohistochemical staining plot of [gamma][delta] TCR + cells in cervical cancer tissues under 200×high-power field (HPF). (B) Example immunohistochemical staining plots with high or low abundance of [gamma][delta] TCR + cells in the 400×field of view. (C) Bar graph showing the frequency distribution of the number of [gamma][delta] TCR + cells per 200×HPF. (D) Box plots demonstrating the number of [gamma][delta] TCR + cells per 200×HPF in patients with different survival statuses. (E) Kaplan-Meier curves show overall survival (OS) for patients with different [gamma][delta] TCR + cell counts per 200×HPF. (F) t-distributed stochastic neighbor embedding (t-SNE) plots of the [gamma][delta] T-cells reference data set, colored by V[delta]1/V[delta]2 T-cell subtypes (left) or donors (right). (G) t-SNE-based feature plots show the expression of typical [alpha][beta]/[gamma][delta] T cells marker genes. (H) Kaplan-Meier curves show OS (left) and progression-free interval (right) for patients with high or low [gamma][delta] T-cell infiltration absolute scores in the The Cancer Genome Atlas-cervical squamous cell carcinoma and endocervical adenocarcinoma cohort. TCR, T-cell receptor.

High abundance of γδ T cells correlates with a favorable prognosis in cervical cancer

Assessing the condition of γδ T cells in a sufficient sample size of patients was a prerequisite for further investigating their effect. To measure the abundance of γδ T cells in cervical carcinoma tissue in a high-throughput way, scRNA-seq data of 8,140 γδ T cells (including 4,243 Vδ1 T cells and 3,897 Vδ2 T cells) purified from adults’ peripheral blood by fluorescence-activated cell sorting were recruited as a reference set (figure 2F). The cell purity and data quality of this reference set were satisfactory (figure 2G, online supplemental figure 1A). Based on this, digital cytometry of γδ T cells was performed based on RNA-seq data of each primary tumor sample from the TCGA-CESC cohort (n=304, online supplemental figure 2) using the CIBERSORTx deconvolution algorithm. Due to the extensive similarity between γδ T cells and CD8⁺ αβ T cells at the transcriptome level, to avoid confusion during deconvolution, we verified that the infiltration levels of the two subtypes of T cells were not correlated (Pearson’s correlation coefficient, r=−0.006, p=0.92) as expected, consistent with a previous study⁴⁴ (online supplemental figure 3). Next, patients with follow-up information (n=291, online supplemental figure 2) were classified into two groups according to the median of absolute score of γδ T-cell abundance. The patients were thus divided into a low γδ T-cell infiltration level group (n=146) and a high infiltration level group (n=145). There was no significant difference in the clinical characteristics of the two groups, including age, race, the International Federation of Gynecology and Obstetrics (FIGO) stage, histological type, prior or synchronous malignant tumor, smoking, etc (online supplemental table 1). Survival analysis showed that patients in the high infiltration group exhibited significantly improved OS (HR=0.605 (95% CI: 0.378 to 0.970); log-rank p=0.037) and progression-free interval (PFI) (HR=0.614 (95% CI: 0.384 to 0.981); log-rank p=0.041) compared with those in the low infiltration group (figure 2H). The 3-year OS rate (79.7% vs 64.4%) and PFI rate (75.4% vs 62.3%) were significantly higher in the high infiltration group than in the low infiltration group. Thus, we verified that γδ T cells have a heterogeneous infiltration level in cervical cancer tissues and that their infiltration level was positively correlated with the prognosis of patients.

Somatic mutation landscapes associated with the abundance of γδ T cells

In the above analysis, we found that clinical characteristics were not determinants of the level of γδ T-cell infiltration (online supplemental table 1). Next, to investigate whether the somatic mutations carried by tumor cells affect their infiltration levels, all patients in the TCGA-CESC cohort with RNA-seq data and SNV or CNA information (n=300, online supplemental figure 2) were redistributed into two γδ T-cell abundance groups according to the median of γδ T-cell score. The TNN gene exhibited the highest frequency of SNVs in the higher infiltration group (34% in the higher infiltration group; 24% in the lower infiltration group), whereas PIK3CA displayed the highest SNV frequency in the lower infiltration group (24% in the higher infiltration group; 31% in the lower infiltration group) (online supplemental figure 4A,B). Genes such as GTF3C1, AFF3, MYH9, FAT3, and CELSR1 showed significantly higher (p<0.01) somatic variation frequencies in the high γδ T-cell abundance group than in the other group (figure 3A, online supplemental table 2). The co-mutation patterns in the high γδ T-cell abundance group were more plentiful than those in the other group (figure 3B). The genes DPP6, LRGUK, and TMEM209 exhibited higher CNA frequencies in the high infiltration level group, whereas SBNO2 and GPX4 showed higher CNA frequencies in the low infiltration level group (figure 3C). Generally, the CNA landscape was similar between the two groups, but amplifications at 8q24.21 and 9p24.1 were more common in the high γδ T-cell abundance group, while deletion at 3p14.1 was more common in the other group (figure 3D, online supplemental figure 4C,D). The majority of low-frequency somatic variations in the mentioned genes did not show a significant prognostic impact. However, the variations of GPX4 and SBNO2, which were more frequently altered in patients with low levels of γδ T-cell infiltration, showed significantly adverse prognostic effects (online supplemental figure 5A,B). These findings highlight notable discrepancies in the somatic variation profiles of the two groups, suggesting a nuanced relationship between somatic variations and γδ T-cell infiltration levels.

Figure 3. The focal somatic mutations and tumor microenvironment characteristics of cervical carcinoma tissue were associated with the infiltration level of [gamma][delta] T cells. (A) Forest plot demonstrating focal somatic single nucleotide variant (SNV) and copy-number alteration (CNA) frequencies that are significantly different between patients in the [gamma][delta] T-cell high and low infiltration level groups. OR>1, higher SNV frequency for patients in the [gamma][delta] T-cell high infiltration level group. OR<1, higher SNV frequency in low infiltration level group. OR=1, the frequency of SNVs was consistent in both groups. (B) Landscape of gene co-mutations in patients with low (left) and high (right) abundance of [gamma][delta] T cells. (C) Bar diagram showing the frequency of SNVs and CNAs of the genes featured in panel A across the two groups of cases. (D) Visualization of CNAs in patients stratified by low (left) or high (right) infiltration levels of [gamma][delta] T cells in chromosomal context. Significantly amplified or deleted genomic regions with FDR<0.05 are indicated in red or blue. The height of the peaks represents the G-score computed by GISTIC2. (E) Differentially expressed genes in tumor tissues of patients with high [gamma][delta] T-cell infiltration. The representative upregulated or downregulated genes are labeled. (F) Sample-wise gene set enrichment scores of the top 10 upregulated and downregulated biological process or molecular function gene sets of Gene Ontology. Both samples and gene sets are ordered by complete clustering. The annotated color bar at the top demonstrates the [gamma][delta] T-cell infiltration grouping for each sample. (G) Top 10 upregulated and downregulated pathways in patients with a high abundance of [gamma][delta] T cells identified by gene set enrichment analysis (GSEA) based on hallmark gene sets in the MSigDB database. (H) Representative immune-related or malignancy-related pathways significantly enriched in patients with a high abundance of [gamma][delta] T cells identified by GSEA based on the Kyoto Encyclopedia of Genes and Genomes database. FDR, false discovery rate. Additional enriched pathways can be found in the online supplemental figure 6 .

γδ T cells play a critical antitumor role in the cervical cancer microenvironment

We have demonstrated that γδ T cells provide a prognostic benefit to patients, and we next explored how this benefit arises. A correlation analysis between the infiltration level of γδ T cells and other immune cells was performed (online supplemental figure 3). The results showed that the infiltration level of γδ T cells negatively correlated with regulatory T cells, but was positively associated with memory CD4⁺ T cells in the resting state. Moreover, the increase in γδ T cells was associated with a decrease in memory B cells and plasma cells. Interestingly, the infiltration level of γδ T cells was significantly positively associated with M1-type macrophages. In addition, as previously mentioned, the abundance of γδ T cells was not associated with CD8⁺ αβ T cells and natural killer (NK) cells, which is consistent with a previous study.⁴⁴ These results implied that infiltration of γδ T cells was associated with a decrease in immunosuppressive cells and an increase in cells associated with antitumor immunity.

Next, all patients in the TCGA-CESC cohort who underwent RNA-seq (n=304, online supplemental figure 2) were redivided into two groups, according to the median of γδ T-cell score. Differential analysis of gene expression indicated extensive differences between the two groups (figure 3E, online supplemental table 3). Meanwhile, the patients from the two groups could be distinctly clustered based on GSVA, indicating a significant difference in TME between the two groups (figure 3F). Interestingly, GSEA revealed that in the TME of patients from the γδ T-cell high-abundance group, antitumor immune-related pathways were generally upregulated, while pathways associated with tumor malignancy were downregulated (online supplemental table 4 and 5). For instance, pathways involving IFN-γ, IFN-α, tumor necrosis factor (TNF)-α, inflammatory response, p53 tumor suppressor protein, cell apoptosis, cytokine–cytokine receptor interaction, and JAK-STAT signaling were significantly upregulated; on the other hand, pathways related to oxidative phosphorylation, hedgehog signaling, angiogenesis, glycolysis, KRAS signaling, and motility proteins were significantly downregulated. These results suggested that high infiltration levels of γδ T cells correlated with activation of antitumor immunity and suppression of tumor development (figure 3G,H, online supplemental figure 6).

Characteristic of γδ T cells within the TIME of cervical cancer

To characterize γδ T cells in the real immune microenvironment of cervical cancer at the single-cell level in more detail, we performed scRNA-seq on tumor tissues obtained from three patients with cervical cancer. A total of 16,667 cells that passed QC (online supplemental figure 1B) were subjected to dimensionality reduction (online supplemental figure 7A–D), clustering (online supplemental figure 7E,F), and annotation into 22 distinct cellular subtypes based on their respective signature genes (figure 4A). Subsequently, each subtype of CD8⁺ αβ T cells and γδ T cells was further classified into multiple clusters based on the expression levels of their highly specific marker genes (figure 4B). As expected, γδ T cells exhibited a transcriptome landscape comparable to that of CD8⁺ αβ T cells (particularly effector memory CD8⁺ T cells and tissue-resident memory CD8⁺ T cells), as well as NK cells, demonstrating similar proximity to them in both the t-distributed stochastic neighbor embedding and Uniform Manifold Approximation and Projection plots. Nonetheless, the marker genes employed for target cell selection exhibited robust discriminatory capabilities (online supplemental figure 8A). γδ T cells were predominantly present in primary foci, while they were hardly observed in bone metastases (figure 4A, online supplemental figure 7C,D). Gene differential expression analysis was performed on γδ T cells versus all other cells, and the 30 genes with the highest cell-averaged differential fold change were used as the signature of γδ T cells in the cervical cancer immune microenvironment (online supplemental table 6). Based on this signature, the abundance of γδ T cells in cervical cancerous tissue from the TCGA cohort (n=306), and normal cervix tissue from the Genotype-Tissue Expression cohort and the TCGA cohort (n=13) was compared. The abundance of γδ T cells in cervical cancerous tissues was significantly higher than that in normal cervix (online supplemental figure 8B). This may indicate that γδ T cells are recruited during tumorigenesis. Thus, γδ T cells were identified from scRNA-seq data of cervical carcinoma tissues, and studies on their function in the TME can follow.

Figure 4. Chasing down [gamma][delta] T cells from single-cell transcriptome data of cervical cancer tumor tissues. (A) Uniform Manifold Approximation and Projection (UMAP)-based (left) or t-SNE-based (right) dimensionality reduction map colored by cell cluster annotations. [gamma][delta] T cells are circled through the black curve for emphasis. (B) Dot plots demonstrating the expression levels of marker genes for each annotated cell cluster. Identities of cell clusters on the Y-axis are ordered by hierarchical clusters based on features shown. The size of the dot represents the proportion of cells expressing the gene in the corresponding cell cluster. The color of the dot represents the average expression level of the gene in the cell cluster. The annotations in the top gray box represent the interpretation of the corresponding marker genes on the X-axis. cDC, conventional dendritic cell; Endo/endothelial, endothelial cells; Epi/epithelial, epithelial cells; ICP, immune checkpoint; Macro, macrophage; Mast, mast cell; mFB, matrix fibroblasts; Mono, monocyte; Mye, myeloid cells; Neutro, neutrophil; NK, natural killer cell; pDC, plasmacytoid dendritic cell; Plasma, plasma cells; Prolif, proliferating markers; Tem, effector memory T cells; Texh, exhausted T cells; Treg, regulatory T cells; Trm, tissue-resident memory T cells; t-SNE, t-distributed stochastic neighbor embedding; vFB, vascular fibroblasts.

γδ T cells exert direct cytotoxic effects in the cervical cancer TME through the evolution of dichotomous states

In the in vivo immune microenvironment of cervical cancer, γδ T cells were found to exhibit two states. They can be distinguished based on the differential expression of marker genes, such as GNLY, resulting in their respective annotations as γδ T cells GNLY⁺ and γδ T cells GNLY^– (figure 4B). On removal of the between-sample batch effect, these two states could be clearly distinguished on the dimensional reduction plot (figure 5A). Cells in the GNLY⁺ state exhibited high expression levels of cytotoxicity-associated genes including GNLY, NKG7, GZMA, GZMB, GZMH, GZMK, and PRF1. Additionally, they also demonstrated elevated expression of immune signaling-related genes such as CCL4 and TYROBP. In contrast, cells in the GNLY⁻ state highly expressed TCR-related genes such as TRDC, TRGC1, TRGC2, etc (figure 5B, online supplemental table 7). Cells in the GNLY⁻ state had a higher differentiation potential (p=0.030) than those in the GNLY⁺ state (figure 5C), implying that the direction of evolution was from the GNLY⁻ state toward the GNLY⁺ state (figure 5D). Along this trajectory, γδ T cells undergo dynamic changes in their cellular state (figure 5E). In addition to the immune-related genes that were highly expressed in each of the two different states of cells described above, the expression of genes such as VCAM1, PTPN22, CLEC2D, ZNF80, and S100A4 was gradually downregulated with evolutionary progression, whereas the expression of genes such as AREG, KLRB1, CD247, CCL4L2, and TNFRSF18 was gradually upregulated, and they may be evolutionary driver genes (figure 5F, online supplemental table 8). Therefore, it is evident that the evolution of the state of γδ T cells present in the TME of cervical cancer and thus may exert a direct antitumor effect by increasing the expression levels of cytotoxicity-related genes.

Figure 5. Dynamic transformation of [gamma][delta] T cells and their communication with other cells in the tumor microenvironment of cervical cancer. (A) UMAP plot showing [gamma][delta] T cells in the two states after removal of the between-sample batch effect. (B) Volcano plot demonstrating differences in the gene expression of [gamma][delta] T cells in the two states. (C) The differentiation potential scores of [gamma][delta] T cells in GNLY - status were significantly higher than those in GNLY + status. (D) Evolutionary trajectories of [gamma][delta] T cells depicted on a UMAP plot. Cells are stained in pseudotemporal order. (E) Density plots of [gamma][delta] T cells in two states along the evolutionary trajectory in panel D. (F) Heatmap of smooth expression values for the 50 genes with the largest expression changes in [gamma][delta] T cells along the evolutionary trajectory of panel D. Rows (ie, genes) were clustered using the Ward D method. Expression values were transformed using log normalization and variance stabilization prior to rendering. (G) Dot plot showing cell-cell contacts from cDCs and macrophages to CD8 + [alpha][beta] T cells and [gamma][delta] T cells. Colors represent the probability of contact. The sizes of the dots represent p values. (H) Dot plot showing secreted signaling communications from [gamma][delta] T cells to CD8 + [alpha][beta] T cells. (I) Circle plot showing the secreted signaling communications between cell subpopulations. The thickness of the lines represents the weights (ie, strength) of the interactions. Interactions from a certain cell subpopulation are depicted by lines of the same color. cDC, conventional dendritic cell; mFB, matrix fibroblasts; NK, natural killer cell; pDC, plasmacytoid dendritic cell; Plasma, plasma cells; Prolif, proliferating markers; Tem, effector memory T cells; Texh, exhausted T cells; Treg, regulatory T cells; Trm, tissue-resident memory T cells; UMAP, Uniform Manifold Approximation and Projection; vFB, vascular fibroblasts.

Cell–cell communication between γδ T cells and other components within the TME of cervical cancer

As a specific T-cell subtype, in addition to direct tumor-killing effects, we hypothesized that γδ T cells may also exert antitumor effects through their mutual regulation with other cell components. To explore the regulatory relationship, cell–cell contacts and secreted signaling communications networks within the TME of cervical cancer were established (online supplemental figure 9A–D). Cell–cell contacts associated with human leukocyte antigen (HLA) class I molecules (HLA-A, HLA-B, HLA-C, HLA-F) from professional antigen-presenting cells (eg, conventional dendritic cells (cDCs) and macrophages) to effector-memory, tissue-resident-memory and exhausted CD8⁺ αβ T cells were observed. Notably, communication through the classical HLA class I pathway or costimulatory factors, such as CD6 and CD80, from cDCs or macrophages to γδ T cells was not observed. However, it was noted that γδ T cells received HLA-E signaling via KLRC1/2 and NKG2A/C (figure 5G). For secreted signaling communications, γδ T cells showed extensive interactions with CD8⁺ αβ T cells, including communications through the CCL5-CCR1/5 ligand receptor pair (figure 5H,I, online supplemental figure 8C). However, the regulation of other immune components such as myeloid cells by γδ T cells was not clearly observed. In summary, the antitumor effect exerted by γδ T cells in the TME of cervical cancer was possibly dominated by direct cytotoxicity while possibly exerting a regulatory effect on the function of CD8⁺ αβ T cells.

Survival advantage from high γδ T-cell abundance was attenuated in patients receiving chemoradiotherapy

Based on the above analyses, the favorable prognostic effect of γδ T cells on cervical cancer and their antitumor function have been largely clarified, but whether these effects differ for patients receiving different therapies is not clear. To investigate the impact of radiotherapy and chemotherapy on the survival advantage conferred by γδ T cells, we conducted survival analyses on three subgroups of patients: (1) those who received neither radiotherapy nor chemotherapy (n=69); (2) those who received radiotherapy alone (n=35); and (3) those who received both radiotherapy and chemotherapy (n=152). Three subgroups of patients were categorized into high-scoring and low-scoring groups by the same threshold according to the median for the whole group of patients with survival data. The results demonstrated a significant correlation between γδ T-cell infiltration levels and OS in patients who did not receive radiotherapy or chemotherapy (HR=0.167 (95% CI: 0.047 to 0.590); log-rank p=0.006) (figure 6A, left). Then, there was a suggestive correlation with OS in patients who received radiotherapy alone (HR=0.216 (95% CI: 0.043 to 1.087); log-rank p=0.063) (figure 6A, center). However, no obvious association was observed between γδ T-cell infiltration levels and OS in patients who received both radiotherapy and chemotherapy (HR=0.880 (95% CI: 0.493 to 1.574); log-rank p=0.667) (figure 6A, right). The changes in survival benefit implied that the survival advantage conferred by γδ T cells may be attenuated in patients receiving chemoradiotherapy. Therefore, we believe that it makes sense to develop a personalized treatment strategy for patients with a high abundance of γδ T cells.

Figure 6. Patients with high levels of [gamma][delta] T-cell infiltration may present as more suitable candidates for immune checkpoint inhibitor therapy and autologous tumor-infiltrating lymphocyte (TIL) treatment in comparison to chemoradiotherapy. (A) Kaplan-Meier curves show overall survival for patients receiving neither radiotherapy nor chemotherapy (left), radiotherapy alone (center), or both radiotherapy and chemotherapy (right). (B) Scatter plot demonstrating a statistically significant positive correlation between the expression level of the top 30 genes signature of [gamma][delta] T cells in the tumor microenvironment of cervical cancer and the expression level of PD-L1 (left) and CTLA-4 (right). (C) Immunotherapy response scores predicted by the EaSIeR algorithm. Welch’s t-test was performed, and the effect size Hedge’s g was computed. (D) t-SNE plots for single-cell RNA sequencing data of TILs in the infused product, colored by patients and their condition responding or not responding to treatment. Patient numbers are consistent with the original study. (E) Feature plots show the distribution pattern of the AUCell score for the [gamma][delta] TCR constant region. (F) AUCell score for the [gamma][delta] TCR constant region in TILs from the infused product for patients responding or not responding to the treatment. TCR, T-cell receptor; TPM, transcripts per million; t-SNE, t-distributed stochastic neighbor embedding.

Patients with a high abundance of γδ T cells may be more amenable to ICI treatment

The addition of immunotherapy, represented by ICIs, has been a breakthrough in the comprehensive treatment of cervical cancer in recent years. We attempted to predict the therapeutic response to ICIs in two groups of patients with high or low abundance of γδ T cells. First, we observed that the expression levels of PD-L1 (Wilcoxon test, p=6.1×10⁻⁴), CTLA-4 (p=7.6×10⁻⁴), and IFN-γ (p=0.037) were significantly higher in patients with a high abundance of γδ T cells (online supplemental figure 10A). The expression of the aforementioned signature of γδ T cells was also correlated with the expression of PD-L1 (r=0.52, p=3.3×10⁻²²), CTLA-4 (r=0.75, p=2.2×10⁻⁵⁷), and IFN-γ (r=0.77, p=5.9×10⁻⁶²) based on Spearman’s correlation (figure 6B, online supplemental figure 10B). Next, the ICI treatment response scores were predicted for each patient by the EaSIeR algorithm. The ICI treatment response scores of patients with high infiltration levels of γδ T cells were significantly higher (p=2.65×10⁻³) (figure 6C). Therefore, the above results suggested that patients with cervical cancer with higher levels of γδ T-cell infiltration may be more amenable to treatment strategies with ICIs.

Infused γδ T-cell proportion correlated with treatment response in autologous TIL immunotherapy

In addition to ICIs, autologous TIL infusion in the clinical trial stage was also an individualized immunotherapy strategy available for patients with cervical cancer. To analyze the relationship between the response to autologous TIL infusion therapy and the γδ T-cell percentage in the infusion products, we obtained in vitro expanded TIL scRNA-seq data from patients with cervical cancer receiving this treatment in a clinical trial. A total of 25,689 TILs passed QC (online supplemental figure 1C). These TILs were derived from eight patients, including five patients who responded well to treatment and three patients who did not respond to treatment (figure 6D). The expression levels of γδ TCR constant region markers, CD4⁺ αβ T cells markers and CD8⁺ αβ T cells markers in TILs were evaluated by the AUCell algorithm. The expression levels of the three types of markers in TILs from each patient are visually shown in feature plots (figure 6E, online supplemental figure 10C). There was a significant difference in the proportion of cells with high expression of γδ TCR in the TILs of five responders and three non-responders, suggesting that the treatment response was associated with the proportion of infused γδ T cells (figure 6F, online supplemental figure 10D). In contrast, the proportions of CD4⁺ αβ T cells and CD8⁺ αβ T cells were not distinctly different between the two groups of patients (online supplemental figure 10E,F). The above results implied that γδ T cells play an important role in the efficacy of both immunotherapy strategies including ICI and autologous TIL infusion therapy.

γδ T-cell abundance in cervical carcinoma tissue could be measured by an MRI-based radiomic model

To assess the impact of γδ T cells on the TME from a macroscopic perspective and to explore non-invasive strategies for assessing their infiltration levels, radiomics analysis was performed on T2-weighted MRI sequences from patients in the TCIA database, whose transcriptome data were previously described (n=53) (figure 7A). Volume of interest was manually segmented by an experienced radiation oncologist and subsequently reviewed by a senior physician (figure 7B). Out of 824 radiomics features, 12 features were found to exhibit significant correlations with the γδ T-cell abundance of patients in the training set measured by the transcriptome based on Pearson’s correlation (r>0.3, p<0.05). Through LASSO regression for dimension reduction (figure 7C,D), a final set of eight radiomics features was used to establish a linear regression model for assessing γδ T-cell infiltration levels. The formula for measuring γδ T-cell abundance based on this model is attached in the online supplemental material. The multiple R² value of the model was 0.625, the adjusted R² value was 0.528, and the p value was 6.067×10⁻⁵ (online supplemental table 9). The γδ T-cell abundance predicted by the radiomics model showed significantly positive correlations with its abundance measured by the transcriptome in both the training set (r=0.79, p=1.3×10^–9) and the test set (r=0.56, p=0.048) (figure 7E). Survival analyses showed a tendency for radiomics-predicted γδ T-cell abundance to stratify the OS of patients with cervical cancer (HR=0.531 (95% CI: 0.218 to 1.289); log-rank p=0.162) (figure 7F). Thus, we not only found that the macroscopic effects of γδ T cells on the TME can be observed by MRI, but also obtained a non-invasive assessment method with potential value for predicting the infiltration level of γδ T cells in cervical carcinoma tissues.

Figure 7. MRI-based radiomics analysis established a linear regression model predicting [gamma][delta] T-cell abundance in cervical carcinoma tissue. (A) Schematic diagram of the radiomics analysis process. (B) An example of segmentation for gross tumor volume on T2 MRI. (C) Parameter tuning plot in least absolute shrinkage and selection operator (LASSO) regression. The topmost values represent the number of radiomics features screened by Pearson’s correlation that are incorporated into the LASSO regression model. [lamda], the weight of the L1 norm. The two vertical dashed lines represent [lamda] min and [lamda] 1se . [lamda] min , the value of [lamda] that results in the smallest mean value of the squared error of the regression model. [lamda] 1se , the value of [lamda] that results in the simplest model within an SE of [lamda] min . (D) Distribution of coefficients for variables in the LASSO regression. Each curve represents a radiomics feature filtered by Pearson’s correlation. The y-axis is their corresponding coefficient in the LASSO regression. The topmost values carry the same meaning as depicted in panel C. (E) Pearson’s correlations between z-score normalized [gamma][delta] T-cell abundance measured by the transcriptome and the fitted value (ie, the predicted [gamma][delta] T-cell abundance) by the linear regression radiomics model in the training set and test set. The gray area represents the range of 95% CIs. (F) Kaplan-Meier curves show overall survival for patients with different fitted values (ie, the predicted [gamma][delta] T-cell abundance) by the linear regression radiomics model. ROI, region of interest; T2WI, T2 weighted image.

Discussion

To the best of our knowledge, this study represents a pioneering exploration of the favorable role played by human γδ T cells in the TIME of cervical cancer at single-cell resolution through a multiomics approach and is the first to elucidate the impact of γδ T cells on the efficacy of chemoradiotherapy, ICIs and autologous TIL treatment in cervical cancer. Furthermore, based on the systematic analysis of abundant data from genome, transcriptome, protein phenotype, and MRI-based radiomics, this study comprehensively investigated the factors influencing the infiltration level of γδ T cells, their direct antitumor function, and their interaction with TME components in cervical carcinoma tissues.

The views of previous studies on the role of γδ T cells in the prognosis of cervical cancer are controversial. On the one hand, some researchers believe that γδ T cells can exert a killing effect on cervical cancer cells. For instance, a study based on cervical cancer cell lines, including HeLa, SiHa and CaSki, found that γδ T cells induced apoptosis in tumor cells after treatment with pamidronate.⁴⁵ However, γδ T cells have also been found to have a role in promoting tumor development. In a study based on a transgenic mouse model of human papillomavirus type 16 oncoprotein-induced carcinogenesis, local epithelial-associated γδ T cells were thought to be involved in the promotion of angiogenesis and cancer development.⁴⁶ Experiments with mixed conclusions in vitro or in animal models created confusion, and unraveling the role played by γδ T cells in patients with actual cervical cancer was necessary. In this study, the abundance of γδ T cells in cervical carcinoma tissues was quantified by cell counting based on immunohistochemistry and digital cytometry based on transcriptomic data, and the resulting prognostic analyses showed that γδ T cells play a favorable role overall.

Despite the lack of previous studies, what factors influence the level of γδ T-cell infiltration and how γδ T cells play their beneficial role in the actual cervical cancer TME are also questions worth exploring. We did not observe the effect of clinical factors such as age, race, pathology pattern and FIGO stage on γδ T-cell infiltration, but found that patients harboring certain somatic mutations were more likely to have high γδ T-cell abundance and that they possessed more plentiful co-mutational events. This may imply that further molecular typing of cervical cancer is relevant for patient stratification. As an unconventional T-cell subtype involved in intrinsic immunity, γδ T cells are thought to potentially kill tumor cells through direct cytotoxicity or through the indirect modulation of other immune components. The results of the present study showed that during tumorigenesis, γδ T cells in the TME increased and the expression levels of granulysin, granule protein, granzyme, and perforin were elevated through the evolution of the cellular state thereby exerting a direct cytotoxic antitumor effect independent of HLA molecules. Furthermore, γδ T cells also contribute to the immune microenvironment of antitumor immune activation. The results of the present study suggest that they can induce a decrease in regulatory T cells, a decrease in B-cell function, and a shift of macrophages to the M1 subtype in the TME of cervical cancer. In the TME rich in γδ T cells, immune-related pathways such as IFN-γ, IFN-α, TNF-α, cytokine−cytokine receptor interaction and JAK-STAT signaling were significantly activated. This may lead to the activation of antitumor-related pathways such as p53 tumor suppressor protein, cell apoptosis, etc, while oxidative phosphorylation, hedgehog signaling, angiogenesis, glycolysis, KRAS signaling, motility proteins and other pathways related to the malignant biological behavior of tumors were suppressed, resulting in the inhibition of tumor development.

Depending on the stage and pathology of the disease, the comprehensive treatment options for cervical cancer recommended by the NCCN Guidelines include cone excision, trachelectomy with lymphadenectomy, external beam radiotherapy and brachytherapy with concurrent platinum-containing chemotherapy, etc. Systemic therapy options include platinum plus paclitaxel (etoposide for small cell neuroendocrine carcinoma), pembrolizumab/nivolumab and bevacizumab. Therefore, chemoradiotherapy is now a cornerstone in the management of many patients. However, our results imply that the prognostic-improvement effect of γδ T cells may have a more pronounced advantage in patients who have received neither chemotherapy nor radiotherapy. The findings of this study also suggest that patients with a high abundance of γδ T cells may be more amenable to immunotherapy (ICIs or autologous TIL therapy). Therefore, the infiltration levels of γδ T cells may serve as a potential biomarker for evaluating the immune status of patients when formulating personalized treatment strategies.

Radiomics provides a way to non-invasively assess the biological characteristics of tumors. Previous studies have developed numerous CT-based⁴⁷ or MRI-based⁴⁸ radiomic biomarkers of tumor-infiltrating CD8⁺ T cells, although the radiomic signatures they reported were heterogeneous with limited reproducibility.⁴⁹ Here, we sought to establish the first radiomic biomarker for tumor-infiltrating γδ T cells, thus demonstrating its broad impact on the TIME. Without an external validation set, we do not have sufficient confidence to ensure that the obtained radiomics model robustly predicts the level of γδ T-cell infiltration. Under the condition of a limited number of patients and data quality, we used a conservative linear correlation-based and L1 norm-based algorithm for feature selection to reduce the risk of overfitting. In addition, this section aimed to explore whether the effect of γδ T cells on the TME could have a macroscopic impact to the extent that it could be captured by imaging tools such as MRI, and establishing a predictive model was not a primary aim of the study. The findings suggest that the selected MRI radiomic features can serve as indicators of γδ T-cell infiltration levels, thereby implying the significant influence of γδ T cells on the tumor and its microenvironment. This attempt is risky, but we expect it to lay the groundwork for the next step of in-depth research.

It should be noted that the present study was an in silico exploratory study based on high-throughput multiomics data and had the following limitations. First, the small sample size of cervical cancer tumor tissue scRNA-seq might affect the generalizability of the results. To reduce the inaccuracy in estimating the level of infiltration brought about by interpatient heterogeneity, we adopted an integrated strategy of using high-quality γδ T-cell scRNA-seq data purified from the peripheral blood of adult volunteers¹² as a reference set for quantity estimation, and then using the scRNA-seq data for cervical cancer tissues to explore their status and function in the cervical cancer TIME. Second, the application of deconvolutional algorithms is often considered to require the fulfillment of the consistency assumption, which poses a challenge for the estimation of rare cells. The good news is that a recent benchmark study⁵⁰ of RNA-seq deconvolutional analysis in a dynamic test environment showed that under the test conditions of large unknown biological content spike-ins, on an absolute measurement scale, the CIBERSORTx²⁹ algorithm obtains fine correlation with the ground truth in both the “orthog” weight matrix (r=0.93) and the “real” weight matrix (r=0.63). Furthermore, the immunohistochemistry-based assessment of γδ T-cell counts yielded similar results demonstrating prognostic benefit. Thus, based on such a comprehensive process, we believe the assessment of the abundance of γδ T cells and the functions they perform is reliable, although further investigations are needed to validate it.

On the basis of the immune activation and antitumor effects of γδ T cells on the TME of cervical cancer presented in this study, it is promising to delve deeper into the mechanisms by which γδ T cells play a unique role in the TME of humans. However, the function and TCR repertoire of γδ T cells exhibit limited conservation between humans and mice, and their abundance is scarce.⁹ Thus, we anticipate challenges in explaining the molecular mechanisms of our findings. Larger sample sizes of high-quality, single-cell level transcriptomic, immune repertoire, and proteomic data based on human cervical cancer samples may be able to play a pivotal role in addressing this issue.

Conclusion

γδ T cells, as a crucial component of the immune microenvironment in cervical cancer tumors, exert a pivotal influence. Patients exhibiting high levels of γδ T-cell infiltration demonstrate significantly improved prognosis and are more amenable to receiving immune-based therapies such as ICI and TIL treatments. Collectively, γδ T cells contribute favorably to antitumor immunity in cervical cancer.

The results published here are in part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga. We are grateful to the research groups making their data publicly available. Figure 1 was created with BioRender.com.

Data availability statement

Data are available in a public, open access repository. The scRNA-seq data generated in this study are available in the Genome Sequence Archive of the National Genomics Data Center at https://ngdc.cncb.ac.cn/?lang=en, reference number HRA005178. The data that support the findings of this study are as follows. The scRNA-seq data are available in NCBI-GEO at https://www.ncbi.nlm.nih.gov/geo/, reference numbers GSE128223 and GSE190075. The transcriptome profiling, variation and clinical data are available from the TCGA Research Network at https://portal.gdc.cancer.gov/projects/TCGA-CESC. The image data are available from the Cancer Imaging Archive (TCIA) at https://doi.org/10.7937/K9/TCIA.2016.SQ4M8YP4 . Additional data related to this paper may be requested from the authors.

Ethics statements

Patient consent for publication

Not applicable.

Ethics approval

This study was approved by the Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (No: E2020605). Participants gave informed consent to participate in the study before taking part.

Footnote

JL, YC, YL and LY contributed equally.

Contributors Conceptualization: JL. Data curation: ZY, JL, YC, YL. Formal analysis: JL, YL, LY, ZZ. Funding acquisition: ZY, JC, ZW. Investigation: JL, YC, YL, LY, ZZ, CLu. Methodology: JL, ZZ. Project administration: ZY, JC. Resources: ZY, JC, ZW, YC, CLi. Software: JL, ZZ. Supervision: ZY, JC, ZW, ZT, YG. Validation: YC, YL, LY. Visualization: JL. Writing—original draft: JL. Writing—review and editing: ZY, YC, YL, LY, ZZ, XW, HB, YZ, SL, MG. ZY is responsible for the overall content as the guarantor.

Funding This work was supported by the National Natural Science Foundation of China under Grant 82172674; Cancer Precision Radiotherapy Spark Program of China International Medical Foundation under Grant 2019-N-11-26; Tianjin Key Medical Discipline (Specialty) Construction Project under Grant TJYXZDXK-009A; Construction Project of Cancer Precision Diagnosis and Drug Treatment Technology (TMUCIH) under Grant ZLJZZDYYWZL08; Clinical Trial Project (TMUCIH) under Grant TZ3582022-005 and TZ3582022-011; and Wu Jieping Medical Foundation under Grant 320.6750.2020-10-96.

Competing interests None declared.

Provenance and peer review Not commissioned; externally peer reviewed.

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.

References

¹ Brenner MB, McLean J, Dialynas DP, et al. Identification of a putative second T-cell receptor. Nature 1986; 322: 145–9. doi:10.1038/322145a0

² Sun L, Su Y, Jiao A, et al. T cells in health and disease. Signal Transduct Target Ther 2023; 8: 235. doi:10.1038/s41392-023-01471-y

³ Silva-Santos B, Serre K, Norell H. γδ T cells in cancer. Nat Rev Immunol 2015; 15: 683–91. doi:10.1038/nri3904

⁴ Park JH, Lee HK. Function of γδ T cells in tumor immunology and their application to cancer therapy. Exp Mol Med 2021; 53: 318–27. doi:10.1038/s12276-021-00576-0

⁵ Christmas SE, Brew R, Deniz G, et al. T-cell receptor heterogeneity of gamma delta T-cell clones from human female reproductive tissues. Immunology 1993; 78: 436–43.

⁶ Sung H, Ferlay J, Siegel RL, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2021; 71: 209–49. doi:10.3322/caac.21660

⁷ Small W, Bacon MA, Bajaj A, et al. Cervical cancer: a global health crisis. Cancer 2017; 123: 2404–12. doi:10.1002/cncr.30667

⁸ Mabuchi S, Kawano M, Takahashi R, et al. Management of metastatic and recurrent cervical cancer. In: Shoupe D, ed. Handbook of gynecology. Cham: Springer International Publishing 2016:1-11, doi:10.1007/978-3-319-17002-2

⁹ Papadopoulou M, Sanchez Sanchez G, Vermijlen D. Innate and adaptive γδ T cells: how, when, and why. Immunol Rev 2020; 298: 99–116. doi:10.1111/imr.12926

¹⁰ Liu J, Lichtenberg T, Hoadley KA, et al. An integrated TCGA pan-cancer clinical data resource to drive high-quality survival outcome analytics. Cell 2018; 173: 400–16. doi:10.1016/j.cell.2018.02.052

¹¹ Cancer Genome Atlas Research Network, Weinstein JN, Collisson EA, et al. The cancer genome atlas pan-cancer analysis project. Nat Genet 2013; 45: 1113–20. doi:10.1038/ng.2764

¹² Pizzolato G, Kaminski H, Tosolini M, 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 U S A 2019; 116: 11906–15. doi:10.1073/pnas.1818488116

¹³ Huang H, Nie C-P, Liu X-F, et al. Phase I study of adjuvant immunotherapy with autologous tumor-infiltrating lymphocytes in locally advanced cervical cancer. J Clin Invest 2022; 132: e157726. doi:10.1172/JCI157726

¹⁴ Clark K, Vendt B, Smith K, et al. The cancer imaging archive (TCIA): maintaining and operating a public information repository. J Digit Imaging 2013; 26: 1045–57. doi:10.1007/s10278-013-9622-7

¹⁵ Lucchesi FR, Aredes ND. The cancer genome atlas cervical squamous cell carcinoma and endocervical adenocarcinoma collection (TCGA-CESC). The Cancer Imaging Archive 2020. doi:10.7937/K9/TCIA.2016.SQ4M8YP4

¹⁶ Chen T, Chen X, Zhang S, et al. The genome sequence archive family: toward explosive data growth and diverse data types. Genomics, Proteomics & Bioinformatics 2021; 19: 578–83. doi:10.1016/j.gpb.2021.08.001

¹⁷ CNCB-NGDC Members and Partners, Xue Y, Bao Y, et al. Database resources of the national genomics data center, China national center for bioinformation in 2022. Nucleic Acids Res 2022; 50: D27–38. doi:10.1093/nar/gkab951

¹⁸ Yuan Z. Data from: unraveling the influence of the cervical cancer tumor immune microenvironment on treatment response and prognosis: insights from single-cell transcriptome sequencing. Genome Sequence Archive 2023. Available: https://ngdc.cncb.ac.cn/gsa-human/browse/HRA005178

¹⁹ Hao Y, Hao S, Andersen-Nissen E, et al. Integrated analysis of multimodal single-cell data. Cell 2021; 184: 3573–87. doi:10.1016/j.cell.2021.04.048

²⁰ McGinnis CS, Murrow LM, Gartner ZJ. Doubletfinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst 2019; 8: 329–37. doi:10.1016/j.cels.2019.03.003

²¹ Liu C, Li X, Huang Q, et al. Single-cell RNA-sequencing reveals radiochemotherapy-induced innate immune activation and MHC-II upregulation in cervical cancer. Signal Transduct Target Ther 2023; 8: 44. doi:10.1038/s41392-022-01264-9

²² Hu C, Li T, Xu Y, et al. Cellmarker 2.0: an updated database of manually curated cell markers in human/mouse and web tools based on scRNA-Seq data. Nucleic Acids Res 2023; 51: D870–6. doi:10.1093/nar/gkac947

²³ Franzén O, Gan L-M, Björkegren JLM. Panglaodb: a web server for exploration of mouse and human single-cell RNA sequencing data. Database (Oxford) 2019; 2019: baz046. doi:10.1093/database/baz046

²⁴ Korsunsky I, Millard N, Fan J, et al. Fast, sensitive and accurate integration of single-cell data with harmony. Nat Methods 2019; 16: 1289–96. doi:10.1038/s41592-019-0619-0

²⁵ Street K, Risso D, Fletcher RB, et al. Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics 2018; 19: 477. doi:10.1186/s12864-018-4772-0

²⁶ Trapnell C, Cacchiarelli D, Grimsby J, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol 2014; 32: 381–6. doi:10.1038/nbt.2859

²⁷ Gulati GS, Sikandar SS, Wesche DJ, et al. Single-cell transcriptional diversity is a hallmark of developmental potential. Science 2020; 367: 405–11. doi:10.1126/science.aax0249

²⁸ Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun 2021; 12: 1088. doi:10.1038/s41467-021-21246-9

²⁹ Newman AM, Steen CB, Liu CL, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol 2019; 37: 773–82. doi:10.1038/s41587-019-0114-2

³⁰ Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015; 12: 453–7. doi:10.1038/nmeth.3337

³¹ Ritchie ME, Phipson B, Wu D, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015; 43: e47. doi:10.1093/nar/gkv007

³² Aibar S, González-Blas CB, Moerman T, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods 2017; 14: 1083–6. doi:10.1038/nmeth.4463

³³ Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-Seq data. BMC Bioinformatics 2013; 14: 7. doi:10.1186/1471-2105-14-7

³⁴ Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005; 102: 15545–50. doi:10.1073/pnas.0506580102

³⁵ Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb) 2021; 2: 100141. doi:10.1016/j.xinn.2021.100141

³⁶ Mermel CH, Schumacher SE, Hill B, et al. Gistic2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol 2011; 12: R41. doi:10.1186/gb-2011-12-4-r41

³⁷ Mayakonda A, Lin D-C, Assenov Y, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018; 28: 1747–56. doi:10.1101/gr.239244.118

³⁸ Lapuente-Santana Ó, van Genderen M, Hilbers PAJ, et al. Interpretable systems biomarkers predict response to immune-checkpoint inhibitors. Patterns (N Y) 2021; 2: 100293. doi:10.1016/j.patter.2021.100293

³⁹ Shi Z, Traverso A, van Soest J, et al. Technical note: ontology-guided radiomics analysis workflow (O-RAW). Med Phys 2019; 46: 5677–84. doi:10.1002/mp.13844

⁴⁰ van Griethuysen JJM, Fedorov A, Parmar C, et al. Computational radiomics system to decode the radiographic phenotype. Cancer 2017; 77: e104–7. doi:10.1158/0008-5472.CAN-17-0339

⁴¹ Zwanenburg A, Vallières M, Abdalah MA, et al. The image biomarker standardization initiative: standardized quantitative radiomics for high-throughput image-based phenotyping. Radiology 2020; 295: 328–38. doi:10.1148/radiol.2020191145

⁴² Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw 2010; 33: 1–22.

⁴³ Patil I. “Visualizations with statistical details: the 'ggstatsplot' approach”. JOSS 2021; 6: 3167. doi:10.21105/joss.03167

⁴⁴ Tosolini M, Pont F, Poupot M, et al. Assessment of tumor-infiltrating TCRVγ9Vδ2 γδ lymphocyte abundance by Deconvolution of human cancers microarrays. Oncoimmunology 2017; 6: e1284723. doi:10.1080/2162402X.2017.1284723

⁴⁵ Lertworapreecha M, Patumraj S, Niruthisard S, et al. Cytotoxic function of gamma delta (gamma/delta) T cells against pamidronate-treated cervical cancer cells. Indian J Exp Biol 2013; 51: 597–605.

⁴⁶ Van Hede D, Polese B, Humblet C, et al. Human papillomavirus oncoproteins induce a reorganization of epithelial-associated γδ T cells promoting tumor formation. Proc Natl Acad Sci U S A 2017; 114: E9056–65. doi:10.1073/pnas.1712883114

⁴⁷ Sun R, Limkin EJ, Vakalopoulou M, et al. A radiomics approach to assess tumour-infiltrating CD8 cells and response to anti-PD-1 or anti-PD-L1 immunotherapy: an imaging biomarker, retrospective multicohort study. Lancet Oncol 2018; 19: 1180–91. doi:10.1016/S1470-2045(18)30413-3

⁴⁸ Su G-H, Xiao Y, Jiang L, et al. Radiomics features for assessing tumor-infiltrating lymphocytes correlate with molecular traits of triple-negative breast cancer. J Transl Med 2022; 20: 471. doi:10.1186/s12967-022-03688-x

⁴⁹ Ramlee S, Hulse D, Bernatowicz K, et al. Radiomic signatures associated with CD8+ tumour-infiltrating lymphocytes: a systematic review and quality assessment study. Cancers (Basel) 2022; 14: 3656. doi:10.3390/cancers14153656

⁵⁰ Jin H, Liu Z. A benchmark for RNA-Seq deconvolution analysis under dynamic testing environments. Genome Biol 2021; 22: 102. doi:10.1186/s13059-021-02290-6