Immune characteristics and genetic markers of esophageal cancer by single-cell analysis: Implications for immunotherapy

doi:10.21203/rs.3.rs-1536299/v1

Download PDF

Research Article

Immune characteristics and genetic markers of esophageal cancer by single-cell analysis: Implications for immunotherapy

https://doi.org/10.21203/rs.3.rs-1536299/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Esophageal cancer (EC) is one of the most common cancers worldwide. The prognosis for patients with the same stage of EC can vary substantially. Progress of single-cell analysis technology has furthered the understanding of tumor heterogeneity. This paper aimed to apply single-cell analysis to explore the characteristics of the tumor environment of EC and provide inspiration for personalized treatment.

Methods: The latest gene expression data and clinical follow-up information of single-cell sequencing results of EC samples were downloaded using from The Cancer Genome Atlas (TCGA) Genomic Data Commons (GDC) Application Programming Interface (API). A differential gene function analysis of the immune infiltration signature agents in the tumor microenvironment was undertaken using bioinformatics analytic methods, and potential molecular targets were sought.

Results: We identified specific cell subsets in the EC and paracancerous samples, including panel cells, natural killer (NK) cells, exhausted CD8+ T cells, CD8+ memory T (Tcm) cells, and effector memory T (Tem) cells, including B cell enrichment in the cancer samples. Differences were detected between B cells and monocytes in Stage II and Stage III tumors, and these differences may have related to RNA transcription and degradation. The CXCL8 protein was identified as a valid potential prognostic marker.

Conclusions: Cell groups with homogenous cell surface markers exhibit intercellular variations that exert a considerable effect on cell function. Our study will contribute to the understanding of the tumor immune microenvironment (TME) and cellular heterogeneity in EC patients and serve as a valuable resource for in-depth exploration of the pathogenesis of EC and the identification of potential therapeutic targets in the future.

esophageal cancer

genetics

immune

single cell analysis

Esophageal cancer (EC) is the 7th most common malignancy in the world, with more than 500,000 new cases diagnosed annually [1]. It is the 6th leading cause of cancer death, causing more than 500,000 cancer deaths globally each year, and the 5-year overall survival (OS) rate for EC patients is less than 20% [2]. In recent years, the complexity of the processes involved in cancer formation has been recognized. Even patients with the same tumor pathological type and stage of may have different prognoses; it is necessary to decipher the complexity of the tumor and the tumor microenvironment. In recent years, the awareness of immunotherapy has gradually increased. Immunotherapy has been proven advantageous in advanced or metastatic EC in both the first-line and second-line therapeutic setting [3–5], and has been established as standard care[6].

The tumor microenvironment (TME) consists of malignant cells and immune cells, stromal cells, blood cells, endothelial cells, and fat cells[7]. Immune and stromal cells are essential components of the TME and have proven valuable for the survival of cancer patients[8]. For example, checkpoint-blocking immunotherapy can benefit a limited subset of cancer patients for whom cancer interacts with immunity[9, 10]. The failure of several anti-angiogenic cancer therapies is thought to be due to metabolic adaptation and reprogramming of cancer cells and abnormalities of endothelial cells and their interaction with pericyte cells[11].

Single-cell transcriptomics has developed rapidly since 2009, when the first single-cell transcriptome profile was described. Cell populations with homogenous cell surface markers harbor cell-to-cell variations, with a considerable impact on cell function. This heterogeneity can be resolved using single-cell RNA sequencing (scRNA-seq)[12, 13]. EC transcriptomes have been based mostly on traditional RNA sequencing techniques in previous studies. The traditional RNA sequencing technology has led to deeper study of the transcription group characteristics of EC and a deeper understanding of the origin and development of EC. However, studies based on traditional RNA sequencing have focused on the overall level of the sample, and more research is needed to understand the heterogeneity of a cluster of cells and individual cells. Single-cell sequencing can provide clues to the development of disease at the individual cellular level. Researchers used scRNA-seq to dissect the complex immune responses of peripheral blood mononuclear cells in acute Kawasaki disease[14]. In addition, single-cell sequencing analysis has also been applied to a variety of cancers, such as multiple myeloma, lung adenocarcinoma, and stomach adenocarcinoma.

In this study, we downloaded public data from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases to perform scRNA-seq on large-scale cells derived from ESCC tumor samples. We discovered some interesting facts that may help deepen our understanding of the immune microenvironment and immunotherapy in EC and provide some new inspiration for EC patients.

Data download and preprocessing.

