Article Text

Original research
Hypoxia inhibits the iMo/cDC2/CD8+ TRMs immune axis in the tumor microenvironment of human esophageal cancer
  1. Chuanqiang Wu1,2,
  2. Huan Yu1,
  3. Fuxiang Liang1,
  4. Xiancong Huang1,3,
  5. Bin Jiang4,
  6. Zhiling Lou1,
  7. Yafei Liu1,
  8. Zixiang Wu1,
  9. Qi Wang1,
  10. Hong Shen5,
  11. Ming Chen6,
  12. Pin Wu1,7 and
  13. Ming Wu1,2
  1. 1 Department of Thoracic Surgery, The Second Affiliated Hospital Zhejiang University School of Medicine,Zhejiang University, Hangzhou, Zhejiang Province, People's Republic of China
  2. 2 Laboratory of Clinical Research Center of Zhejiang Province, The Second Affiliated Hospital Zhejiang University School of Medicine, Zhejiang University, Hangzhou, Zhejiang Province, People's Republic of China
  3. 3 Department of Thoracic Surgery, Zhejiang Cancer Hospital, Hangzhou, Zhejiang Province, People's Republic of China
  4. 4 Department of Thoracic Surgery, Shandong Provincial Hospital, Jinan, Shandong Province, People's Republic of China
  5. 5 Department of Medical Oncology, The Second Affiliated Hospital Zhejiang University School of Medicine, Zhejiang University, Hangzhou, Zhejiang Province, People's Republic of China
  6. 6 Department of Bioinformatics, College of Life Sciences, Zhejiang University, Hangzhou, Zhejiang Province, People's Republic of China
  7. 7 Key Laboratory of Tumor Microenvironment and Immune Therapy of Zhejiang Province, The Second Affiliated Hospital Zhejiang University School of Medicine, Zhejiang University, Hangzhou, Zhejiang, People's Republic of China
  1. Correspondence to Dr Pin Wu; pinwu{at}zju.edu.cn; Professor Ming Wu; iwuming22{at}zju.edu.cn

Abstract

Background Esophageal cancer (ESCA) is a form of malignant tumor associated with chronic inflammation and immune dysregulation. However, the specific immune status and key mechanisms of immune regulation in this disease require further exploration.

Methods To investigate the features of the human ESCA tumor immune microenvironment and its possible regulation, we performed mass cytometry by time of flight, single-cell RNA sequencing, multicolor fluorescence staining of tissue, and flow cytometry analyses on tumor and paracancerous tissue from treatment-naïve patients.

