Article Text
PDF
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.
- Genital Neoplasms, Female
- Tumor Microenvironment
- T-Lymphocytes
- Immunotherapy
- Biostatistics
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.
This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See http://creativecommons.org/licenses/by-nc/4.0/.
Statistics from Altmetric.com
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.1 They mainly reside in peripheral mucosal barriers. These cells play crucial roles in maintaining physiological homeostasis and orchestrating immune responses during pathological conditions.2 γδ 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.2
The successful isolation of γδ T cells from the human cervix was accomplished three decades ago.5 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.6 Among these, approximately 90% of new cases and fatalities related to cervical cancer are concentrated in less developed regions.7 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.8 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.9 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.
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 GSE12822312 and GSE19007513 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 Archive16 of the National Genomics Data Center17 with reference number HRA005178.18
Supplemental material
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)19 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).20 Linear and graph-based dimensionality reduction was performed, and cells were clustered using the Louvain algorithm as described in the online supplemental material.
Supplemental material
Annotation of cell clusters
Cell clusters were annotated based on a previous study21 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 Harmony24 (V.0.1.1). Trajectory analysis was conducted using Slingshot25 (V.2.8). Monocle26 (V.2.28) was also used for trajectory exploration and visualization. CytoTRACE27 (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 CellChat28 (V.1.6.1). Cell–cell contact and secreted signaling were observed separately.
Estimation of cell infiltration from bulk transcriptome data
The CIBERSORTx29 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,12 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 set30 using a similar approach.
Differential expression and gene set analyses
Differential expression analyses were performed using the limma31 package (V.3.56.2) or a Wilcoxon rank-sum test. Gene set analyses were performed based on the AUCell package,32 gene set variation analysis (GSVA)33 or gene set enrichment analysis (GSEA).34 The clusterProfiler35 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). GISTIC36 (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 Maftools37 (V.2.16).
In silico estimation of response to ICIs
The EaSIeR38 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,39 which operates based on PyRadiomics40 (V.3.7). The majority of the extracted features were in accordance with the guidelines set forth by the Image Biomarker Standardization Initiative.41 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 package42 (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 χ2 test or Fisher’s exact test were performed in R (V.4.3.1) or MedCalc Statistical Software (V.19.6.1). The ggstatsplot43 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.
The infiltration level of γδT cells in tumor tissue was associated with an overall favorable prognosis for patients with cervical cancer. (A) Example immunohistochemical staining plot of γδ TCR+ cells in cervical cancer tissues under 200×high-power field (HPF). (B) Example immunohistochemical staining plots with high or low abundance of γδ TCR+ cells in the 400×field of view. (C) Bar graph showing the frequency distribution of the number of γδ TCR+ cells per 200×HPF. (D) Box plots demonstrating the number of γδ TCR+ cells per 200×HPF in patients with different survival statuses. (E) Kaplan-Meier curves show overall survival (OS) for patients with different γδ TCR+ cell counts per 200×HPF. (F) t-distributed stochastic neighbor embedding (t-SNE) plots of the γδ T-cells reference data set, colored by Vδ1/Vδ2 T-cell subtypes (left) or donors (right). (G) t-SNE-based feature plots show the expression of typical αβ/γδ T cells marker genes. (H) Kaplan-Meier curves show OS (left) and progression-free interval (right) for patients with high or low γδ 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 study44 (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.
Supplemental material
Supplemental material
Supplemental material
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.
Supplemental material
Supplemental material
The focal somatic mutations and tumor microenvironment characteristics of cervical carcinoma tissue were associated with the infiltration level of γδ 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 γδ T-cell high and low infiltration level groups. OR>1, higher SNV frequency for patients in the γδ 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 γδ 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 γδ 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 γδ 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 γδ T-cell infiltration grouping for each sample. (G) Top 10 upregulated and downregulated pathways in patients with a high abundance of γδ 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 γδ 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.
Supplemental material
γδ 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.44 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.
Supplemental material
Supplemental material
Chasing down γδ 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. γδ 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.
Dynamic transformation of γδ T cells and their communication with other cells in the tumor microenvironment of cervical cancer. (A) UMAP plot showing γδ T cells in the two states after removal of the between-sample batch effect. (B) Volcano plot demonstrating differences in the gene expression of γδ T cells in the two states. (C) The differentiation potential scores of γδ T cells in GNLY− status were significantly higher than those in GNLY+ status. (D) Evolutionary trajectories of γδ T cells depicted on a UMAP plot. Cells are stained in pseudotemporal order. (E) Density plots of γδ 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 γδ 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+ αβ T cells and γδ T cells. Colors represent the probability of contact. The sizes of the dots represent p values. (H) Dot plot showing secreted signaling communications from γδ T cells to CD8+ αβ 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.
Supplemental material
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.
Patients with high levels of γδ 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 γδ 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 γδ TCR constant region. (F) AUCell score for the γδ 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−4), CTLA-4 (p=7.6×10−4), 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−22), CTLA-4 (r=0.75, p=2.2×10−57), and IFN-γ (r=0.77, p=5.9×10−62) 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−3) (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.
Supplemental material
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 R2 value of the model was 0.625, the adjusted R2 value was 0.528, and the p value was 6.067×10−5 (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.
MRI-based radiomics analysis established a linear regression model predicting γδ 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. λ, the weight of the L1 norm. The two vertical dashed lines represent λmin and λ1se. λmin, the value of λ that results in the smallest mean value of the squared error of the regression model. λ1se, the value of λ that results in the simplest model within an SE of λ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 γδ T-cell abundance measured by the transcriptome and the fitted value (ie, the predicted γδ 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 γδ 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.45 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.46 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-based47 or MRI-based48 radiomic biomarkers of tumor-infiltrating CD8+ T cells, although the radiomic signatures they reported were heterogeneous with limited reproducibility.49 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 volunteers12 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 study50 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 CIBERSORTx29 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.9 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.
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
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.
Acknowledgments
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.
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
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.