We used TCGA Genomic Data Commons (GDC) Application Programming Interface (API) to download the latest gene expression data and clinical follow-up information for EC on 2 February, 2022. The GSE145370 dataset was downloaded from the GEO database as GSE145370_RAW.tar. In processing TCGA data, only normal and primary tumor samples were retained. For GSE145370 data processing, we used R (version 3.6.4, https://www.r-project.org/) and the Seurat R software package (version 3.6, https://satijalab.org/seurat/) to combine single-cell sequencing results for each sample from the GSE145370_RAW dataset to build a single-cell dataset.

Quality control of single-cell data.

As above, we used R and the Seurat R package for the analyses. The filter criteria of cells were determined after referencing previous studies[15, 16]. According to the median number of genes and the percentage of mitochondrial genes in the samples, cells with < 200 and > 4,000 genes (potential cell duplets) and a mitochondrial gene percentage of < 20% were filtered. After quality control (QC), 103,223 high-quality cells were obtained. The relationship between the percentage of mitochondrial genes and the messenger RNA (mRNA) reads was detected and visualized as the relationship between the number of mRNAs and the reads of mRNA. After data normalization, all highly variable genes in single cells were identified after controlling for the relationship between average expression and dispersion. Subsequently, we used principal component analysis (PCA) with variable genes as the input and identified significant principal components (PCs) based on the JackStraw function. A total of 20 PCs were selected as the input for uniform manifold approximation and projection (UMAP) and t-distributed stochastic neighbor embedding (t-SNE) when statistically significant. With a resolution of 0.5, cells were clustered by the FindClusters function and classified into 24 different cell types (clusters 0–23). Next, we used the FindAllMarkers function to find differentially expressed genes (DEGs) between each type of cell.

Cell cycle bias analysis.

Cell cycle analysis was performed using the Seurat program (CellCycleScoring function). We used a previously defined core set of 43 G1/S and 54 G2/M cell cycle genes. Cells were classified according to the maximal average expression (‘cycle score’) in these 2 gene sets. When the cycle scores of G1/S and G2/M were both less than 2, we considered these cells noncycling; otherwise, we considered the cells proliferative. After the cell cycle analysis, no bias induced by cell cycle genes was observed.

Notes for cell subgroups.

Cell type assignment was performed based on the marker genes reported in previous studies[17], the CellMaker database, and the SingleR package. Subclustering analysis of B cells, monocytes, effector CD8 + memory T (Tem) cells, CD8 + T cells, CD4 + T cells, and naive CD8 + T cells was performed using the SubsetData function in the Seurat program. Following data normalization and scaling, all highly variable genes in single cells were identified. Subsequently, 20 PCs were selected as the input for t-SNE when statistically significant. The FindClusters function was used to classify cells into different clusters. Next, the FindAllMarkers function was used to find DEGs between each subtype. Next, to assign pathway activity estimates to individual cells, the R package gene set variation analysis (GSVA) with the standard settings was implemented in the single-sample gene set enrichment analysis (ssGSEA), and a gene set containing all human pathways was downloaded from the Molecular Signatures Database (MsigDB), including the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. We identified significantly enriched pathways defined by the Benjamini–Hochberg-corrected P-value < 0.05.

Feature-rich analysis.

The KEGG pathway enrichment analysis systems were applied to the R package, clusterProfiler[18], to identify over-represented KEGG pathway genes. In addition, we used the Signal Pathway Impact Analysis (SPIA) software package to activate and inhibit the KEGG pathway based on the different genes, and used GSEA to analyze the KEGG pathway enrichment of cells and samples from different sources, using the gene expression spectrum. For both analyses, a q-value of < 0.05 denoted statistical significance.

Single-cell data set quality control.

We downloaded the GEO dataset GSE145370 and analyzed it using the Seurat R package. The number of minimal information units (MIUs) was not significantly correlated with mitochondrial gene percentage (Fig. 1A), but positively correlated with the number of mRNAs (Fig. 1B). We also analyzed the number of mRNAs, the readings of mRNA, and the distribution of mitochondrial genes (Fig. 1C), which revealed that most genes were distributed between 0 and 2,500 bp, mRNA readings were distributed between 0 and 1,000 bp, and mitochondrial percentages were below 15%. Further, we removed cells with more than 20% expression originating from mitochondrial genes and those with fewer than 200 genes. The expression spectrum of 33,538 genes from 103,223 cells was obtained, and the number of cells in each sample was counted (Table 1), including the number of filtered mRNAs, the readings of mRNA, and the distribution of mitochondrial genes (Fig. 1D). Highly variable genes (top 2,000) were identified for downstream analysis (Fig. 1E) after quality control.

Identification of cell subgroups by t-SNE analysis.

Considering that cells of the same type may be classified into different clusters when clustering because of different cell cycles, we conducted cell cycle regression analysis to observe such effects. We found that in each region, 3 cycles of cells were randomly distributed, and there was no difference in cell cycle status, so there was no need to eliminate the effects of cell cycles (Fig. 2A). Further, we used main component analysis to extract cell characteristics and the distance matrix, to identify 24 cell groups (Fig. 2B). There were different lymph node (LN) ratios in these cell groups, of which clusters 0, 1, 6, and 4 were mainly from tumor tissue, and clusters 5, 8, 2, and 3 were mainly from tissue peripheral to cancer (Fig. 2C). Further, applying rank and testing methods to identify differences in genes using the cell type annotation based on the CellMarker software package (diff_genes_wilcox_fdr.txt per cell group), these 24 clusters were eventually annotated to 16 cells (all.marker.anno .txt), which we defined as: Cell Cluster 0, cell monocytes from tumor tissue; Cell Cluster 1, Paneth cells; Cell Cluster 2, B cells; Cell Cluster 3, monocytes 1; Cell Cluster 4, exhausted CD8 plus T cells; Cell Cluster 5, CD8 + plus effector memory T (Tem) cells; Cell Cluster 6, natural killer (NK) cells; Cell Cluster 8, monocytes 2. In addition, we observed different proportions of tumor samples in different monocyte populations, with B cells derived mainly from normal samples and Paneth cells mainly from tumor samples (Fig. 2D), with significant individual differences observed in these cells in different patients.

Analysis of specific cell subgroups and functions in cancer and cancer side samples.

We observed the concentration of Paneth cells, NK cells, and exhausted CD8 + plus T cells in the tumor samples, and the concentration of CD8 + effector memory T (Tem) cells in the cancer side samples, which suggested the presence of specific cell subgroups in the cancer and cancer-side samples. We analyzed the differences in gene expression in each cell subgroup of cancer and cancer side samples, and we observed the differences in gene expression. Among the Panel cell, Natural killer cell, exhausted CD8 + T cell, the Natural killer cell has the most differential genes, and most of the difference genes in the Exhausted CD8 + T cell are contained by the Natural killer cell As shown in Fig. 3A. We selected the genes common to the three types of cells to include in the R software package clusterProfiler for KEGG Functional enrichment analysis which used a total of 29 pathways (step3_A_KEGG_enrich.txt), with the most significant top 10 shown in Fig. 3B. In addition to our analysis of the genetic differences between the CD8 + effector memory T (Tem) cell, B cell cancer, and cancer side sample tissues, we observed 2 types of genetic differences with a large intersection (Fig. 3C). Similarly, we identified common differences in genes using the R software package, clusterProfiler, for the KEGG functional enrichment analysis, which used a total of 22 pathways (step3_B_KEGG_enrich.txt), with the top 10 shown in Fig. 3D. We observed a significant overlap between the different gene functions in cells in the cancer and cancer side sample tissues (Fig. 3E). For example, genes with differences in the specific cell subgroups of cancer and cancer side samples were highly concentrated in the Ribosome and Parkinson's disease pathways, which showed that there were similarities of dysfunction in the specific cell subgroups of cancer and the cancer side samples.

Cell subgroups and disorder pathways for different clinically phased patients.

To identify prognostically related cell subgroups, we compared the differences between Stage II and Stage III patients across cell types. We observed significant differences in cell sources in all cell groups (step4.stage.count.txt), suggesting significant differences in tumor microenvironmental characteristics between samples with different clinical stages, the most significant of which were among the B cell and monocyte groups. The B cell count was significantly higher in Stage II than in Stage III, and the monocyte count was significantly higher in Stage II than in Stage III, which suggested unique immune cell immersion in the different stages. We compared the differences in gene expression between cancer and cancer side samples of the B cell and the monocyte, and observed that most of the degenerated genes in the monocyte intersected with 324 genes in the B cell (Fig. 4A). Further, using the R software package SNP Panel Identification Assay (SPIA), we analyzed differences in the pathways of gene participation between the 2 types of cells. In B cells, we observed a total of 40 pathways with disorders (e.g., step4_spia. A .txt), of which 24 were activation pathways (Fig. 4B) and 16 pathways were inhibition pathways (Fig. 4C). In monocytes, we observed a total of 42 offset pathways (e.g., step4_spia. B .txt), of which 25 were pathway activations (Fig. 4D) and 17 were inhibition pathways (Fig. 4E). Among them, 14 of the activated pathways are consistent, and 10 of the suppressed pathways are consistent, which indicates that the signaling pathways involved have some consistency overall. However, there are 6 pathways to the contrary, 5 of which are activated in tumor samples in B cells but inhibited in monocytes, including protein processing in the endoplasmic reticulum, RNA transport, influenza A, the mitogen-activated protein kinase (MAPK) signaling pathway, toxoplasmosis, and RNA degradation, mainly related to RNA transcription and degradation.

Transcription of classical mononucleosis.

We also observed the presence of 5 different mononucleocytes in the monocyte1 marker-based cell subgroups due to CD302, SLC43A2, BACH1, QKI, monocyte1-based RIN3, ODF3B, H CK, LAT2, FGR, SFT2D1, IQGAP1, and CLEC7A, the monocyte2 marker-based cell groups due to IL1B, S100A12, SERPINB2, VCAN, EREG, F3, THBS1, PTGS2, GCA, PROK2, GK, EHD1, SLC25A37, and JMJD1C, the monocyte3 marker-based cell groups due to G0S2, VIM, SH3 BGRL3, FBP1, CTSB, DUSP1, IL1RN, LGALS3, and S100A11, the monocyte4 marker-based cell groups due to IFI30, LST1, RNF144B, ACSL1, FPR2, C5AR1, CD55, CUX1, CD300E, C19orf38, IRAK3, and. NLRP3. These monocytes exhibited different cancer sample ratios, with high proportions of tumor samples in the Monoocyte3 group and high proportions of arranged samples in the Monoocyte1, Monocyte2, and Monocyte3 groups, which showed different mononucleosis components in the cancer and cancer side samples. Further, in our subgroup clustering of these monocytes, we obtained 6 cell subgroups (Fig. 5A), which included monocytes, CD8 + T cells, and microglial cells. We observed that microglial cells mainly originated from tumor samples, and Monocyte2 originated mainly from cancer side samples (Fig. 5B). Further comparison of the pathway enrichment levels between these cell clusters, such as by step5.spia .txt using the SPIA method showed that CD8 + T cell abnormalities had the least pathway impacts, and may represent a class of mononucleocytes that make up a small number of normal samples, in addition to several other classes that showed a large number of different path abnormalities (Fig. 5C), We observed multiple inhibition pathways in microglial cells with significant activation pathways in other cells, such as RNA transport and T cell receptor signaling pathways, which suggested the presence of a unique class of dysfunctional monocyte clusters in the tumor samples. In addition, clusters of microglial cells, CD8 + T cells, monocyte2 cells, and monocyte3 cells represented a higher level of major tissue compatibility complexes (MHCs).

mRNA expression and differences between B-cell subclasses.

The B lymphocytic dysfunction, an important component of adaptive immunity, is thought to be important in the pathogenesis of tumors, where 2 B cell groups are detected. There were 11,104 cells in the B cell group, labeled as CD69, FCER2, BANK1, FCMR, ADAM28, SIDT1, CXCR5, LY9, NGLY1, PARP15, BLK, and ARHGAP24. Of the cells in the B cell group, 1,312 were marker-based groups due to LRMP, TCL1A, BCAS4, CD22, VPREB3, SWAP70, and EZR. These B cells exhibited different proportions of tumor samples, and through t-SNE analysis, the B cells were divided into 6 clusters (Fig. 6A), where the CD8 + T cells, NK cells, and microglial cells were mainly located in the tumor samples and the microglial cells and neutrophils were similarly located in the tumor and cancer side samples. This suggested that there was a similar group of cells (microglial cell, neutrophil) in the cancer and cancer side samples. To compare the functional differences between the tumor samples and the cancer side samples in the microglial cells and neutrophils, we calculated the difference in the KEGG pathway enrichment of the tumors in the 2 cell groups using the GSEA method. We observed that a total of 106 pathways were enriched in microglial cells. These pathways were enriched in the cancer side samples (e.g., step6.gsea.kegg.txt), the most significant of which were CARDIAC_MUSCLE_CONTRACTION, OOCYTE_MEIOSIS, and HUNTINGTONS_DISEASE (Fig. 6B). Interestingly, the neutrophils had a total of 55 pathways (step6.gsea.kegg_ B.txt) which were concentrated in the cancer side samples, the most significant of which were the first 3 pathways ALZHEIMERS_DISEASE, SELENOAMINO_ACID_METABOLISM, and PARKINSONS_DISEASE (Fig. 6C), all associated with geriatric diseases. In addition, there was a significant intersection between the enriched pathways in the neutrophil and the rich pathways of microglial cells (Fig. 6D). This suggested that there were some similarities in the disordered pathways of the tumor samples at the pathway level in the different types of cells.

Identification of the genetic signature associated with prognosis.

To identify prognostically related genetic markers, we selected the cell group with the closest proportion of normal and tumor samples. We extracted the gene marker from the source cell subgroup, obtaining a total of 3,154 gene markers. First, we calculated that the 4,211 gene expression differences in Stage III, Stage II, and normal samples contained a total of 460 genes, (scRNAsub.exp.0 edgp .txt). We defined the sample feature vector as follows: normal 0, Stage II 2, and Stage III 3. We calculated the Spearman’s rank correlation coefficient of the expression of these genes in each cell and sample feature vector, selecting a false discovery rate (FDR) < 0.05 to obtain a total of 177 significantly related genes (scRNAsub.exp.edgp.cr .txt). In the functional enrichment analysis of these genes using the R software package, clusterProfiler, we observed that they were mainly concentrated in the KEGG pathways, such as the ribosome, leukocyte transendotal migration, and viral myocarditis pathways (Fig. 7A). In addition, they enriched the molecular functions of the structural constituents of the ribosome (Fig. 7B), and we observed that these genes were enriched into a large number of biological processes and cellular components, such as nuclear-transcribed mRNA catabolic processes and nonsense-mediated decay (Fig. 7C). In cellular components, such as cell-substrate-adhesin junctions and focal adhesions (Fig. 7D), these biological pathways are closely related to tumor proliferation and metastasis. We mapped these genes to the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) version 11.0b (www.string-db.org) database to obtain protein interoperability networks (ppi.network.txt) of these genes, and analyze the interrelationships of these genes to the network degree distribution (ppi.attr.table.txt) (Fig. 7E) and the Closeness centrality distribution (Fig. 7F). The Betweenness centrality distribution curve (Fig. 7E). The proximity to centrality, and the greater the amount of centrality indicate the importance of the network. We selected the genes of the top 5 degrees, Closeness centrality, and Betweenness centrality, and then selected the intersection of the three to identify two important genes, CXCL8 and CCL2.Finally, we observed that the CXCL8 gene had higher expressions in multiple cell groups (Fig. 7H) and that the CCL2 gene expressed only weak signals in each cell population; hence, we considered CXCL8 an effective potential prognostic marker.