Results We depicted the immune landscape of the ESCA and revealed that CD8+ (tissue-resident memory CD8+ T cells (CD8+ TRMs) were closely related to disease progression. We also revealed the heterogeneity of CD8+ TRMs in the ESCA tumor microenvironment (TME), which was associated with their differentiation and function. Moreover, the subset of CD8+ TRMs in tumor (called tTRMs) that expressed high levels of granzyme B and immune checkpoints was markedly decreased in the TME of advanced ESCA. We showed that tTRMs are tumor effector cells preactivated in the TME. We then demonstrated that conventional dendritic cells (cDC2s) derived from intermediate monocytes (iMos) are essential for maintaining the proliferation of CD8+ TRMs in the TME. Our preliminary study showed that hypoxia can promote the apoptosis of iMos and impede the maturation of cDC2s, which in turn reduces the proliferative capacity of CD8+ TRMs, thereby contributing to the progression of cancer.

Conclusions Our study revealed the essential antitumor roles of CD8+ TRMs and preliminarily explored the regulation of the iMo/cDC2/CD8+ TRM immune axis in the human ESCA TME.

  • Esophageal Cancer
  • Tumor infiltrating lymphocyte - TIL
  • Dendritic
  • Immunosuppression

Data availability statement

Data are available in a public, open access repository. Data are available upon reasonable request. All data relevant to the study are included in the article or uploaded as supplementary information. The transcriptome sequence data of ESCA were obtained from the Genome Sequence Bank of Beijing Institute of Genomics, Chinese Academy of Sciences (https://ngdc.cncb.ac.cn/search/all?&q=HRA004010). The single-cell RNA sequencing data of ESCA were downloaded from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE160269).

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

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

Statistics from Altmetric.com

Request Permissions

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

WHAT IS ALREADY KNOWN ON THIS TOPIC

  • As a malignant tumor, the treatment and prognosis of patients with esophageal cancer (ESCA) are still not ideal. The occurrence and development of ESCA are closely related to immune dysregulation, however, our understanding of its specific regulatory mechanisms is still very limited.

WHAT THIS STUDY ADDS

  • We have depicted the single-cell immune landscape of ESCA ecosystem and identified the key immune cell subpopulation CD8+ tissue resident memory T cells (TRMs) related to disease progression. We also revealed the heterogeneity of CD8+ TRMs in function and development. In addition, we found that in the tumor microenvironment of ESCA, the hypoxic environment impairs the intermediate monocytes/conventional dendritic cells/CD8+ TRMs immune axis, thereby reducing the antitumor ability of CD8+ TRMs and further promoting the progression of cancer.

HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY

  • These findings may influence the immunotherapeutic strategies for ESCA. To improve the immunotherapeutic effect of ESCA, we need to further study how to achieve this goal by regulating the immune environment. These studies provide a theoretical basis for the development of new treatment methods.

Background

Despite advances in the management and treatment of esophageal cancer (ESCA), the overall prognosis for patients with this disease remains poor.1 A comprehensive understanding of the differences in immune cell subsets between paracancerous tissue and tumor tissues, and an exploration of the potential functional mechanism of key subgroups of immune cells, will help us to guide ongoing clinical studies and construct new immunotherapy evaluation models

Smoking, alcoholism, prolonged intake of hot beverages and nitrite products, Helicobacter pylori infection, poor oral hygiene, and genetic mutations are common risk factors for ESCA, which can all lead to chronic inflammation.1–4 Recent studies have highlighted the crucial role of chronic inflammation in the progression of ESCA.5 6 Myeloid immune cells are the body’s primary response cells to inflammation7 and have the potential to directly or indirectly influence the progression of tumors.8–10 Insights gleaned from numerous studies on the tumor microenvironment (TME) underscore the pivotal role of the myeloid immune cells.11 12 The TME plays a pivotal role not only in determining the fates of T cells, but also in shaping the overall structure of the TME.13 However, our understanding of the specific mechanisms through which they operate remains limited.

CD8+ memory T (Tm) cells can be broadly classified into three groups: circulating CD8+ T central memory cells (CD8+ Tcm), CD8+ T effector memory cells (CD8+ Tem), and CD8+ tissue-resident memory T cells (CD8+ TRMs).14 CD8+ TRM cells, often found in barrier tissues such as the skin and mucosa, play a pivotal role in maintaining immune balance, defending against pathogenic infections, and monitoring potential malignancies. This role has been substantiated by previous research, which offers compelling evidence of their critical function.15–18 However, the functional characteristics and regulatory mechanisms of CD8+TRMs remain to be explored.

Here, we employed cytometry by time of flight (CyTOF) to characterize the immune profiles of tumors and corresponding paracancerous tissue samples from 20 patients with ESCA. Specifically, we identified CD8+ TRMs as the primary subset of tumor-killing cells in ESCA that are associated with disease progression. Additionally, we revealed significant heterogeneity in the ontogeny and function of CD8+ TRMs. By analyzing the intrinsic mechanisms that contribute to the formation of the ESCA TME, we identified the iMo/cDC2/CD8+ TRMs axis, which plays a crucial role in supporting the proliferation of CD8+ TRMs within the ESCA TME. Finally, our findings demonstrated that hypoxia is a principal modulator of disease progression in ESCA through the perturbation of the iMo/cDC2/CD8+ TRMs immune axis. In conclusion, our study highlights fundamental changes in the TME of patients with human ESCA and preliminarily explores the critical role of hypoxia-induced dysfunction of the iMo/cDC2/CD8+ TRMs axis in tumor progression.

Material and methods

Surgical sample collection

Human specimens were obtained from the Second Affiliated Hospital of Zhejiang University (SAHZU) and Zhejiang Cancer Hospital (ZJCH). The tumor tissues and corresponding normal tissues of patients with ESCA were obtained from surgical resections for analysis using CyTOF, flow cytometry and immunohistochemistry. The specimens were collected within 5 min of surgical resection and placed in a cold tissue preservation solution (Miltenyi Biotec, CAT#130-100-008). Ethics approval was obtained from the Medical Ethics Committee of SAHZU and ZJCH, and all individual patients provided written informed consent. Patient information is summarized in online supplemental table 1.

Supplemental material

Supplemental material

CyTOF sample preparation

Fresh tumor and normal tissues were mechanically cut into small pieces and enzymatically digested in RPMI 1640 (Gibco, CAT#11875093) supplemented with 1 mg/mL collagenase type I (Solarbio, CAT# C8140), 1 mg/mL collagenase type IV (Solarbio, CAT# C8160), 1% penicillin–streptomycin (Procell, CAT# PB180120), and 4% fetal bovine serum (Procell, CAT# 164210–50) using the GentleMACS system (Miltenyi Biotec, CAT#130-096-427). After complete tissue dissociation, the digestion solution was filtered into a centrifuge tube using a 70 µm sieve. The filtrate was centrifuged, the pellet was washed twice, and the pellets was resuspended in 1×phosphate-buffered saline (PBS) (Procell, CAT# PB180327). From the obtained single-cell suspensions, 3×106 cells per sample were added to a 1.5 mL centrifuge tube for live and dead cell staining and mass spectrometry staining experiments.19

CyTOF sample staining and data acquisition

After sample preparation, single-cell suspensions were stained with metal-conjugated antibodies. The antibodies used for mass cytometry are shown in online supplemental table 3. A total of 3×106 cells per sample were stained for cell surface markers at 4°C for 30 min, stained with 1 mM cisplatin (Solarbio) for 5 min at room temperature, and fixed using transcription factor buffer (BD Biosciences), followed by intracellular staining at 4°C for 60 min. Then, the cells were resuspended in deionized water containing Ce140, Eu151, Eu153, Ho165, and Lu175 (Fluidigm). After normalizing and randomizing values near zero using Helios software, the FCS files were uploaded for analysis.

Supplemental material

CyTOF data preprocessing and quality control

Three preprocessing steps are needed for CyTOF data: bead normalization to eliminate errors caused by instrument signal drift, debarcoding to obtain cell data from a single cell, and manual gating to isolate effective target single cells (removing cell debris, dead cells, cells, non-target cells, etc). After preprocessing, the sample data were used for subsequent analysis.

Data analysis and visualization

FlowJo software (BD Biosciences, V.10.6.2) was used to obtain the target cell population, with an average of 20,000 cells randomly selected from each sample, resulting in a total of 800,000 cells. R was used to perform PhenoGraph clustering analysis20 based on machine learning. High-dimensional data were visualized using dimensionality reduction methods such as t-distributed stochastic neighbor embedding, heatmaps, density plots, and bar plots to intuitively capture practical conclusions with biological significance.

Preprocessing of single-cell RNA-sequencing data

The Unique Molecular Identifiers (UMI) data for T cells and myeloid cells were downloaded from GSE160269.21 Doublets in the original data were detected and removed using artificial nearest neighbors based on DoubletFinder (V.2.0.3).22 The processed data were imported into the Seurat package (V.4.1.0) for quality filtering and downstream analysis. For quality filtering, low-quality cells (<400 genes/cell and >10% mitochondrial genes) were removed. Ultimately, 57,215 cells with a median of 1,170 detected genes per cell were retained for downstream analysis.

PCA, clustering, and annotation

We identified genes with highly variable expression (HVG) based on average expression and dispersion level thresholds using the “FindVariableGenes” function with default settings. The top 2,000 HVGs were used for further data processing. The normalized gene expression was converted into Z-scores (values centered at 0 and variance 1) using the “ScaleData” function, and downstream Principal Component Analysis (PCA) dimensionality reduction was performed using the “RunPCA” function. Cell clusters were identified using a clustering algorithm based on shared nearest neighbor modular optimization with “FindClusters”. We set the resolution parameter (Res) to 0.5 and k to 15 for the k-nearest neighbor algorithm for most clustering analyses. Differential expression analysis between clusters was performed using the Wilcoxon test implemented in “FindAllMarkers” with Bonferroni correction of p values. We annotated the cell clusters according to these marker genes and visualized the clustering results and gene expression in a Uniform Manifold Approximation and Projection (UMAP) scatter plot using “RunUMAP”. Among these cells, the 29,340 CD8+ T cells were categorized into three clusters with Res=0.15. The 10,027 tissue-resident T cells were grouped into five clusters with Res=0.25. The 9,992 myeloid cells were grouped into nine clusters with Res=0.25.

Gene set variation analysis

Pathway activity scores were calculated for individual CD8+ TRMs clusters using the gene set variation analysis (GSVA) method in the R package GSVA. The GSVA method estimates gene set scores by computing enrichment statistics for each gene set in each cluster. The gene sets were loaded from Bioconductor (https://www.bioconductor.org/) using the GSVA package (V.1.42.0). Pathways with significant differences in activity scores were selected using the LIMMA package (V.3.50).

SCENIC analysis

The regulon transcriptional activity analysis was performed by following the default parameters of the SCENIC pipeline in SCENIC (V.1.2.4),23 following default parameters SCENIC pipeline, and raw expression matrices were used as input to analyze activated regulons in each CD8+ TRMs subset. Briefly, we identified regulons using the RcisTarget package (V.1.14.0). The activity of the regulons per cell was then quantified by the AUCell package (V.1.16.0). The regulon area under the curve thresholds for determining binary activity were manually assessed and adjusted based on the expressed regulon. Moreover, differentially activated regulons in each subset of CD8+ TRMs were identified by the Wilcoxon test using the other subsets as controls. Finally, the corresponding p values were corrected for False Discovery Rate (FDR) using the Benjamini-Hochberg method.

Pseudotime trajectory analysis

Monocle (V.2)24 was used to determine possible lineage differentiation of the identified CD8+ T cells and tissue-resident T cells using default parameters. First, cell identification information and gene expression data for each cell were extracted and integrated for subsequent analysis. Then, the “estimateSizeFactors”, “estimateDispersions”, and “FindAllMarkers” functions with default parameters were used to construct statistical models to characterize the data. Genes with a p value<0.05 were filtered out for further analysis and plotting. We used the DDRTree method to reduce the dimensionality via the “reduceDimension” function.

Cell–cell interaction network

The intercellular communication inference, visualization and analysis of the single-cell RNA sequencing (scRNA-seq) data were performed using the CellChat package25 of R software (V.1.1.3). The crosstalk between cells was analyzed via “CellChat” as follows: the “CellChat” object was constructed through the “createCellChat” function; cellular communication network inference was performed through the “computeCommunProb”, “computeCommunProbPathway” and “aggregate net” functions; and visualization was performed through the “netVisual_aggregate”, “netVisual_signalingRole”, and “netVisual_bubble” functions.

Weighted gene coexpression network analysis and Gene Ontology enrichment analysis

Complementary DNA library construction, purification, and RNA sequencing (RNA-seq) were performed on the BGISEQ platform at BGI Technology (Wuhan, China), following standard protocols. Standard bioinformatics analyses were performed using the BGI platform, and expression matrices were used for further analysis. Signed coexpression networks were constructed for bulk ESCA RNA-seq data using the weighted gene coexpression network analysis (WGCNA, V.1.71) package in R according to the “WGCNA” manual. A correlation heatmap of modules and sample features and a hierarchical clustering map of modules and sample features were constructed to determine the critical modules with the highest correlation coefficients with the occurrence of ESCA. Gene Ontology enrichment analysis for selected genes was carried out using the DAVID functional annotation database (https://david.ncifcrf.gov).

T-cell receptor data analysis

T-cell receptor (TCR) data were processed by running the cellranger vdj (10x Genomics, V.2.1.1) with–reference cellranger-vdj-GRCh38-alts-ensembl-2.0.0 TCR for α-chain and β-chain assembly and clonetyping. Only TCRs with full α and β chain sequences were included in the analysis. We obtained a paired transcriptome and TCR sequence for 8,773 CD8+ TRMs (online supplemental table 8). To assign two cells to the same clonotype, they must share α and β sequences. If a clonotype was present in at least two cells, the cells carrying that clonotype were considered clonal, and the number of cells with this α-β pair indicated the degree of clonality of that clonotype.

Supplemental material

Definition of naïvety, tissue residence, exhaustion, cytotoxicity, and proliferation scores

To construct the cell characteristic evaluation system, we obtained previously identified tissue residence, exhaustion, and cytotoxicity signature gene lists (online supplemental table 6) and calculated signature scores. We defined the naïve gene score of T cells using the Z-score of the mean expression of four well-defined naïvety markers (CCR7, TCF7, LEF1, and SELL). We assessed the proliferative properties of the clusters by the expression of different cell cycle genes (online supplemental table 6).16 26 27

Supplemental material

Myeloid cell definition of CyToF

Neutrophils (CD45+CD3CD56CD20CD16+CD66b+), eosinophils (CD45+CD3CD56CD20CD16CD66b+), iMo (CD45+CD3CD56CD20CD14++CD16+), cMo (CD45+CD3CD56CD20CD14+ CD16), DN-DC (CD45+CD3CD56CD20CD11c+CD141CD1c), cDC1 (CD45+CD3CD56CD20CD141+CD1cCD11c+), cDC2 (CD45+CD3CD56CD20CD11c+CD1c+CD141), HLA-DR+CD141+CD11ccells, HLA-DR+CD11cCD141- cells, PD-L1loPDL2loCells, PD-L1PDL2lo cells, HLA-DRCD14CD16+ cells, CD11C+PDL2lo cells and unknown cells.

Multicolor immunofluorescence

Human tissue microarray (TMA) specimens from 77 patients (detailed information is provided in online supplemental table 2) were purchased from a commercial tissue distributor (Superbiotek, CAT#ESC77). In brief, the TMA sections were first placed in an oven at 70°C for 1 hour and then deparaffinization, rehydration, and subjected to antigen retrieval. The TMA was preincubated with immunofluorescence Blocking Buffer (Cell Signaling Technology, CAT#12411) to block non-specific sites. Sections were then incubated with anti-CD3 (Biolynx, CAT# BX50022-C3) antibody at room temperature for 60 min, horseradish peroxidase (HRP) secondary antibody and Tyramide Signal Amplification (TSA) fluorescence dye were added, and then single-color staining for CD8 (Biolynx, CAT# BX50036-C3, CD103 (Abcam, CAT#ab129202) or CD11C (Newmarket Scientific, CAT#GB11059), CD1C (Abcam, CAT#ab246520) expression were performed. Finally, the sections were subjected to 4',6-diamidino-2-phenylindole (DAPI) nuclear staining, mounting and image acquisition.28

Supplemental material

Flow cytometry and cell sorting

Minced tissue was subjected to enzymatic dissociation as described above. For surface staining, the cells were stained in PBS containing 0.5% Bovine Serum Albumin (BSA). Following surface staining, transcription factor staining buffer was used to fix and permeabilize the cells, and then the cells were incubated with the Ki67 antibody at room temperature for 50 min. Samples were analyzed on an LSRII Fortessa (BD Biosciences) or sorted to 95% purity on MoFlo Astrios EQs (BECKMAN). Flow cytometry analysis was performed using FlowJo software (Tree Star). For detailed antibody information, see online supplemental table 3.29

Co-culture, hypoxic culture and induction differentiation experiments

Following sorting, CD8+ TRMs and cDC2s were co-cultured for 48 hours at a 1:5 ratio in Roswell Park Memorial Institute) RPMI 1640 medium containing 10% FBS and 1% penicillin–streptomycin. For the hypoxic culture experiment, the treatment-naïve ESCA tumor single-cell suspensions were treated with or without 100 µM cobalt chloride for 48 hours in RPMI 1640 medium supplemented with 10% FBS and 1% penicillin–streptomycin. CD14+CD16+ intermediate monocytes (iMos) were sorted from whole blood collected from patients with treatment-naive ESCA. After sorting, the iMos were cultured for 2 days in the presence or absence of 50 ng/mL recombinant human interleukin (IL)-4 and 100 ng/mL recombinant human Granulocyte-Macrophage Colony-Stimulating Factor(GM-CSF) in complete RPMI 1640 medium. Differentiable CD8+ TRMs (dTRMs) and proliferative CD8+ TRMs (pTRMs) were sorted from tumor tissues collected from patients with treatment-naive ESCA during surgery. After sorting, the dTRMs and pTRMs were cultured for 3 days in the presence or absence of 20 ng/mL recombinant human TGF-β (Abcam, CAT#ab50036) or 100 ng/mL recombinant human IL-2 (Abcam, CAT#ab269149) in complete RPMI 1640 medium. Following induction, the sample was stained with antibodies and analyzed on an LSRII Fortessa.

Apoptosis assays

Cell apoptosis was detected according to the manual provided with the reagent kit (LinkBio, CAT#AT101). In brief, harvested cells were cultured under normoxic or hypoxic conditions and washed with prechilled 1×PBS solution. After centrifugation, the PBS was discarded, and the cells were resuspended in 1×binding buffer. Then, 5 µL of annexin V-FITC and 10 µL of Propidium Iodide (PI) were added to each sample, which was then gently vortexed and incubated at room temperature for 5 min and analyzed on an LSRII Fortessa.

Immunohistochemical staining

The paraffinized ESCA sections were baked in an oven at a constant temperature for 1 hour, followed by sequential immersion in different concentrations of ethanol for deparaffinization. After deparaffinization, antigen retrieval was performed in a pressure cooker, and the sections were washed with PBS. The samples were completely covered with 5% goat serum and incubated at 37°C in a humidified chamber for 30 min. The blocking solution was removed, and IL-4 (Bioss, bs-0581R) or GM-CSF (Bioss, bs-3790R) antibodies were added, followed by overnight incubation at 4°C. After removing the working antibody solution, the sections were washed with Tris Buffered Saline with Tween 20 (TBST) buffer. Diluted secondary fluorescence antibodies were then added and incubated at 37°C for 1 hour. Following the removal of the secondary antibody working solution, the samples were washed again with TBST buffer. DAPI working solution was applied to the washed samples, which were then incubated at room temperature in the dark for 10 min. After adding an antifade mounting medium, the samples were observed under a fluorescence microscope, and images were captured.

Flow cytometry cell type definition

CD8+ TRMs: CD45+CD3+CD8a+CD69+CD103+

cDC2: LinCD45+HLA-DR+CD11C+CD1C+

iMo: LinCD45+HLA-DR+CD14++CD16+

Animal studies

Animal studies were conducted in compliance with protocols approved by the Zhejiang University Institutional Animal Care and Use Committee. The CD103 conditional knockout C57 (C57BL/6)/JGpt-Itgaeem1Cflox/Gpt mice and CD8 tool C57/JGpt-Tg(Cd8a-iCre)19/Gpt mice were were purchased from GemPharmatech (Nanjing, China). First, these two lines of mice were crossed to generate a CD8+ T-cell-specific CD103 conditional knockout mouse model (CD8creTg/Tg; ITGAEfl/fl). Subsequently, 6-week-old C57 mice were randomly divided into two groups: a control group and a tamoxifen-induced group. To obtain ITGAE conditional knockout CD8creTg/Tg; ITGAE△/△C57 mice, 2 mg/25 g (body weight) tamoxifen (Sigma-Aldrich, CAT#T2859) was injected into the mice every other day.30 The mice in the control group were treated in the same in, except that the volume of tamoxifen was replaced with corn oil. For the subcutaneous tumor formation assay CD8creTg/Tg; ITGAEfl/fl C57 mice or CD8creTg/Tg; ITGAE△/△ C57 mice received 5×105 AKR cells per mouse. 10 days later, tumor volume growth was monitored by measuring tumor size three times per week with a caliper, and body weight was assessed. Six weeks after tumor inoculation, the mice were sacrificed, and the tumors were dissected and weighed.

Statistical analysis

Statistical analysis was performed using the R, GraphPad Prism V.9.0 software and Microsoft Excel computer programs. The analysis mainly used a statistical significance test method to compare the percentages of cell subpopulations obtained by clustering in different sample groups. The methods used in this study included a two-sided t-test or Wilcoxon rank-sum test (applied when data did not meet the assumptions of normality) for testing. When p<0.05, the data between different groups were considered to be significantly differences. Data with significant differences are represented by asterisks: *indicates p<0.05, **indicates p<0.01, ***indicates p<0.001, and ****indicates p<0.0001; “no significant difference” is denoted by “ns”.

Results

The human ESCA immune landscape

To determine the characteristics of the Tumor Immune Microenvironment (TIME) and its potential regulatory mechanisms in human ESCA, we conducted several analyses including CyTOF, scRNA-seq, and flow cytometry, as well as in vitro co-culture experiments on a cohort of patients with treatment-naïve ESCA (figure 1A). We discovered significant differences in the expression of immunologically-related molecules between tumors and adjacent cancer tissues through CyTOF analysis (online supplemental figure S1A, B). Using immune cell lineage markers, we identified 42 immune cell clusters in the human ESCA TIME (online supplemental figure S1C, D), which were then grouped into seven immune cell subgroups based on the expression characteristics of cell markers (figure 1B, online supplemental figure S1E, F). We found that the largest subgroup was T cells, and the second largest subgroup was myeloid immune cells (figure 1C), with variability in immune cell infiltration observed among different patients (online supplemental figure S1G). We also observed a significant alteration in the TIME cell composition, where the proportion of granulocytes and natural killer-like cells increased in tumor samples (p<0.05), while the proportions of other myeloid cells (p<0.001) and dendritic cells (DCs) (p<0.05) were decreased compared with those in normal tissues (figure 1D). The frequencies of T cells, B cells, and Mons/Macs (monocytes/macrophages) in tumors were found to be similar to those in paired normal tissues (figure 1D). We evaluated the relationships between immune cell subgroups and clinicopathologic features but found no significant correlations (figure 1E).

Supplemental material

Supplemental material

Figure 1

Single-cell immune landscapes of ESCA tumor ecosystems. (A) Study overview. (B) t-SNE visualization of 800,000 merged events from adjacent normal tissues (normal, n=20) and tumor tissues (tumor, n=20). Each dot represents one cell, and each color represents one cluster. (C) The bar graph represents the average proportions of immune cell subgroups in normal tissue (normal) and tumor tissue (tumor). The color corresponds to the cluster type. The three-line table below the bar graph contains the specific average proportion values. (D) Comparison of the frequencies of various immune cell subgroups in normal and tumor tissues. Each point represents a sample (n=20), tumors are indicated in red, and normal tissues are indicated in blue. The error bars represent the mean±SD. (*p<0.05, **p<0.01, ***p<0.001, ****p<0.0001, ns (no significance)). (E) Relationships between the frequencies of various immune cell subgroups and clinicopathological parameters in patients with ESCA. Different colors indicate different groups. Box plot: whiskers represent the minimum to maximum values. Bounds of boxes represent the 25th and 75th percentiles, and centers represent the mean±SD (n=20). CyTOF, cytometry by time of flight; DC, dendritic cell; ESCA, esophageal cancer; NK, natural killer; scRNA-seq, single-cell RNA sequencing; TMA, tissue microarray; TNM, tumor/node/metastasis; t-SNE, t-distributed stochastic neighbor embedding.

CD8+ TRMs negatively correlate with ESCA progression

To identify cell clusters associated with the progression of ESCA, the clustering algorithm PhenoGraph was applied to divide the T cell into 28 populations, and t-distributed stochastic neighbor embedding (t-SNE) maps were generated to visualize the T-cell cluster distribution (figure 2A). Among these clusters, we identified 3 different DNT (CD4CD8T) cell clusters, 12 unique CD8+ T-cell clusters, and 13 distinct CD4+ T-cell clusters. The t-SNE plot revealed that, while there was considerable overlap among these T-cell clusters, there were also substantial differences between tumor and normal tissues (figure 2B). Paired comparative analysis revealed differences in CD4+ and CD8+ T-cell clusters between normal tissues and tumor tissues. Specifically, concerning CD4+ T-cell subsets within the tumor, we observed an increased infiltration of cluster 05 (C05) and C06, as well as decreased infiltration of C03, C07, C08, C10, and C11 (online supplemental figure S2A). Moreover, C05 was the largest differentially infiltrated cluster (figure 2C). Among various CD8+ T-cell clusters, our study revealed a significant increase in the infiltration of C25, C26, and C27 and a decrease in the infiltration of C17, C18, C19, and C22 in tumor tissues (online supplemental figure S2B). C25 was the most highly infiltrated cluster in tumor tissues, whereas C19 was the dominant cluster in normal tissues among these CD8+ T-cell clusters (online supplemental figure S2D).

Supplemental material

Figure 2

CD8+ TRMs infiltration is associated with tumor progression and overall survival of patients with ESCA. (A) t-SNE plot of T cells showing 28 major clusters, including 3 distinct DNT cell clusters, 13 distinct CD4+ T-cell clusters and 12 distinct CD8+ T-cell clusters. Each dot represents a cell and is colored based on cluster type. (B) t-SNE plots overlapping different T-cell clusters in tumor and normal tissues. Tumors are shown in blue, and normal tissues are shown in red. (C and D) t-SNE plots of CD8a, CD103, granzyme B (C) CTLA-4 and PD-1 (D) expression in T cells from all samples. The dots circled in red are CD8+TRMs. Each dot represents one cell. The cells are colored according to density. High expression is indicated in red; yellow indicates moderate expression; low expression is shown in blue. (E) Relationships between the frequencies of CD8+ T cells, circulating CD8+ T cells, and CD8+ TRMs and clinicopathological parameters in patients with ESCA. Different colors indicate different groups. The error bars represent the mean±SD. (F) Multicolor immunohistochemistry staining was performed using anti-CD3, anti-CD8a, and anti-CD103 antibodies and DAPI to validate the presence of CD8+ TRMs in ESCA tumors and their correlation with TNM stage. The following stains were used: DAPI (blue), CD3 (yellow), CD8A (pink), CD103 (red). (G) The T cells, CD8+ T cells, CD103+CD8+T cells infiltration scores in different TNM stages of ESCA. Different colors indicate different groups. The error bars represent SDs (T cells infiltration scores: T cells/total cells; CD8+T cells scores: CD8+ T cells/total T cells; CD103+CD8+T cells scores: CD103+CD8+ T cells/total T cells). (H) Kaplan-Meier survival curves of patients with ESCA grouped by the proportion of infiltrating CD8+ TRMs. The ratio of 1% CD8+ TRMs infiltration (CD103+CD8+ T cells/ total T cells) was used to delineate the groups with high and low immune infiltration. (I) Growth curves of tumors from C57 (C57BL/6) mice inoculated with AKR cells. The blue line represents the control group (CD8creTg/Tg;ITGAEfl/fl), and the red line represents the CD8+TRMs depletion group (CD8creTg/Tg;ITGAE△/△). ****p<0.0001. (J) Tumor weight was monitored after the mice were sacrificed. The blue points represent the control group (CD8creTg/Tg;ITGAEfl/fl), and the red points represent the CD8+ TRMs depletion group (CD8creTg/Tg;ITGAE△/△). Each point represents a sample (N=6). ****p<0.0001. ESCA, esophageal cancer; TNM, tumor/node/metastasis; t-SNE, t-distributed stochastic neighbor embedding; TRMs, tissue resident memory T cells; DAPI, 4’,6-diamidino-2-phenylindole; PD-1, Programmed cell death protein 1;CTLA-4, Cytotoxic T-Lymphocyte Antigen 4.

Notably, among the seven distinct infiltrating cell clusters, four clusters consisted of CD103+CD8+ T cells, which were previously defined as CD8+ TRMs.14 18 Granzyme B is an important molecule that enables T cells to kill pathogens and cancer cells. We observed that in the CD8+CD103+ T-cell subgroup, the frequency of granzyme B-secreting cells significantly increased (figure 2C). Furthermore, we noted greater Programmed Cell Death protein-1 (PD-1) and Cytotoxic T-Lymphocyte Antigen 4 (CTLA-4) expression in intratumoral CD8+ TRMs than in normal tissues (figure 2D). This observation suggests that this particular cellular population could be a potential target for immune checkpoint blockades (ICB). CD8+ TRMs have been implicated in a variety of clinicopathological features and patient survival outcomes across several solid tumor types.15 26 31 Consistent with these findings, our study revealed a negative correlation between CD8+ TRMs infiltration in tumors and tumor stage as well as tumor/node/metastasis (TNM) stage; conversely, this correlation was not observed for total CD8+ T cells or circulating CD8+ T cells (figure 2E). To further investigate the potential clinical and prognostic implications of CD8+ TRMs infiltration in patients with ESCA, we performed multicolor immunofluorescence staining on paraffin-embedded TMA specimens derived from 77 patients with ESCA. In this study, we employed anti-CD3, anti-CD8A, and anti-CD103 antibodies to visualize and quantify CD8+ TRMs within the TME. The numbers of infiltrating T cells (CD3+), CD8+ T cells (CD3+CD8A+) and CD8+ TRMs (CD3+CD8A+CD103+) were then quantified in different patients. The details of the TMA analysis of the ESCA samples are displayed in online supplemental table 2. We found that the percentages of infiltrating CD8+ TRMs (CD3+CD8+CD103+) and CD3+ T cells were negatively correlated with TNM stage (figure 2F,G, table 1); the infiltration of CD3+ T cells and CD8+ T cells was related to lymph node metastasis; and the infiltration of these cells did not show a clear correlation with histological type, differentiation, age, or sex (table 1). Having established a correlation between cell infiltration and tumor stage, we next sought to investigate whether these infiltrating cell populations are associated with the overall survival of patients. Notably, the frequency of tumor-infiltrating CD8+ TRMs was significantly positively correlated with the overall survival of patients with ESCA (figure 2H), but the frequencies of CD3+ T cells and CD8+ T cells were not significantly correlated (online supplemental figure S2E, F). To corroborate the potential antitumor effects of CD8+ TRMs cells in the context of ESCA, we conducted an in vivo experiment in C57 mice with conditional deletion of CD103+CD8+ T cells. We engrafted the mouse ESCA AKR cell line to construct a subcutaneous xenograft tumor model in CD8creTg/Tg; ITGAEfl/fl C57 mice and CD8creTg/Tg; ITGAE△/△ C57 mice. Notably, the depletion of CD8+ TRMs accelerated tumor growth (figure 2I) and significantly increased tumor burden (figure 2J).

Table 1

Correlation of immune cell infiltration and clinical parameter in tissue microarray

The cytotoxic and exhausted programs tumor-infiltrating CD8+ TRMs

To further confirm the above CyTOF results, especially those for CD8+ T cells, we analyzed scRNA-seq data from 60 patients with ESCA. All the data for each patient were merged, and clustering was performed using Seurat. We generated scRNA-seq libraries for CD3+ T cells, which provided transcriptomes for 54,466 cells. Then, we subjected CD3D+CD8A+ T cells to reclustering and projected them onto a two-dimensional space using UMAP analysis (online supplemental figure S3A). Subclustering of CD8+ T cells revealed three subsets annotated based on published signatures and cell markers.31–34 Differential gene expression analysis of the scRNA-seq data revealed CD8+ T cells subsets with different differentiation statuses and functional states (online supplemental figure S3B, online supplemental table S4). CD8A-C0-GZMK cluster exhibited high expression of cytotoxicity-related transcripts (GZMK, LMNA, GPR183, and FOS), and memory-related transcripts (IL7R, GPR183, and ANXA1), indicating a heterogeneous population of CD8+ Tm cells. CD8A-C1-CXCL13 cluster exhibited expression patterns suggestive of the CD8+ TRMs phenotype, including high expression of CD69, ITGAE, and ITGA1 and low expression of the tissue efflux marker S1PR1(online supplemental figure S3C), similar to CD8+ TRMs described in humans and mice.31 32 35 36 CD8A-C2-CCR7 cluster exhibited a naïve signature, including expression of CCR7 and SELL, and low expression of cytokine and effector genes. To better define cell subtypes, we validated the features of naïvety,27 tissue residency,16 cytotoxicity,26 and exhaustion26 based on previously reported scoring systems (online supplemental figure S3F). Consistent with previous results, CD69, ITGAE (CD103), and ITGA1 (CD49a), especially ITGAE, could be used to define CD8+ TRMs in ESCA (online supplemental figure S3C, D). We also analyzed the pseudotime trajectory of CD8+ T cells using the Monocle2 algorithm to reflect the cell lineage relationships of CD8+ T cells. The results showed that CD8+ TRMs and circulating CD8+ Tm cells were located at the two ends of the pseudotime trajectory (online supplemental figure S3E), indicating a clear demarcation between CD8+ TRMs and other memory CD8+ T cells. Compared with other CD8+ T cells subsets, CD8+ TRMs in ESCA secreted a greater number of cytotoxic molecules (online supplemental figure S3F). In addition, we discovered that CD8+ TRMs displayed a more pronounced exhaustion signature than other subsets (online supplemental figure S3D). As expected, we noted that the expression of immune checkpoints was primarily observed in the CD8+ TRMs population (online supplemental figure S3G). These results indicate that CD8+ TRMs are the central tumor-killing cells and targets of ICBs in ESCA.

Supplemental material

Supplemental material

The tTRMs cluster is the main cluster that mediates antitumor immunity within tumor-infiltrating CD8+ TRMs

Notably, CD8+ TRMs, which possess superior antitumor properties and are clinically relevant in ESCA, did not exhibit significant differences between tumor and normal tissues, similar to circulating CD8+ T cells (figure 3A,B). We speculate that not all CD8+ TRMs in ESCA are tumor-reactive T cells. Since we found no significant differences in total CD8+ TRMs, we performed a detailed analysis of the subset. Using CD49a, CD69, and CD103 expression as markers, we identified five populations of CD8+ TRMs (C19, C20, C25, C26 and C27) in ESCA (figure 3C). C19 and C25 were the most common CD8+ TRMs subgroups in normal tissues and tumor tissues, respectively (figure 3D). To further characterize these CD8+ TRMs clusters in ESCA, we assessed the expression of additional markers including CD45RO, CD39, and CD28 (figure 3C). C25 was nearly six times more abundant in tumor tissues than in normal tissues (figure 3D, online supplemental figure S2B), with high expression of CD39 (figure 3C); the cells were phenotypically characterized as CD39+ CD8+ TRMs, named tumor CD8+ TRMs (tTRMs). In contrast, C19 was three times more common in normal tissues than in tumor tissues (figure 3D, online supplemental figure S2B), and these cells were phenotypically characterized as CD39CD45ROloCD28 CD8+TRMs (figure 3C), named normal CD8+ TRMs (nTRMs). C26 highly expressed CD45RO andCD28 but not CD39 (figure 3C). These cells were characterized as CD28+CD45ROhiCD8+ TRMs and were named pTRMs. C27 exhibited elevated expression of CD28 (figure 3C), a molecule known to play pivotal roles in T cells activation and differentiation. Our analysis revealed that C27 was positively correlated with C19 (nTRMs), C20, C25 (tTRMs), and C26 (pTRMs) (figure 3E). Based on these results, it can be hypothesized that C27 cells have the potential to differentiate into other CD8+ TRMs, and thus being considered as dTRMs, defined as CD28+CD45ROintCD8+ TRMs. To verify the potential differentiation ability of dTRMs, we sorted dTRMs from tumor tissues using flow cytometry; the Transforming Growth Factor-β (TGF-β) signaling pathway, which has been reported to be essential for the formation and maintenance of CD8+ TRMs,37 therefore, we stimulated dTRMs with TGF-β in vitro and detected their phenotypic changes. The results showed that dTRMs in tumors can be induced to differentiate into the CD39+CD8+ TRMs. This finding confirms the potential differentiation ability of dTRMs (figure 3F,G). Studies have shown that CD101 has an inhibitory effect on T-cell proliferation38 and function.39 Our research revealed that CD101 was expressed in C19 and C20 (figure 3H), particularly in C20, suggesting that these clusters contained cells in a relatively suppressed state; therefore, we designated the cells in C20 to be inhibitory CD8+ TRMs (iTRMs). Correlation analysis revealed a significantly greater correlation between dTRMs and nTRMs and a significantly higher correlation between pTRMs and tTRMs in tumors (figure 3E), suggesting that certain factors in the TME change the relationships between CD8+ TRM clusters.

Figure 3

The tTRMs is the main CD8+ TRMs cluster mediating antitumor immunity. (A and B) Box plot showing the percentages of circulating CD8+ T cells (A) and CD8+ TRMs (B) in normal and tumor tissues of the esophagus. The bounds of the boxes represent the 25th and 75th percentiles, and centers represent median. Cyan indicates normal tissue; dark blue indicates tumor tissue. Each dot represents one cell, and each color represents one cell cluster. (C) Heatmap showing the immune molecule expression of distinct T-cell clusters. Red indicates high expression, purple indicates low expression, yellow indicates medium expression. (D) The bar graph represents the average proportion of different CD8+ TRMs in normal tissue (normal) and tumor tissue (tumor). Green represents C19, orange represents C20, pink represents C25, yellow represents C26, light orange represents C27. The number below the bar graph contains the specific average proportion values. (E) Correlations of the five different CD8+ TRMs clusters in normal and tumor tissues. The correlation coefficients of CD8+ TRMs clusters are represented by a coloring scheme from red (positive correlation) to white (negative correlation), while blue represents an absence of correlation. (F) Flow cytometry representative plot of differentiation induction experiments detected the conversion of dTRMs to tTRMs after 72 hours of stimulation with TGF-β (20 ng/mL). (G) The bar graphs show the percentages of CD39+TRMs in the control (black) and TGF-β-stimulated groups (pink) of differentiation induction experiments (n=4). The error bars represent the mean±SD. (H and K) Density map of the expression of CD101 (H) and CD38 (K) in distinct CD8+ TRMs clusters. Blue indicates normal tissue, and orange indicates tumor tissue. (I and J) Scatter plot showing the expression of granzyme B (I) and PD-1 and CTLA-4 (I) in five CD8+ TRMs clusters. Each dot represents one patient. Blue indicates normal tissue, and red indicates tumor tissue. The bars represent the median. (L) Relationships between tTRMs infiltration frequency and clinicopathological parameters in esophageal cancer. The different colors represent different groups; the error bars represent the mean±SD. dTRM, differentiable CD8+ TRMs; iTRM, inhibitory CD8+ TRM; nTRM, normal CD8+ TRMs; TNM, tumor/node/metastasis; TRMs, tissue resident memory T cells; tTRM, tumor CD8+ TRMs; TGF, Transforming Growth Factor;PD-1, Programmed cell death protein 1;CTLA-4, Cytotoxic T-Lymphocyte Antigen 4.

It has been reported that neoantigen-specific tumor-infiltrating lymphocytes (TILs) are almost exclusively CD39+ cells, whereas CD39 TILs represent bystander cells.40 41 Our research found that tTRMs specifically expressed CD39 (figure 3C). In addition, our analyses revealed that both tTRMs and pTRMs exhibited higher levels of the effector molecule granzyme B than did other CD8+ TRMs subsets. Compared with normal tissues, tumor tissues showed significantly increased expression of granzyme B in tTRMs, but not in pTRMs (figure 3I). Furthermore, the proportion of infiltrating tTRMs in tumors was substantially greater than the proportions of other cell clusters (figure 2D). Taken together, our findings suggest that tTRMs is the primary subtype of CD8+ TRMs, play a crucial role in antitumor immunity. We next investigated the expression of immune checkpoint molecules on CD8+ TRMs in tumor tissues compared with matched normal tissues. Notably, tTRMs exhibited the highest levels of immune checkpoint expression, with significant upregulation of PD-1 and CTLA-4 expression compared with other CD8+ TRMs clusters (figure 3J). These results suggest that tTRMs may represent a promising target for ICB therapy for treating ESCA. The therapeutic effect of ICBs varies significantly between individual patients42; hence, the reasons for this phenomenon were preliminarily explored. Several previous studies have reported that the upregulation of CD38 expression can impede the efficacy of antibodies targeting PD-1/(Programmed Cell Death Ligand-1) PD-L1.43 We found widespread higher expression of CD38 in CD8+ TRMs populations in tumors, especially in C25 (figure 3K). This result may be associated with the clinical outcomes of ICBs in patients with ESCA, and combination therapy targeting this molecule may enhance the therapeutic efficacy. Finally, we conducted an analysis to determine the correlation between various subsets of CD8+TRMs and the clinical characteristics of patients with ESCA. Our findings revealed that tTRMs infiltration was inversely correlated with tumor stage (figure 3L). However, no such correlation was observed between the tumor stage and other clusters of CD8+TRMs (online supplemental figure S4A–D).

Supplemental material

Transcriptome and TCR profiles of tumor-infiltrating CD8+ TRMs

To examine the diversity of CD8+ TRMs at the transcriptional level, we isolated triple-positive T cells (CD69+ITGAE+ITGA1+) from CD3D+CD3E+CD8A+ T cells and analyzed their transcriptome profiles (online supplemental figure S5A). Unsupervised clustering of 10,027 CD8+ TRMs by Seurat revealed five clusters (figure 4A), and the gene expression profile of each CD8+ TRMs cluster is shown in the heatmap (online supplemental figure S5B, online supplemental table 5). C0-CD8A-TRM-CXCL13 was the most abundant cluster among the CD8+ TRMs (figure 4B). In C0-CD8A-TRM-CXCL13, there was higher expression of the ITGAE, and the immune checkpoints PD-1 and CTLA-4 (figure 4C), which were defined in the above studies as tTRMs. C1-CD8A-TRM-IL-7R was the second-largest subgroup of CD8+ TRMs (figure 4B) in the tumor; this cluster had relatively lower expression of immune checkpoints than the tTRM (figure 4C, online supplemental figure S5C), corresponding to the nTRMs identified by CyTOF. C3-CD8A-TRMs-STMN1 exhibited higher expression of genes associated with proliferation, such as MKI67 and STMN1 (figure 4C, online supplemental figure S5B), and constituted a very small percentage of CD8+ TRMs (figure 4B), akin to the pTRMs previously described. Next, we scored the proliferative capacity of different CD8+ TRMs and found that C4-CD8A-TRMs-LTB had a lower proliferative capacity (online supplemental figure S5D), consistent with the high CD101 expression of iTRMs (figure 3H). The clusters showed varying resident and exhaustion scores (online supplemental figure S5E), indicating that they may differ in their developmental and functional states. Thus, we used the Monocle2 package to infer the developmental trajectories of various CD8+ TRMs. As shown in figure 4D,E and C2-CD8A-TRM-HSPA6 can evolve into other CD8+ TRMs, corresponding to the aforementioned dTRMs. Moreover, pTRMs were adjacent to tTRMs and displayed a similar state in the pseudotime path analysis (figure 4D), and CyTOF analysis showed that in both tumor and normal tissues, the frequencies of pTRMs had a positive correlation with those of tTRMs (figure 3E). Therefore, we hypothesized that pTRMs may serve as reserve cells for tTRMs and can replenish tTRMs through proliferation and differentiation. To validate this hypothesis, we used flow cytometry to sort pTRMs (CD28+CD45ROhiCD8+TRMs) in tumors; IL-2 can stimulate the proliferation and differentiation of T cells activated by specific antigens; therefore, we cultured the sorted cells under IL-2-containing conditions. Subsequent flow cytometry analysis showed that the cultured pTRMs cells could differentiate into CD39+TRMs under the stimulation with IL-2 (figure 4F,G).

Supplemental material

Supplemental material

Figure 4

Transcriptional and TCR profiles of heterogeneous CD8+ TRMs. (A) UMAP plot of heterogeneous CD8+ TRMs clusters. Each point represents a single cell, colored according to the respective cell cluster: C0: red; C1, pink; C2: green; C3: blue; C4: yellow. (B) The bar graph represents the average proportion of different CD8+ TRMs in tumor tissue. C0: red; C1, pink; C2: green; C3: blue; C4: yellow. (C) Dot plot illustrating the Z-score-normalized mean expression (color scale) and percentage of expressing cells (size scale) for specified genes in each CD8+ TRMs cluster. (D) Monocle2 state plot. Each point in the graph represents a cell, and the color represents the state of the cell. The trace shows the cell state change, and the arrow direction indicates the direction of the cell state change. (E) Pseudotime trajectories for different CD8+ TRMs clusters were constructed using Monocle2. Each dot signifies one cell, and each color denotes one cell cluster. The direction of the arrow indicates the progression of the cell state change. (F) Flow cytometry representative plot of differentiation induction experiments detected the conversion of pTRMs to tTRMs after 72 hours of stimulation with IL-2 (100 ng/mL). Above: control group; below: IL-2 treatment group. (G) The bar graphs show the percentages of CD39+ TRMs cells in the control (black) and IL-2 treatment groups (pink) of differentiation induction experiments (n=4). The error bars represent the mean±SD. (H) The heatmap displays area under the curve score t values for the regulation of expression by transcription factors, computed using the SCENIC workflow. Red indicates enriched gene expression; blue indicates decreased gene expression. (I and J) Bar plots showing the numbers of clonotype (I) and clonal cells (J) in each CD8+ TRMs cluster. Different colors represent different clonotypes. Clonotypes were divided into unique (n=1) and clonal (n=2 and n≥3) categories based on cell number. Clonal cells contained at least two cells of the same clonotype. (K) The pie chart showing the proportion of cloned cells and non-cloned cells in each CD8+ TRMs cluster. Different colors represent different groups. Dark yellow represents cloned cells, blue represents non-cloned cells. (L) The pie chart showing the percentage of TCR clones shared between distinct TRMs clusters. Different colors represent different groups. Orange represents proportion of shared clonotypes, blue represents proportion of unshared clonotypes. Numbers indicate the specific proportion of shared clonotypes. (M) The number of TCR clones shared between distinct CD8+TRMs clusters. The intersection relationship between different clusters is represented by points connected by lines; the bar chart in the lower left represents the size of each cluster clone type; the lower right shows the intersection situation of each cluster, and the bar chart above represents the number of TCR intersections in the included clusters. dTRMs, differentiable CD8+ TRMs; IL, interleukin; iTRM, inhibitory CD8+ TRM; nTRM, normal CD8+ TRMs; TCR, T-cell receptor; TRMs, tissue resident memory T cells; tTRM, tumor CD8+ TRMs; UMAP, Uniform Manifold Approximation and Projection.

To further clarify the regulatory mechanisms of distinct CD8+ TRMs types, the SCENIC (single-cell regulatory network inference and clustering) method was used to construct a transcription factor (TF) regulatory network to study the mechanism of cellular heterogeneity in CD8+ TRMs. The SCENIC analysis revealed six transcriptional regulatory modules consisting of 18 highlighted TFs (figure 4H). To explore the functions of the modules, we extracted the TFs and their regulators from each module and performed Gene Ontology (GO) enrichment analysis (online supplemental figure S6A, online supplemental table 7). Module 1 encompassed major biological processes such as DNA replication and cell division, consistent with the hyperproliferative phenotype of pTRMs (online supplemental figure S6A). Module 2, characterized by differentiation regulation, expressed mainly in dTRMs; the TGF-β signaling pathway is crucial for the formation and maintenance of CD8+ TRMs and is also enriched in this cluster, which indicates that dTRMs are a developing cluster of CD8+ TRMs (online supplemental figure S6B), which is consistent with their differentiable phenotype. Module 2 was also present in pTRMs, consistent with their differentiation ability. Module 3 exhibited high enrichment of the interferon and tumor necrosis factor signaling pathways, indicating effective regulation of the killing potential of the module; the module was significantly more common in tTRMs than in other types (online supplemental figure S6C). Module 6 showed enrichment of pathways related to the inhibition of cell proliferation, which may be associated with the low proliferative capacity of nTRMs and iTRMs and is consistent with high CD101 expression and a low proliferation score. The above results show that the different phenotypes and functional states of CD8+ TRMs are regulated by their intrinsic transcriptional programs. Further research on the roles of different types of CD8+TRMs in immune responses and how to optimize immune responses by regulating these TFs is highly important for optimizing the function of CD8+ TRMs.

Supplemental material

Supplemental material

The TCR repertoire of different CD8+ TRM clusters

Next, we analyzed TCR sequences to explore the lineage structures and relationships of CD8+ TRMs. Following a series of quality control and filtration steps, we identified a total of 8,773 cells in five distinct clusters, all of which exhibited full-length TCR α and β chains s (online supplemental table 8). Our analysis identified 1,034 unique clonotypes, each comprised of a pair of productive α and β chains of TCRs. Notably, 540 of these clonotypes were detected in two or more cells, resulting in a total of 7,740 clonal cells. We observed clonal expansion in each cluster, with clonal sizes ranging from 2 to 665 (online supplemental figure S7A); almost 90% of the CD8+ TRMs were clonal (online supplemental figure S7B), demonstrating the clonal advantage of CD8+ TRMs. However, the TCR clonotype composition of CD8+ TRMs varied widely among patients; some patients showed more than 80% clonal expansion, while some patients showed less than 10% clonal expansion and some even showed no clonal expansion (online supplemental figure S7C). By analyzing the clonal types of different clusters, we found that nTRMs had the richest clone types (figure 4I); and that tTRMs had the greatest number of clonal cells (figure 4J). Further analysis of the clonal degrees of different clusters revealed that tTRMs and dTRMs had the highest clonal proportion (figure 4K). We found that there was a wide range of shared clonal types between different clusters (online supplemental figure S7D); tTRMs and the other four clusters had the most extensively shared clonal types (figure 4L). Among these shared clonal types, 28 TCR sequences coexisted in five clusters (figure 4M). This indicates that CD8+ TRMs with different functional states have a common clonal origin.

Supplemental material

Composition and phenotypes of myeloid cells in ESCA

In the TME of ESCA, infiltrating inflammatory immune cells are known to play a critical role in regulating tumor formation, progression, and immunotherapy response.10 12 44 Thus, we determined the differences in these cell populations between the tumor tissues and matched normal tissues. We employed the PhenoGraph algorithm to cluster the myeloid immune cell (CD45+CD3CD20CD56)based on the expression patterns of immune-related molecules(figure 5B) and subsequently visualized the results in a t-SNE plot (figure 5A). The group was phenotypically subdivided into neutrophils, eosinophils, iMos, classical monocytes, double-negative (DN) DCs, type 1 conventional dendritic cells (cDC1s), cDC2s, other non-classical myeloid cells and undefined cells (figure 5C). We then compared and investigated the composition of myeloid cell subsets in tumor and adjacent tissues. Classic monocytes (cMos) were the most common myeloid cluster in tumor and paracancerous tissues, while iMos was the smallest cluster (figure 5D). We observed a significant increase in the frequency of eosinophils and iMos in tumor tissues compared with paracancerous tissues (figure 5E). In contrast, the abundance of cDC2s was reduced in tumors, whereas the numbers of neutrophils, cMos, DN-DCs, and cDC1s were not significantly different between tumor and paracancerous tissues (figure 5E). We then examined interpatient variability in the infiltrating myeloid cell populations (figure 5F). Notable differences were found in the infiltration of myeloid cell subsets between different patients. For example, in most tumor tissues, the predominant infiltrating myeloid cells were cMos, but in some patients, the largest infiltrating population was granulocytes (figure 5F).

Figure 5

Composition, immunophenotype and characteristics of myeloid cells in esophageal cancer. (A) Cytometry by time of flight t-distributed stochastic neighbor embedding plot of myeloid cells showing the compositions of the 14 primary clusters. Each point represents a single cell, colored according to the respective cluster. (B) Cell markers used to characterize myeloid cell phenotypes. Ag: antigen. (C) The heatmap delineates the expression of immune molecules in distinct myeloid immune cell clusters. Red indicates high expression, yellow indicates moderate expression, and blue indicates low expression. (D) The bar graph represents the average proportion of myeloid cell subgroups in normal tissues (normal) and tumor tissues (tumor). The color corresponds to the cluster type. The three-line table below the bar graph contains the specific average proportion values. (E) Comparison of the frequencies of various myeloid cell populations between normal and tumor tissues. Each point represents a sample (n=20). The error bars represent the mean±SD. (*p<0.05, **p<0.01, ***p<0.001, ****p<0.0001, ns (no significance)). (F) Proportions of myeloid cell subsets in normal and tumor tissues of each patient. The colors indicate different subsets. Each bar represents a patient. (G) Single-cell RNA sequencing UMAP plot of myeloid cells showing the composition of nine primary clusters. Different colors represent different myeloid subsets. Each point represents a single cell. (H) Dot plot of selected marker genes in distinct myeloid cell clusters. Yellow represents high expression, and purple represents low expression. The color scale indicates the average expression level. The size scale indicates the percentage of expressing cells. DC, dendritic cells; NK, natural killer; t-SNE, t-distributed stochastic neighbor embedding; UMAP, Uniform Manifold Approximation and Projection; HLA-DR, Human Leukocyte Antigen-DR; PD-L1, Programmed Cell Death Ligand-1.

We next investigated the intrinsic changes in these cell populations by scRNA-seq. Unsupervised clustering was used to assess the populations of infiltrating myeloid cells, which were sorted into nine distinct groups (figure 5G). The differentially expressed transcripts among the nine clusters are shown in online supplemental table 9. The heatmap displays the molecular signatures of specific subpopulations of myeloid cells (online supplemental figure S8A). We identified nine cell populations based on their gene expression profiles, which included cMos, iMos,45 cDC1s,46 47 cDC2s,11 plasmacytoid dendritic cells (pDCs),48 Mast cells,49 and Macro (macrophage)-1, Macro-2, and Macro-3 populations (figure 5H, online supplemental figure S8A). Considering the significant difference in the infiltration of cDC2s and iMos in the ESCA CyTOF results, we performed GO analysis to identify key biological functions of the two myeloid subsets. The main pathways affected by iMos were associated with receptor-mediated endocytosis, viral entry into host cells, the inflammatory response, and proteolysis (online supplemental figure S8B), suggesting that these cells may be related to the processing of protein antigens in ESCA. Peptide antigen assembly with the major histocompatibility complex (MHC) class ll (MHC-II) protein complex, positive regulation of T cells activation, positive regulation of T cells-mediated cytotoxicity and T cells proliferation were enriched in cDC2s, demonstrating the vital role of cDC2s in the control of T cells function in ESCA (online supplemental figure S8C).

Supplemental material

Supplemental material

The iMo/cDC2/CD8+ TRMs axis sustains the generation and activation of CD8+ TRMs in the TME

Studies have demonstrated that myeloid cells play a critical role in the development and maintenance of CD8+ TRMs,50 however, the specific interactions between myeloid cells and CD8+ TRMs require further investigation. We used CellChat to analyze the intercellular communication among myeloid cells and CD8+ TRMs clusters. There were numerous shared receptor-ligand interactions between immune cells, highlighting the close communication between these immune cell subsets (figure 6A). Considering the observed differences in cDC2 infiltration between normal and tumor tissues, the GO analysis of cDC2s revealed enrichment of pathways related to T cells activation, cytotoxicity, and proliferation pathways. Additionally, network analysis of cell-to-cell communication showed that cDC2s have regulatory effects on various CD8+ TRMs (figure 6B); thus, our focus was on examining the regulation of CD8+ TRMs by cDC2s. To test this hypothesis and determine whether the interaction between cDC2s and CD8+ TRMs is direct or indirect, we performed Immunofluorescence (IF) staining, and our IF analysis of ESCA tissues confirmed a direct interaction between cDC2 and CD8+ TRMs (figure 6C). Furthermore, to explore the effect of cDC2s on CD8+ TRMs, we separated cDC2s and CD8+ TRMs from ESCA tissues by flow cytometry-based sorting and found that co-culture of the two subsets significantly increased the percentage of Ki67-positive CD8+ TRMs (figure 6D, online supplemental figure S9A), indicating the important function of cDC2s in promoting the proliferation of CD8+ TRMs. Typically, antigens are processed and presented to CD4+ T cells via the MHC-II pathway and to CD8+ T cells via the MHC class I (MHC-I) pathway. To further explore the regulation of CD8+ TRMs by cDC2s, we examined the expression of MHC-I and MHC-II in cDC2s in ESCA. Interestingly, we observed MHCI+ cDC2s in both tumor and normal tissues, although the frequency was significantly lower in tumors (figure 6E, online supplemental figure S9B). In contrast, we did not observe a significant difference in the expression of Human leukocyte antigen-DR (HLA-DR, an MHC class II molecule) between tumor and normal tissues (online supplemental figure S9C). In addition, MHC-I expression is an important feature of DC maturation,51 52 and low infiltration of MHC1+ cDC2s in tumors suggests possible arrest of cDC2s development in tumors. The costimulatory molecule CD80 is regarded as a maturation marker for DCs52 53; therefore, we next examined the levels of infiltrating CD80+ cDC2s in tumors, and consistent with our speculation, the level of infiltrating CD80+cDC2s in tumors was decreased (figure 6F, online supplemental figure S9D), indicating that the development of cDC2s in the TME was restricted. To verify whether the function of cDC2 in tumors is impaired, we sorted cDC2 in paracancerous tissues and tumor tissues throughflow cytometry and co-cultured them with CD8+ TRMs. This study revealed that compared with the ability of cDC2s from paracancerous tissues, the ability of cDC2s from tumors to promote the proliferation of CD8+ TRMs cells is significantly decreased (figure 6G, online supplemental figure S9E) This observation implies that while cDC2s may have the potential to directly interact with CD8+ TRMs, this ability appears to be compromised within the TME.

Supplemental material

Figure 6

The iMo/cDC2/CD8+ TRMs axis sustains the generation and activation of CD8+ TRMs in the tumor microenvironment. (A) Interaction plot showing the number and strength of predicted intercellular interactions between myeloid cells and distinct CD8+ TRMs clusters. The thickness of the lines represents the number and strength of the interactions. Thick lines represent stronger interactions. (B) Detailed view of the strengths of interactions between cDC2s and other immune cell subgroups. The thickness of the lines represents the strength. (C) Immunofluorescence staining of cDC2 and CD8+ TRMs in paraffin sections of ESCA tissues. Blue: DAPI; red: CD8; orange: CD103; green: CD11c; pink: CD1c. (D) Flow cytometry was used to detect the proliferation of CD8+ TRMs. After the flow sorting of CD8+ TRMs (CD45+CD3+CD8a+CD69+CD103+) in the tumor, typically divided into two groups. One group was cultured alone, and the other group is co-cultured with cDC2s (CD45+CD3CD56CD20CD11C+CD1C+) that are sorted from the tumor at a ratio of 1:5 (n=5 per group). Flow cytometry detection is performed 48 hours later. Different colors represent different groups. The black represents cultured alone; the pink represents the co-cultured group. The error bars represent the mean±SD. P values are shown above the bar graph. (E and F) Flow cytometry analysis was used to identify the frequencies of MHC-I+ cDC2 (E) and CD80+ cDC2 (F) collected from ESCA tumor tissue and tumor-adjacent tissue (n=5). The black color represents tumor-adjacent tissue; the pink represents tumor tissue. The error bars represent the mean±SD. P values are shown above the bar graph. (G) Flow cytometry detects the proliferation properties of CD8+ TRMs (n=4). After the flow sorting of CD8+TRMs in the tumor, divided into two groups. One group is co-cultured with cDC2 sorted from paracancerous tissues, and the other group is co-cultured with cDC2 sorted from tumor tissues. cDC2 and CD8+TRMs are cultured at a ratio of 1:5, and flow cytometry detection is performed after 48 hours of culture. The black represents co-cultured with cDC2 sorted from paracancerous tissues; the pink represents the co-cultured with cDC2 sorted from tumor tissues. The error bars represent the mean±SD. P values are shown above the bar graph. (H) Direct comparison of cDC2 marker (CD11c and CD1c) and iMo marker (CD14 and CD16) expression in iMo between normal and tumor tissues. Each point represents a sample; tumor samples are indicated in red, and normal samples are indicated in green. Box plot: whiskers represent min to max. Bounds of boxes represent the 25th and 75th percentiles, and centers represent median. (I) Immunofluorescence detection of GM-CSF and IL-4 expression in ESCA tissues. Above: GM-CSF staining; below: IL-4 staining. (J) Differentiation induction experiments detected the conversion of iMos into cDC2s (n=5). The iMos in peripheral blood were sorted by flow cytometry and divided into two groups. The black represents cultured alone; the pink represents the cultured with IL-4 (10 ng/mL) and GM-CSF (20 ng/mL) for 48 hours. The error bars represent the mean±SD. P values are shown above the bar graph. cDC, conventional dendritic cells; cMos, classic monocytes; ESCA, esophageal cancer; IL, interleukin; iMos, intermediate monocytes; MHC, major histocompatibility complex; TRM, tissue resident memory T cells; DAPI, 4’,6-diamidino-2-phenylindole; pDC, Plasmacytoid dendritic cells; GM-CSF, Granulocyte-Macrophage Colony-Stimulating Factor.

iMos are of interest because, relative to other monocyte subsets, these cells exhibit proinflammatory properties and are associated with the risk of inflammatory disease.54 55 CyTOF analysis revealed increased iMo infiltration (figure 5E) and CD66b expression of iMos (figure 5B) in tumor tissues compared with normal tissues. Previous studies have suggested that CD66b+ monocytes are a proinflammatory myeloid cell subset in cancer,56 and the enrichment of this subset in ESCA is consistent with the reports that ESCA is an inflammatory tumor.5 6 Through CyTOF analysis, we found that iMos in tumors simultaneously expressed higher levels of the cDC2 markers CD11c and CD1c (figures 5C and 6H). The monocyte lineage is known to be highly adaptable to changes within the TME and can differentiate into various functional subsets of macrophages or DCs in response to different signals and stimuli.57 Therefore, we hypothesized that iMos have the ability to differentiate into cDC2s. Monocytes usually differentiate into DCs cells in the presence of IL-4 and GM-CSF58; therefore, we detected the expression of IL-4 and GM-CSF in ESCA and found that these two cytokines are present in ESCA (figure 6I). Next, we conducted an experimental study, in which we separated iMos from peripheral blood throughflow cytometry and induced their differentiation to verify whether iMo can differentiate into cDC2s. After 2 days of culture with recombinant human GM-CSF and IL-4, iMos differentiated into cDC2s (figure 6J, online supplemental figure S9F), suggesting that iMos are one of the sources of cDC2s in ESCA.

Hypoxia disrupted the iMo/cDC2/ CD8+ TRMs axis in ESCA

To further elucidate the factors influencing ESCA TIME formation, we performed a bulk transcriptomic analysis of patient with ESCA samples. WGCNA was used to identify modules highly correlated with ESCA development (figure 7A). Interestingly, we found that the turquoise module had a pronounced effect on tumor development (figure 7B). GO enrichment analysis revealed enrichment of genes related to the hypoxia and the immune response in the turquoise module (figure 7C). This observation implies that the ESCA TIME may be closely related to hypoxia. Both experimental research and clinical research have revealed that hypoxia changes the immune status of the TME,59 60 but the mechanism by which increased hypoxia modulates tumor immunity remains unclear. Given these results, we next studied the effect of hypoxia on the iMo/cDC2/CD8+ TRMs axis. As shown in figure 6, we found that cDC2 maturation and function were restricted to tumors. This raises the question of whether this phenomenon is caused by hypoxia. To explore the potential reasons for this observation, a single-cell suspension derived from tumors was incubated under hypoxic conditions. Flow cytometry results revealed that the number of CD80+ cDC2s CD80+ cDC2s decreased after incubation in a hypoxic environment compared with after incubation in a normoxic environment (figure 7D, online supplemental figure S10A), indicating that a hypoxic environment inhibits cDC2 maturation. Since hypoxia inhibits the development of cDC2, whether hypoxia affects the transformation of iMos into cDC2s was tested. Cytokine-induced transformation of iMos into cDC2s was significantly inhibited after culture under hypoxic conditions (figure 7E, online supplemental figure S10B). To explore the possible reasons why hypoxia inhibiting the differentiation of iMos into cDC2s, we sorted iMos in peripheral blood and detected their apoptosis and differentiation into macrophages under hypoxic and normoxic conditions. The research results show that the hypoxic environment promotes the apoptosis of iMo (figure 7F,G) but does not change its ability to differentiate into macrophages (figure 7H,I). These results demonstrated that the hypoxic TME of ESCA may prevent CD8+ TRMs from exerting an antitumor immune response by inhibiting the generation and maturation of cDC2s.

Supplemental material

Figure 7

Hypoxia disrupts the iMo/cDC2/CD8+ TRM axis in ESCA. (A) Cluster dendrogram of ESCA samples with bulk transcriptome sequencing data. WGCNA, weighted gene coexpression network analysis. (B) Module-trait-relationship plot. The x-axis represents phenotypic traits, while the y-axis represents corresponding modules, with different colors indicating different modules. The square indicates the correlation and p value between the phenotypic traits of the x-axis and each module of the y-axis. The legend on the right indicates the correlation. As the value approaches 1, the color deepens to red, indicating a positive correlation. As the value approaches −1, the color deepens to blue, indicating a negative correlation. (C) Enriched Gene Ontology biological process terms of the MEturquoise module. The x-axis represents fold enrichment. The y-axis represents the cellular biological process. The size of the point represents the number of genes. The color of the point represents the significance of enriched genes. The depth of color reflects the degree of significance. (D) Determination of CD80+ cDC2 frequency by flow cytometry under normoxic and hypoxic conditions (n=5). A hypoxic environment was induced by treating with 100 µmol/L of cobalt chloride (CoCl2) for 48 hours. Flow cytometry was used to detect the changes in the single-cell suspension with or without CoCl2.The black color represents the group without CoCl2 treatment; the pink color represents the CoCl2 treatment group. The error bars represent the mean±SD. P values are shown above the bar graph. (E) Detection of cDC2s by flow cytometry after IL-4 and GM-CSF stimulation under normoxic and hypoxic conditions (n=6). Single-cell suspension is detected by flow cytometry after 48 hours treatment of CoCl2. Different colors represent different groups. The black color represents the group without CoCl2, IL-4, and GM-CSF treatment; the pink color represents the IL-4 and GM-CSF treatment group; the cyan color represents the CoCl2, IL-4, and GM-CSF treatment group. The error bars represent the mean±SD. P values are shown above the bar graph. (F) Representative plot of annexin V-FITC/PI double staining flow cytometry. The above is control group (without CoCl2 treatment group) and the below is CoCl2 (100 µmol/L, 48 hours) treatment group. (G) Statistical analysis of iMos apoptosis rate (including early and late apoptotic cells) under normoxic and hypoxic conditions (n=3). Different colors represent different groups. The black color represents without CoCl2 treatment group; the pink color represents CoCl2 treatment group. The error bars represent the mean±SD. P values are shown above the bar graph. (H) Differentiation induction experiments detected the conversion of iMos into macrophage by flow cytometry. The iMos in peripheral blood were sorted by flow cytometry and divided into two groups. The control group represents cultured with IL-4 (10 ng/mL) and GM-CSF (20 ng/mL) for 48 hours; the hypoxia group represents the cultured with IL-4 (10 ng/mL), GM-CSF (20 ng/mL) and CoCl2 (100 µmol/L) for 48 hours. (I) Statistical bar plot of iMos transformed into macrophages by flow cytometry (n=5). The black represents cultured with IL-4 (10 ng/mL) and GM-CSF (20 ng/mL) for 48 hours; the pink represents the cultured with IL-4, GM-CSF and CoCl2 (100 µmol/L) for 48 hours.The error bars represent the mean±SD. P values are shown above the bar graph. (J) Schemata of the hypoxia /iMo/cDC2/CD8+TRMs immune axis in ESCA. cDC, conventional dendritic cells; ESCA, esophageal cancer; IL, interleukin; iMos, intermediate monocytes; MHC, major histocompatibility complex; TRM, tissue resident memory T cells; FSC-H, Forward Scatter Height; GM-CSF, Granulocyte-Macrophage Colony-Stimulating Factor; PI, Propidium Iodide.

In summary, in a normoxic environment, iMos transform into cDC2s under the cytokines stimulation and then activates CD8+ TRMs to kill target cell; in contrast, under hypoxic conditions, the apoptosis of iMos increases, leading to a decrease in their transformation to cDC2s. At the same time, the maturation of cDC2s is inhibited, which prevents the full activation of CD8+ TRMs, resulting in CD8+ TRMs being unable to exert effective tumor-killing effects (figure 7J).

Discussion

The present study first comprehensively characterized the immune landscape of the ESCA by CyTOF analysis, revealing that the TIME differed from a normal cellular environment. Second, we established a relationship between CD8+ TRMs and the clinicopathological characteristics of patients with ESCA, elucidated the heterogeneity of CD8+ TRMs, and identified the relationship between the ontogeny and function of different CD8+ TRM clusters. Third, we revealed the evolutionary relationship between iMos and cDC2s and the interaction between cDC2s and CD8+ TRMs. Finally, we preliminarily explored the inhibitory effects of a hypoxic environment on the iMo/cDC2/CD8+ TRMs axis.

Previous studies have characterized the ESCA TIME using scRNA-seq.21 61 62 However, these studies considered only tens of thousands of cells and usually had large sampling biases caused by differences between different patients, and the small clusters identified were prone to false positives or false negatives.63 More importantly, scRNA-seq defines cell populations at the transcriptional level, and studies have shown that there is no significant correlation between the abundance of RNA and protein expression levels.64 For instance, our study showed that T cells highly expressed the immune checkpoint molecule PD-1 in ESCA, but in these scRNA-seq studies, the expression of PD-1 was lower than that of other immune checkpoints. Therefore, establishing correlations between genotypes and phenotypes from transcriptomic data alone remains a challenge.65 With its ability to accurately and sensitively analyze large numbers of cells in a short period, CyTOF enables high-throughput functional and phenotypic analysis of single-cell suspensions.66 In this study, we integrated scRNA-seq and CyTOF to investigate the TIME and provide more precise data on tumor-infiltrating immune cell composition, lineage, and functional status. We applied CyTOF to profile the landscape of immune cells in ESCA tumors and paracancerous tissues and discovered differences in the infiltration of multiple immune subsets in the tumor. However, although T cells were the largest group of infiltrating immune cells in ESCA, they did not show significant differences in frequency, which is inconsistent with previously published scRNA-seq results.21 61 62 This may be due to the incomplete correlation between protein expression and messenger RNA expression. The infiltration of various immune cells and the large proportion of T cells suggest that ESCA is not an “immune desert” tumor, and the application of immunotherapy in ESCA is likely to benefit patients. Our study examined the expression of multiple immune checkpoints in T cells, and revealed that PD-1 and CTLA-4 were the most highly expressed. This result is consistent with the encouraging antitumor activity of PD-1 blockade in ESCA.67 68

Our research also revealed substantial differences between CD8+ TRMs and other CD8+ Tm cells.36 Moreover, we discovered that CD8+ TRMs infiltration was linked to the ESCA tumor stage and proved to be a better predictor of patient survival than T cells or CD8+ T cells infiltration. Studies have identified various subsets of human CD8+ TRMs in multiple tissues,36 prompting us to investigate whether all CD8+ TRMs possess tumor-killing capabilities. Compared with other CD8+ TRMs, tTRMs with the highest immune checkpoint expression were the most abundant in tumors and exhibited a clear tumor-killing function. These findings are consistent with a recent study on CD8+ TRMs heterogeneity, which identified the TIM-3+IL-7R TRMs subset as the tumor-specific CD8+ TRMs population with the highest PD-1 expression in lung cancer.69 Additionally, we found a cluster of highly proliferative CD8+ TRMs, which had previously been described in breast cancer,36 that displayed a significant positive correlation with tTRMs. The scRNA-seq analysis and differentiation induction assay revealed that the cluster possessed robust proliferation and differentiation capabilities and served as a reserve cell pool for tTRMs. The CD8+ TRMs subset with the lower PD-1 expression represented the primary population in normal esophageal mucosa (nTRM), corresponding to the cells in the C1-TRM-IL-7R cluster identified through scRNA-seq data analysis. This cluster demonstrated limited exhaustion signatures and proliferative capacity, consistent with recently reported graft CD8+ TRMs lacking exhaustion characteristics.17 TCR clonotypes also play roles in determining the developmental trajectories of different clusters.33 61 We found that nTRMs had the most abundant clonotype, which may be related to the recognition of various antigens in the esophagus. The tTRMs contained the largest number of clonal cells, which may be related to the substantial expansion of CD8+ TRMs targeting specific tumor antigens. In addition, our analysis revealed the same TCR clonotypes among various clusters, indicating the common clonal origin of these clusters. These striking results support the existence of functional and developmental heterogeneous clusters within CD8+ TRMs.

In our study, we observed substantial infiltration of DCs and monocytes in tumors, with cDC2s and iMos showing differential infiltration. Previous studies have suggested that cDC2s present antigens to CD4+ T cells in an MHC-II-dependent manner. Recent studies have shown that DCs activate naïve CD8+ T cells and efficiently generate CD8+ TRMs following exposure to foreign antigens in an MHC-I-dependent manner.70 Currently, we have preliminarily shown that cDC2s can interact with and enhance the proliferation of CD8+ TRMs, which mean that cDC2 is vital for maintaining the generation of CD8+ TRMs. To our knowledge, this is the first report showing that cDC2s play a positive regulatory role in the generation of CD8+ TRMs. iMos are another subset of myeloid cells that show differential infiltration in ESCA. It is now generally accepted that cMos are capable of differentiating into monocyte-derived macrophages and monocyte-derived DCs.71 Nevertheless, the development of iMo is poorly understood. Interestingly, our study showed that iMos had the propensity to differentiate into cDC2s on stimulation with GM-CSF and IL-4. Our findings serve as a good supplement to the existing research on iMo function.

Hypoxia is one of the central players in shaping the TIME, and studies have shown positive feedback between hypoxia and inflammation signaling.72 Chronic inflammation leads to local hypoxia, which exacerbates inflammation.73 However, the complex interplay of immune cells infiltrating in the hypoxic TME still needs to be elucidated. Previous studies have shown that lasting deleterious effects on DC maturation after hypoxic challenges, such as reduced MHC-II expression by cDC2s under hypoxia, enable cDC2s to present antigens to regulatory T cells.74 Similarly, our data indicate that cDC2 maturation and MHC-I expression are limited under hypoxic conditions. Nonetheless, the specific mechanism underlying this phenomenon requires further investigation. While numerous studies have shown that exposure to hypoxia can enhance monocyte recruitment75; there is limited research on the effect of hypoxia on monocyte differentiation, particularly with respect to iMo. Our research revealed that iMos have the potential to differentiate into cDC2s, but under hypoxic conditions, the apoptosis rate of iMos increases, and their differentiation into cDC2s is inhibited. However, the underlying mechanisms remain to be elucidated.

Our research provides a theoretical basis for understanding the ESCA immune microenvironment, but there are still some limitations. First, due to technical limitations, our depiction of the immune landscape is not comprehensive, for example, the depiction of immune groups such as pDC and macrophage subgroups by CyTOF still needs to be improved. Second, many of the findings deserve further exploration. For example, the specific mechanism by which tTRMs exert antitumor effects and inhibit cDC2s maturation under hypoxic conditions, both need further research and explanation. Finally, additional in vivo and in vitro experiments are needed to verify the CyTOF and single-cell transcriptome sequencing results.

Code availability

CyTOF analysis follows the general program:

PhenoGraph: https://github.com/jacoblevine/PhenoGraph;

Xshift: https://github.com/ginberg/xshift_operator;

PARC:https://github.com/ShobiStassen/PARC;

Example scripts of single cell sequencing data for processing and analyzing data are available at (https://github.com/junglingfish/EC). Any additional information may be reasonably obtained from the respective author.

Data availability statement

Data are available in a public, open access repository. Data are available upon reasonable request. All data relevant to the study are included in the article or uploaded as supplementary information. The transcriptome sequence data of ESCA were obtained from the Genome Sequence Bank of Beijing Institute of Genomics, Chinese Academy of Sciences (https://ngdc.cncb.ac.cn/search/all?&q=HRA004010). The single-cell RNA sequencing data of ESCA were downloaded from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE160269).

Ethics statements

Patient consent for publication

Ethics approval

This study involves human participants and was approved by the Second Affiliated Hospital of Zhejiang University School of Medicine, IH2020001026, Zhejiang Cancer Hospital MR-33-21-014018. Participants gave informed consent to participate in the study before taking part.

Acknowledgments

We would like to thank Professor Xiannian Zhang of Beijing Advanced Innovation Center, School of Basic Medicine, Capital Medical University, for generously sharing the sequencing data of esophageal samples.

References

Supplementary materials

Footnotes

  • CW, HY, FL and XH contributed equally.

  • Contributors MW and PW designed the study. CW, HY, FL and XH are responsible for sample processing, data analysis and flow cytometry experiments. ZL and YL are jointly accountable for sample collection and help coordinate this study’s aspects. ZW and QW provided advice on the study. MW, PW, CW, HY and FL wrote this paper with comments from all authors. CW, HY, FL and XH contributed to the manuscript equally.The first author and the corresponding author will serve as the guarantors for the overall content of our paper. They accept full responsibility for the work and the conduct of the study, have access to the data.

  • Funding We would like to thank Zhejiang Province Key Technology R&D Program (No. 2023C03064) and the National Natural Science Foundation of China (No. 81272594) for their support of this experiment.

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