TCGA dataset validation.

To verify the prognostic correlation of CXCL8, we downloaded the TCGA-Expression Spectrum Data Set (TCGA-ESDS; TCGA-ESCA.exp.txt) and the corresponding clinical follow-up data (TCGA-ESCA.exp.os .txt), from which the expression spectrum of CXCL8 was extracted, and differences in expression between the cancer and cancer side tissue samples (Fig. 8A) could be observed. In the tumor samples, the expression of CXCL8 was significantly higher than in the cancer side samples. In addition, we compared the relationship between CXCL8 expression and prognosis, using the maxstat package of R, to group CXCL8 expression, and analyzed the prognosis by considering differences of overall survival (OS), disease-free survival (DFS), and progression-free survival (PFS) in patients with high expression and low expression. We observed significant differences in prognosis and the low expression group, which indicated that CXCL8 is a potential prognostic marker in EC.

Many patients are diagnosed with advanced EC at first diagnosis, having already surpassed the optimal treatment period, and thus receiving a relatively poor prognosis[19]. Currently, the most common method of predicting the prognosis for patients with EC is the American Joint Committee on Cancer (AJCC) staging system[20]. However, patients at the same clinical stage may still experience different survival periods, which may be related to different molecular characteristics of the tumor. Therefore, it is important to develop a more sensitive prognostic diagnostic approach based on the molecular characteristics of the patient's tumor by identifying new prognostic markers. We identified cell subgroups using t-SNE analysis and found differences in cell cluster enrichment between tumor tissue and normal tissue, where clusters 0, 1, 6, and 4 were mainly found in tumor tissue and clusters 5, 8, 2, and 3 were mainly found in paracancerous tissue. In addition, we observed different proportions of cell populations emerging in the tumor samples. These results indicate significant individual differences in these cells in different patients.

Tumor heterogeneity has long been one of the important obstacles in the clinical treatment of tumors. Different patients may receive different treatment responses due to complex immunogenicity[21, 22]. Existing immunotherapy is based on immune checkpoint-specific treatment options for T cell antigen receptors; however, not all patients benefit from immunotherapy due to the immune heterogeneity of tumors[23]. Most patients show significant differences in response and resistance, suggesting that further analysis of individual tissues and individual cells in tumor patients is necessary to develop specific immunotherapy strategies[24]. Based on recent reports and sequencing data, immune cells can more systematically describe the molecular signatures of the immune system[25–27].

Genes that differ from the specific cell subgroups of Paneth cell, NK cell, exhausted CD8 + T cell, effector CD8 + effector memory cell, and B cell are highly concentrated in the Ribosome and Parkinson's disease pathways. These results showed that there were some similarities in dysfunction between cancer and paracancerous samples in specific cell subgroups. There were significantly different tumor microenvironmental characteristics in samples of different clinical phases, the most significant differences being in the phases of B cells and monocytes. The pathways of gene enrichment in these 2 types of cells were generally the same, and the differences in the pathways were mainly related to RNA transcription and degradation.

Patients with advanced EC, whether untreated or after multi-treatment, may benefit from this treatment[3, 5, 28]. However, the use of immune checkpoint inhibitors (ICIs) in the neoadjuvant treatment phase is still very limited. In the context of tumor immunotherapy, neoadjuvant immunotherapy may have a more prominent advantage as antigen exposure will enhance the degree of the tumor-specific T cell response and prolong the duration of action[4]. In this study, it was found that there were significantly different tumor microenvironmental characteristics in samples with different clinical stages. The B cells and monocytes exhibited different characteristics in Stage II and Stage III, and the MAPK signaling pathway exhibited activation and inhibition due to different gene participation in the 2 types of cells, respectively. Preclinical studies have shown that modulating the MAPK signaling pathway can improve the sensitivity of melanoma to immunotherapy. Inhibiting MEK activity in melanoma cell lines can increase antigen levels, thereby enhancing the immunity of antitumor T cells[29]. In addition, studies have shown that regulating MEK activity can affect the infiltration of immune cells in the tumor microenvironment, including T cells and B regulatory cells, enhancing the response to anti-PD1 immunotherapy, thereby enhancing anti-tumor immunity[30]. In animal models, it has also been shown that regulating MAPK pathway activity can affect the number of inter-tumor, effector-phenotype CD8 T cells, thereby inhibiting tumor growth. After the combined use of pathway regulators with PD-L1 drugs, a synergistic effect has been shown[31]. Another study found that metformin activated the jak1/2/3/STAT5 and AKT/mTOR pathways in a p38 MAPK-dependent manner, significantly enhancing NK cytotoxicity, which increased the infiltration of NK cells and T cells into tumors[32]. By amalgamating the results of this study and previous studies, it is reasonable to assume that the MAPK signaling pathway has a regulatory effect on immune infiltration of the tumor microenvironment and response to immunotherapy. The pre-trial results of Nivolumab combined with neoadjuvant chemoradiotherapy included 16 patients with stage II-III EC or esophageal-gastric junction cancer, 5 patients with pathological complete response (pCR), 9 patients with pathological hypothesis, and 15 patients with R0 resection. The immune microenvironment in samples of different stages and the functional pathways enriched by different gene concentrations can affect immune efficacy to some extent, suggesting that the application of immunotherapy in the early stage of the clinical staging of cancer patients may be reasonable.

The NK cells have the dual function of direct killing and immune regulation, taking on important tasks in complex immune processes, either by lysing tumor cells directly through antibody-dependent cell-mediated cytotoxicity and removing major histocompatibility complex damaged cells in the primary immune response or by programmed death ligand-1 (PD-L1) induction of tumors and recruitment of dendritic cells and subsequent T cells to initiate the tumor microenvironment for secondary adaptive immune response[33]. In this study, we found that NK cells had the most differential genes among the 3 relatively enriched cell populations in paracarcine signaling samples. It can be assumed that differences in NK cell distribution and gene expression in EC samples are closely related to the efficacy of immunotherapy, and this finding also provides information for subsequent studies to enhance the suicidal effect of NK cells and induce secondary adaptive immune responses.

We selected the marker gene in the secretory cell group that was closest to the normal sample and tumor sample ratio, and found that the CXCL8 gene had a higher expression in multiple cell groups. The CXCL8 gene (also known as interleukin 8) is a versatile, pro-inflammatory chemokine that was originally recognized as a neutrophil chemical attractant, and its important role in the development of tumorigenesis has also been discovered in recent years. In several human cancers, CXCL8 has been found to be upregulated at the tumor invasion front. Studies have shown that the interaction of CXCL8 between tumors and the tumor microenvironment is achieved by influencing angiogenesis, genetic diversity, survival, proliferation, immune escape, metastasis, and multidrug resistance, thereby promoting tumor progression[34].

The role of CXCL8 in gastrointestinal tumors is increasingly being demonstrated. There is evidence that exosome release triggers the production of CXCL8, leading to increased neutrophil recruitment and neutrophil extracellular trap formation in the colon cancer tumor microenvironment, ultimately leading to deterioration of colon cancer [35]. Studies have also confirmed that high expression of interleukin (IL)-8 is associated with poor clinical prognosis in colorectal cancer patients[36]. In addition, serum CXCL8 levels are considered an important risk factor for the development of gastric cancer [37]. The CXCL8 gene plays a role in the association of epigenetic changes with the emergence of immunosuppressive microenvironments[38]. Therefore, deciphering the regulatory/signaling cascade of CXCL8 and its downstream effects may have prognostic and clinical prospects for cancer therapy drugs directed by the tumor microenvironment.

Although significant progress has been made in the development of scRNA-seq technology, there are still some limitations. The scRNA-seq of eukaryotic cells is limited to the detection of large amounts of non-lipoprotein RNA (poly(A)-RNA) expression in mammalian cells with poly(A) tail (poly(A)-RNAs). The advantages of reading rRNA sequencing is avoided by initiating oligonucleotides (dT). However, this approach inevitably excludes other RNA information that does not have poly(A), such as circRNA and microRNAs. In recent years, many studies have shown that non-coding RNAs (such as circRNA and microRNAs) play an important role in the development of various tumors, which has become a hotspot for tumor mechanism research. Therefore, ongoing developments in this field are important for further applications of scRNA-seq in potential clinical significance[39, 40].

In this study, we found specific cell subgroups in cancerous and paracancerous samples, such as the Paneth cell, NK cell, exhausted CD8 + T cell in the tumor sample and the CD8 + effector memory T (Tem) cell, and B cell in the normal tissues. Unlike the monocyte in Stage II and Stage III, these differences may be related to RNA transcription and degradation. We believe that CXCL8 is a valid potential prognostic marker. Several clinical samples and prospective studies are still needed to evaluate the value of the risk model for predicting the prognosis and immune response of EC patients and to determine the optimal cutoff value.

Cell groups with homogenous cell surface markers have intercellular variations that exert a considerable effect on cell function. Our study will contribute to the comprehension of the TME and cellular heterogeneity in EC patients and serve as a valuable resource for in-depth exploration of the pathogenesis of EC and identification of potential therapeutic targets in the future.

Esophageal cancer

natural killer

TME

tumor immune microenvironment

ICB

immune checkpoint blockades

TCGA༚The Cancer Genome Atlas༛GEO༚Gene-Expression Omnibus༛GSVA༚Gene set variation analysis༛DEGs

differentially expressed genes

MAPK༚mitogen-activated protein kinase༛OS༚overall survival༛PFS༚progression-free survival༛DFS༚disease-free survival༛PD-L1༚programmed death ligand-1༛

Ethics approval and consent to participate

Not applicable.

Consent for publication

All authors agreed on the manuscript.

Availability of data and material

The data that support the findings of this study are available on request from the corresponding author.

Competing interests

The authors declare that there are no conflict of interests.

Funding

This work was supported by National Natural Science Foundation of China (82102827, 81872454), Chinese Society of Clinical Oncology (Y-Young2020-0003); Beijing Bethune Charitable Foundation (flzh202119).

Authors' contributions

BX and HCZ were involved in the conception and design of the study. HC and I Hsuan Tseng developed the methodology. KLZ and XLG were involved in data acquisition. HC and DSA analyzed and interpreted the data. The manuscript was written by BX and XMF, and it was critically reviewed by all other authors. Administrative and technical support was provided by HCZ, KLZ and XLG. The overall study was supervised by JYS and QL. All authors have read and approved the final manuscript.

Acknowledgements

Not applicable.

Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A and Bray F (2021) Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 71:209–249. doi: 10.3322/caac.21660
Huang J, Koulaouzidis A, Marlicz W, Lok V, Chu C, Ngai CH, Zhang L, Chen P, Wang S, Yuan J, Lao XQ, Tse SLA, Xu W, Zheng ZJ, Xie SH and Wong MCS (2021) Global Burden, Risk Factors, and Trends of Esophageal Cancer: An Analysis of Cancer Registries from 48 Countries. Cancers (Basel) 13. doi: 10.3390/cancers13010141
Kato K, Cho BC, Takahashi M, Okada M, Lin CY, Chin K, Kadowaki S, Ahn MJ, Hamamoto Y, Doki Y, Yen CC, Kubota Y, Kim SB, Hsu CH, Holtved E, Xynos I, Kodani M and Kitagawa Y (2019) Nivolumab versus chemotherapy in patients with advanced oesophageal squamous cell carcinoma refractory or intolerant to previous chemotherapy (ATTRACTION-3): a multicentre, randomised, open-label, phase 3 trial. Lancet Oncol 20:1506–1517. doi: 10.1016/S1470-2045(19)30626-6
Huang J, Xu J, Chen Y, Zhuang W, Zhang Y, Chen Z, Chen J, Zhang H, Niu Z, Fan Q, Lin L, Gu K, Liu Y, Ba Y, Miao Z, Jiang X, Zeng M, Chen J, Fu Z, Gan L, Wang J, Zhan X, Liu T, Li Z, Shen L, Shu Y, Zhang T, Yang Q, Zou J and Group ES (2020) Camrelizumab versus investigator's choice of chemotherapy as second-line therapy for advanced or metastatic oesophageal squamous cell carcinoma (ESCORT): a multicentre, randomised, open-label, phase 3 study. Lancet Oncol 21:832–842. doi: 10.1016/S1470-2045(20)30110-8
Kojima T, Shah MA, Muro K, Francois E, Adenis A, Hsu CH, Doi T, Moriwaki T, Kim SB, Lee SH, Bennouna J, Kato K, Shen L, Enzinger P, Qin SK, Ferreira P, Chen J, Girotto G, de la Fouchardiere C, Senellart H, Al-Rajabi R, Lordick F, Wang R, Suryawanshi S, Bhagia P, Kang SP, Metges JP and Investigators K- (2020) Randomized Phase III KEYNOTE-181 Study of Pembrolizumab Versus Chemotherapy in Advanced Esophageal Cancer. J Clin Oncol 38:4138–4148. doi: 10.1200/JCO.20.01888
Shah MA, Kennedy EB, Catenacci DV, Deighton DC, Goodman KA, Malhotra NK, Willett C, Stiles B, Sharma P, Tang L, Wijnhoven BPL and Hofstetter WL (2020) Treatment of Locally Advanced Esophageal Carcinoma: ASCO Guideline. J Clin Oncol 38:2677–2694. doi: 10.1200/JCO.20.00866
Hui L and Chen Y (2015) Tumor microenvironment: Sanctuary of the devil. Cancer Lett 368:7–13. doi: 10.1016/j.canlet.2015.07.039
Hirata E and Sahai E (2017) Tumor Microenvironment and Differential Responses to Therapy. Cold Spring Harb Perspect Med 7. doi: 10.1101/cshperspect.a026781
Sharma P, Hu-Lieskovan S, Wargo JA and Ribas A (2017) Primary, Adaptive, and Acquired Resistance to Cancer Immunotherapy. Cell 168:707–723. doi: 10.1016/j.cell.2017.01.017
Su S, Zhao J, Xing Y, Zhang X, Liu J, Ouyang Q, Chen J, Su F, Liu Q and Song E (2018) Immune Checkpoint Inhibition Overcomes ADCP-Induced Immunosuppression by Macrophages. Cell 175:442–457 e23. doi: 10.1016/j.cell.2018.09.007
Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, Berent-Maoz B, Pang J, Chmielowski B, Cherry G, Seja E, Lomeli S, Kong X, Kelley MC, Sosman JA, Johnson DB, Ribas A and Lo RS (2016) Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell 165:35–44. doi: 10.1016/j.cell.2016.02.065
Stein CM, Weiskirchen R, Damm F and Strzelecka PM (2021) Single-cell omics: Overview, analysis, and application in biomedical science. J Cell Biochem 122:1571–1578. doi: 10.1002/jcb.30134
Potter SS (2018) Single-cell RNA sequencing for the study of development, physiology and disease. Nat Rev Nephrol 14:479–492. doi: 10.1038/s41581-018-0021-7
Wang Z, Xie L, Ding G, Song S, Chen L, Li G, Xia M, Han D, Zheng Y, Liu J, Xiao T, Zhang H, Huang Y, Li Y and Huang M (2021) Single-cell RNA sequencing of peripheral blood mononuclear cells from acute Kawasaki disease patients. Nat Commun 12:5444. doi: 10.1038/s41467-021-25771-5
Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, Baglaenko Y, Brenner M, Loh PR and Raychaudhuri S (2019) Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods 16:1289–1296. doi: 10.1038/s41592-019-0619-0
Young MD, Mitchell TJ, Vieira Braga FA, Tran MGB, Stewart BJ, Ferdinand JR, Collord G, Botting RA, Popescu DM, Loudon KW, Vento-Tormo R, Stephenson E, Cagan A, Farndon SJ, Del Castillo Velasco-Herrera M, Guzzo C, Richoz N, Mamanova L, Aho T, Armitage JN, Riddick ACP, Mushtaq I, Farrell S, Rampling D, Nicholson J, Filby A, Burge J, Lisgo S, Maxwell PH, Lindsay S, Warren AY, Stewart GD, Sebire N, Coleman N, Haniffa M, Teichmann SA, Clatworthy M and Behjati S (2018) Single-cell transcriptomes from human kidneys reveal the cellular identity of renal tumors. Science 361:594–599. doi: 10.1126/science.aat1699
Liao J, Yu Z, Chen Y, Bao M, Zou C, Zhang H, Liu D, Li T, Zhang Q, Li J, Cheng J and Mo Z (2020) Single-cell RNA sequencing of human kidney. Sci Data 7:4. doi: 10.1038/s41597-019-0351-8
Yu G, Wang LG, Han Y and He QY (2012) clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16:284–7. doi: 10.1089/omi.2011.0118
Ko HH, Cheng SL, Lee JJ, Chen HM, Wang CW, Cheng SJ and Kok SH (2016) Factors influencing the incidence and prognosis of second primary tumors in patients with oral squamous cell carcinoma. Head Neck 38:1459–66. doi: 10.1002/hed.24457
Hsu PK, Chen HS, Liu CC and Wu SC (2018) Application of the Eighth AJCC TNM Staging System in Patients With Esophageal Squamous Cell Carcinoma. Ann Thorac Surg 105:1516–1522. doi: 10.1016/j.athoracsur.2017.12.032
Seth S, Li CY, Ho IL, Corti D, Loponte S, Sapio L, Del Poggetto E, Yen EY, Robinson FS, Peoples M, Karpinets T, Deem AK, Kumar T, Song X, Jiang S, Kang Y, Fleming J, Kim M, Zhang J, Maitra A, Heffernan TP, Giuliani V, Genovese G, Futreal A, Draetta GF, Carugo A and Viale A (2019) Pre-existing Functional Heterogeneity of Tumorigenic Compartment as the Origin of Chemoresistance in Pancreatic Tumors. Cell Rep 26:1518–1532 e9. doi: 10.1016/j.celrep.2019.01.048
Shibata M and Hoque MO (2019) Targeting Cancer Stem Cells: A Strategy for Effective Eradication of Cancer. Cancers (Basel) 11. doi: 10.3390/cancers11050732
Smith A, Roy A, Karapetis CS, Broadbridge V and Price T (2022) Immunotherapy use in oesophagogastric cancers-a review of the literature. Br J Cancer. doi: 10.1038/s41416-022-01751-4
Servetto A, Salomone F, Di Costanzo F, Iuliano R, Marandino L, Napolitano F, Santaniello A, De Placido P, De Placido S, Di Maio M, Formisano L and Bianco R (2022) Inadequate health-related quality of life assessment and reporting in phase III clinical trials of immune checkpoint inhibitors in solid cancers: a systematic review. Crit Rev Oncol Hematol:103649. doi: 10.1016/j.critrevonc.2022.103649
Li C, Liu F, Sun L, Liu Z and Zeng Y (2022) Natural killer cell-related gene signature predicts malignancy of glioma and the survival of patients. BMC Cancer 22:230. doi: 10.1186/s12885-022-09230-y
Liao P, Wang W, Wang W, Kryczek I, Li X, Bian Y, Sell A, Wei S, Grove S, Johnson JK, Kennedy PD, Gijon M, Shah YM and Zou W (2022) CD8(+) T cells and fatty acids orchestrate tumor ferroptosis and immunity via ACSL4. Cancer Cell. doi: 10.1016/j.ccell.2022.02.003
Masuda K, Kornberg A, Miller J, Lin S, Suek N, Botella T, Secener KA, Baccarella AM, Cheng L, Ingham M, Rosario V, Al-Mazrou AM, Lee-Kong SA, Kiran RP, Stoeckius M, Smibert P, Del Portillo A, Oberstein PE, Sims PA, Yan KS and Han A (2022) Multiplexed single-cell analysis reveals prognostic and non-prognostic T cell types in human colorectal cancer. JCI Insight. doi: 10.1172/jci.insight.154646
Owen D and Chaft JE (2018) Immunotherapy in surgically resectable non-small cell lung cancer. J Thorac Dis 10:S404-S411. doi: 10.21037/jtd.2017.12.93
Kono M, Dunn IS, Durda PJ, Butera D, Rose LB, Haggerty TJ, Benson EM and Kurnick JT (2006) Role of the mitogen-activated protein kinase signaling pathway in the regulation of human melanocytic antigen expression. Mol Cancer Res 4:779–92. doi: 10.1158/1541-7786.MCR-06-0077
Yarchoan M, Mohan AA, Dennison L, Vithayathil T, Ruggieri A, Lesinski GB, Armstrong TD, Azad NS and Jaffee EM (2019) MEK inhibition suppresses B regulatory cells and augments anti-tumor immunity. PLoS One 14:e0224600. doi: 10.1371/journal.pone.0224600
Ebert PJR, Cheung J, Yang Y, McNamara E, Hong R, Moskalenko M, Gould SE, Maecker H, Irving BA, Kim JM, Belvin M and Mellman I (2016) MAP Kinase Inhibition Promotes T Cell and Anti-tumor Activity in Combination with PD-L1 Checkpoint Blockade. Immunity 44:609–621. doi: 10.1016/j.immuni.2016.01.024
Xia W, Qi X, Li M, Wu Y, Sun L, Fan X, Yuan Y and Li J (2021) Metformin promotes anticancer activity of NK cells in a p38 MAPK dependent manner. Oncoimmunology 10:1995999. doi: 10.1080/2162402X.2021.1995999
Guillerey C, Huntington ND and Smyth MJ (2016) Targeting natural killer cells in cancer immunotherapy. Nat Immunol 17:1025–36. doi: 10.1038/ni.3518
Asokan S and Bandapalli OR (2021) CXCL8 Signaling in the Tumor Microenvironment. Adv Exp Med Biol 1302:25–39. doi: 10.1007/978-3-030-62658-7_3
Shang A, Gu C, Zhou C, Yang Y, Chen C, Zeng B, Wu J, Lu W, Wang W, Sun Z and Li D (2020) Exosomal KRAS mutation promotes the formation of tumor-associated neutrophil extracellular traps and causes deterioration of colorectal cancer by inducing IL-8 expression. Cell Commun Signal 18:52. doi: 10.1186/s12964-020-0517-1
Wang Y, Zhao M, He S, Luo Y, Zhao Y, Cheng J, Gong Y, Xie J, Wang Y, Hu B, Tian L, Liu X, Li C and Huang Q (2019) Necroptosis regulates tumor repopulation after radiotherapy via RIP1/RIP3/MLKL/JNK/IL8 pathway. J Exp Clin Cancer Res 38:461. doi: 10.1186/s13046-019-1423-5
Pawluczuk E, Lukaszewicz-Zajac M, Gryko M, Kulczynska-Przybik A and Mroczko B (2021) Serum CXCL8 and Its Specific Receptor (CXCR2) in Gastric Cancer. Cancers (Basel) 13. doi: 10.3390/cancers13205186
Dong F, Qin X, Wang B, Li Q, Hu J, Cheng X, Guo D, Cheng F, Fang C, Tan Y, Yan H, He Y, Sun X, Yuan Y, Liu H, Li T, Zhao Y, Kang C and Wu X (2021) ALKBH5 Facilitates Hypoxia-Induced Paraspeckle Assembly and IL8 Secretion to Generate an Immunosuppressive Tumor Microenvironment. Cancer Res 81:5876–5888. doi: 10.1158/0008-5472.CAN-21-1456
Kharchenko PV (2021) The triumphs and limitations of computational methods for scRNA-seq. Nat Methods 18:723–732. doi: 10.1038/s41592-021-01171-x
Gonzalez-Silva L, Quevedo L and Varela I (2020) Tumor Functional Heterogeneity Unraveled by scRNA-seq Technologies. Trends Cancer 6:13–19. doi: 10.1016/j.trecan.2019.11.010

Table 1

Number of cells in each sample

Tag	beforeCellCount	afterCellCount	percent
GSM4317409	12342	11179	0.91
GSM4317410	8338	7025	0.84
GSM4317411	11018	8759	0.79
GSM4317412	10727	9955	0.93
GSM4317413	15820	14400	0.91
GSM4317414	4986	3944	0.79
GSM4317415	7325	6718	0.92
GSM4317416	6134	5241	0.85
GSM4317417	5057	4797	0.95
GSM4317418	10680	10200	0.96
GSM4317419	4375	4001	0.91
GSM4317420	3502	3039	0.87
GSM4317421	6597	6151	0.93
GSM4317422	8256	7814	0.95

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Immune characteristics and genetic markers of esophageal cancer by single-cell analysis: Implications for immunotherapy

Status:

Version 1

Abstract

Figures

Introduction

Methods

Results

Discussion

Conclusion

Abbreviations

Declarations

References

Tables

Additional Declarations

Status:

Version 1