Shengyi Zhang ¹ Xinyi Zhang ² Xiahou Zhikai ³ Shunqing Zuo ¹ Jialong Xue ¹ Yi Zhang ^1*

• ¹ Department of Thoracic Surgery, Songjiang Hospital Affiliated to Shanghai Jiao Tong University School of Medicine, Shanghai, China
• ² Clinical Medical College, Southwest Medical University, LuZhou, China
• ³ China Institute of Sport and Health Science, Beijing Sport University, Beijing, China

Abstract Esophageal cancer (EC) is a formidable health issue due to its high aggressiveness and poor prognosis. Despite advancements in treatment modalities, the survival rate remains dismal, often attributed to late-stage diagnosis and resistance to therapy. Mast cells (MCs), tissue-resident immune cells, play a significant role in the tumor microenvironment (TME) and have been implicated in both tumor-promoting and tumor-suppressing activities. Recent developments in high-throughput single-cell RNA sequencing (scRNA-seq) have allowed for the detailed dissection of tumor heterogeneity and the complex interactions between cancer cells and their microenvironmental components. However, the application of scRNA-seq in understanding the immune landscape of EC remains limited. In this study, we employed high-throughput scRNA-seq to investigate the TME of EC, focusing on the role of MCs. We analyzed six samples from patients with EC, identifying 11 major cell types and highlighting the significant presence of MCs in pericarcinoma tissues. Further sub-clustering of 5,001 MCs revealed eight distinct subtypes, among which MCs with high SRSF7 expression showing a strong tumor preference and potential tumor-promoting characteristics. Our findings underscore the complex interactions within the TME and suggest that targeting MC subtypes, particularly through the EGFR signaling pathway, could offer new therapeutic avenues. This study provides a comprehensive immune map of EC and lays the groundwork for developing more effective treatment strategies. 1. Introduction Esophageal cancer (EC) is a common malignant tumor of the gastrointestinal system, with the seventh highest incidence and sixth highest mortality rate in the world(1). In China, the incidence and mortality rates of EC rank third and fourth, respectively, among all malignant tumors(2). Despite the development of a multidisciplinary treatment approach, the prognosis remains unfavorable(3). The 5-year survival rate for EC is only 21%, after pancreatic and liver cancers(4). Therefore, EC has been a major malignant tumor threatening the health of Chinese residents. EC consists of two main subtypes, esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma, with ESCC accounting for about 90% of all EC cases worldwide(5). EC is an aggressive cancer with rapid growth and a high rate of lymph node metastasis, usually involving the upper two thirds of the esophagus(6). In retrospective studies in EC, smoking, hot tea consumption, red meat consumption, poor oral health, low intake of fresh fruits and vegetables, and low socioeconomic status were associated with a higher risk of EC(7). Previous studies have shown that chronic inflammation plays a central role in progression from esophageal precancerous lesions (EPL) to esophageal squamous cell carcinoma (ESCC), that dietary inflammatory potential has been linked to both EPL and ESCC, and that inflammatory imbalances promote tumorigenesis, and that the consumption of anti-inflammatory foods may be helpful in the prevention of EPL and ESCC(8-10). Difficulty swallowing and swollen lymph nodes in the neck do not appear until the cancer has progressed to an advanced stage(11), and the treatment of EC patients faces major challenges due to the lack of early symptoms, high malignancy, poor prognosis, and surgical complexity of EC. Although we have made great progress in the treatment of EC in recent years, especially through preoperative radiotherapy combined with immunotherapy, which shows a broad potential in the treatment of EC. However, due to the high rate of post-treatment recurrence and the limitations of drugs and treatment strategies after metastasis, only a small proportion of EC patients can benefit from the available treatments, while the majority of patients respond poorly to the treatments, and therefore, the overall survival rate of EC is still disappointing in China(3,12). In addition, due to the heterogeneity and complexity of tumors, the mechanisms of tumor proliferation, metastasis, drug resistance, and immunosuppression are unknown. Therefore, elucidating the molecular mechanisms of tumorigenesis and tumor progression is crucial for effective control and management of tumor development. Notably, the presence of non-tumor cells within the tumor tissue is also critical for tumor development(13). Therefore, shifting the therapeutic focus to other components of the tumor microenvironment (TME) may become an important strategy for future tumor therapy. The introduction of TME has played a very powerful role in advancing oncology research. TME has had an incredibly important role in the development and evolution of EC(14). The TME consists of multiple cellular components (e.g., fibroblasts, endothelial cells, and immune cells) and extracellular components (including cytokines, hormones, extracellular matrices, and growth factors), which form a complex network that encapsulates EC cells. These cells shape cancer biology and influence the response to treatment(15-17). In TME, mast cells (MCs) are tissue-resident immune cells that are important players in diseases associated with chronic inflammation such as cancer. Because MCs can infiltrate solid tumors and promote or limit tumor growth, MCs may polarize to either pro- or anti-tumor phenotypes and remain a challenging area of research(18). Previous articles have also hypothesized that NRF2 in combination with AC-MCs may be a predictive marker for prognosis and may influence immunotherapy by modulating PD-L1 in EC(19). High-throughput single-cell RNA sequencing (scRNA-seq), developed in recent years, is an effective method that has been shown to dissect heterogeneous tumors and decipher the interactions between cancer cells and their microenvironmental components, and to elucidate the transcriptomic profiles of both the cancer cells and the microenvironmental components(20-22), which is the basis and foundation for furthering the understanding of cancers and the development of effective early diagnostic and therapeutic strategies, previous studies have dissected the esophageal squamous cell carcinoma ecosystem by single-cell transcriptomic analysis(23), but its application in EC is still limited. At the same time, there is still a long way to go for early detection of esophageal cancer(24), and prognostic tools lack the necessary accuracy to facilitate individualized patient management strategies(25). Therefore, in this study, scRNA-seq was used to sequence EC samples in order to decipher the immune microenvironment of EC, reveal the immune map of EC, and provide new insights for the treatment of EC. We also performed knockdown experiments to verify the importance of the EGFR signaling pathway in the crosstalk between MCs and tumor cells, providing new ideas for the targeted treatment of MCs. It has also been shown that the receptor we studied, EGFR, is a strong prognostic indicator in head and neck, ovarian, cervical, bladder and esophageal cancers(26). The functional role of mast cell subtypes in EC and their association with tumor tissues are extensively discussed and summarized in this paper, and a prognostic model is established, which provides a valuable resource for deeper understanding of the causes and progression of EC and helps to improve its therapeutic strategies. 2. Methods and Materials 2.1 Data source The scRNA-seq data of EC were acquired from the GEO website (https://www.ncbi.nlm.nih.gov/geo/) under the accession number GSE196756. Considering the utilization of publicly accessible data derived from databases, it was not required to secure an ethical endorsement for this investigation. 2.2 Single‐Cell Sequencing The gene expression data were imported into the R software and analyzed using the Seurat R package(27). Cells of inferior quality were excluded based on the following criteria: (1) nFeature between 300 and 7,500; (2) nCount between 500 and 100,000; (3) mitochondrial gene expression occupying no more than 25% of the total gene count within the cell; (4) erythrocyte gene expression not surpassing 5% of the total gene count within the cell. The harmony R package(28,29) was employed to alleviate the influence of batch effects among the samples. The dim value was set to 30, while the resolution parameter was configured to 1.2. Subsequently, all gene expression data underwent normalization and scaling using the "Normalize Data" and "Scal Data" functions within the Seurat R package(30). For the purpose of principal component analysis, the "Find Variable Feautres" function(31) was implemented to identify the top 2,000 most variable genes. These cells were then segregated into clusters based on the top 30 principal components (PCs) using the "Find Clusters" function at a resolution of 1.0. Finally, the top 30 significant PCs were selected to dimensionality reduction and visualization of gene expression through uniform manifold approximation and projection (UMAP)(32,33). 2.3 Identification of cell subtypes Cell clusters were initially discerned utilizing the "Find Clusters" and "Find Neighbors" functions within Seurat, employing a default resolution of 0.8. Afterwards, these cell clusters were bestowed with annotations based on the average gene expression of representative markers. In order to evaluate differentially expressed genes (DEGs) across distinct cell clusters, a Wilcoxon rank sum test was employed utilizing Seurat's "Find All Markers" function. The parameters min.pct and min.diff.pct were established at 0.25, while the LogFc threshold was set to 0.25. 2.4 Cancer preferences analysis In order to evaluate the predilection of MCs subtypes for cancer, odds ratios were computed utilizing the calculation methodology(34). 2.5 Trajectory analysis of MCs subtypes The slingshot R package was employed to deduce cellular lineages and pseudotimes. It delineated the structure of lineages through clustering-based minimum spanning trees and employed synchronized master curves to model branching trajectories for these lineages. The "get Curves" function was utilized to acquire refined trajectory curves. The association between gene expression and pseudotime was characterized by modeling the noise distribution of each gene through a generalized additive model with negative binomials. This approach allowed for the simulation of genes exhibiting a gradual alteration in expression throughout the pseudotime continuum(35). 2.6 Assessment of cell stemness AUCell(36) represents a novel approach to discerning cells harboring active genes within single-cell RNA-seq datasets. Given a gene set as input, it provides an evaluation of the "activity" exhibited by that particular gene set in each individual cell. In the context of this study, AUCell was employed to quantitatively assess the level of stemness exhibited by various subtypes of MCs. To hypothesize the temporal trajectory of cell differentiation, the CytoTRACE R package was utilized(37). 2.7 Enrichment analysis of cellular subtypes By leveraging the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Genome Enrichment Analysis (GSEA) tools, available at http://software.broadinstitute.org/gsea/msigdb, within the Cluster Profiler R package(38), we carried out enrichment analysis on the DEGs. To discern the disparities among various risk groups within the bulk data, the DESeq2 R package was applied, employing a threshold of |logFC| > 2 and a p-value threshold below 0.05. 2.8 Cell communication analysis The CellChat R package(39) was used to analyze complex cell-to-cell interactions and develop regulatory networks based on ligand-receptor expression. The "netVisual DiffInteraction" function was applied to depict differences in communication strength among cells, and the "Identify Communication Patterns" function was utilized to estimate the number of communication patterns. A significance threshold of 0.05 was set. Various visualizations, including circle plots, bubble plots, and violin plots, were used to represent the incoming and outgoing signals of all cells 2.9 Scenic analysis In evaluating the transcriptional activity within diverse subtypes of tumor cells, we employed the SCENIC analysis with Python. 2.10 Cell culture Cell lines TE-10 and KYSE-30 were acquired from the American Type Culture Collection. The TE-10 cell line was grown under standard conditions (37°C, 5% CO2, 95% humidity) in RPMI1640 media with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin. KYSE-30 cell line was grown under standard conditions (37°C, 5% CO2, 95% humidity) in RPMI1640 media with 10% FBS, 1% penicillin-streptomycin, and 1% sodium pyruvate. 2.11 Cell transfection EGFR knockdown was accomplished through the use of GenePharma (Suzhou, China) small interfering RNA (siRNA) constructs. According to Lipofectamine 3000 RNAiMAX (Invitrogen, USA) manufacturer's instructions, transfection was carried out. Two knockdown constructs (Si-EGFR-1 and Si-EGFR-2) and a negative control (si-NC) were transfected into cells that had been plated at 50% confluency in six-well plates. Every transfection was carried out using Lipofectamine 3000 RNAiMAX (Invitrogen, USA). 2.12 Cell viability assay Using the CCK-8 assay, the cell viability of transfected AGS and SGC-7901 cells was evaluated. After being cultivated for 24 hours, cells were planted at a density of 5×10³ cells per well in 96-well plates. Following the addition of 10μL of CCK-8 reagent (A311-01, Vazyme) to each well, the plates were incubated for two hours at 37°C in the dark. On days 1, 2, 3, and 4 post-transfections, absorbance at 450 nm was measured using a microplate reader (A33978, Thermo). Plotting of the mean OD values was done. 2.13 5-Ethynyl-2'-deoxyuridine (EdU) proliferation assay In 6-well plates, 5×10³ cells were planted per well with transfected CNE2 and HNE2 cells, and they were grown for an entire night. A 2x EdU working solution was then created by combining serum-free medium with 10 mM EdU. Following two hours of incubation at 37°C, the cells were rinsed with PBS, fixed for thirty minutes with 4% paraformaldehyde, permeabilized for fifteen minutes with a solution of 2 mg/mL glycine and 0.5% Triton X-100, then stained for thirty minutes at room temperature using a solution of 1X Apollo and 1X Hoechst 33342. The measurement of cell proliferation was done by fluorescence microscopy. 2.14 Wound-healing assay In 6-well plates, stabilized transfected cells were plated and allowed to grow to confluence. Each well was scratched with a sterile 200μL pipette tip, and then it was cleaned with PBS to get rid of any remaining cell debris before being incubated in a medium without serum. Using Image-J software, the breadth of the scratches was measured after they were photographed at 0 and 48 hours. 2.15 Transwell assay Before the experiment, cells were fasted for 24 hours in a serum-free medium. The upper chamber of Costar plates was filled with cell suspension after being treated with Matrigel (BD Biosciences, USA), while the lower chamber was filled with media containing serum. In a cell culture incubator, the cells were incubated for 48 hours. To evaluate the cells' ability to invade, they were fixed with 4% paraformaldehyde after incubation and stained with crystal violet. 2.16 Construction and Validation of the Prognostic Model We determined the most important predictive genes using LASSO regression analysis and univariate Cox analysis. The risk coefficients for each prognostic gene were then determined using multivariate Cox regression analysis, allowing for the creation of a risk score model: Risk score = ∑_i^n▒〖Xi×Yi〗. X stands for the coefficient and Y represents the gene expression level. Using the "surv-cutpoint" function to compute the best cutoff value, patients were divided into two groups: low-risk and high-risk. We also used the Survival R package for survival analysis of the created risk score model and the "ggsurvplot" function (30) to depict survival curves in order to observe the prognostic outcomes in various patient cohorts. ROC curves were plotted using the timeROC R package to assess the predictive model's accuracy and calibration. 2.17 Kaplan-Meier survival curve of selected genes We performed a survival analysis utilizing the R packages survminer and survival. The area under the ROC curve (AUC) was calculated after generating ROC curves for 1-year, 3-year, and 5-year survival rates using the Survive and Time ROC R packages. Model validation was conducted through survival analysis and time-dependent ROC analysis. To evaluate the model, we employed a heatmap, a scatter plot of survival status, and a distribution of risk scores. 3. Result 3.1 ScRNA sequencing revealed the main cell types in the EC To identify the major cell types during the progression of EC, we collected pericarcinoma and tumor tissue samples from three EC patients for single-cell RNA sequencing (scRNA-seq). We also checked the quality and completeness of the raw data. This included checking for missing values, outliers, or any anomalies that might affect the analysis. We excluded genes in the sample that did not meet the minimum expression threshold. For example, genes with low counts or low variability were excluded as they may not provide meaningful insights. After performing initial quality control and removing batch effects, we retained a total of 29,719 cells. We categorized these 29,719 cells into 30 cell clusters by dimensionality reduction (Figure 1A). According to the cell gene map and typical markers, 30 cell clusters were finally identified into 11 cell types, including B-Plasma cells (IGKC), T-NK cells (IL32), mast cells (MCs, TPSB2), neutrophils (S100A8), fibroblasts (DCN), myeloid cells (LYZ), epithelial cells (EPCs, KRT5), proliferating-cells (MKI67), endothelial cell (ECs, AQP1), smooth muscle cell (SMCs, MYH11), and neurons (NRXN1). From the pie charts and bar graphs, we could learn that for tissues, T-NK cells accounted for the largest proportion in tumor tissues, followed by B-Plasma cells, while MCs were the most predominant cell type in pericarcinoma tissues; for phases, T-NK accounted for the largest proportion in both the G2M and the S phases, while on the contrary, most of MCs accounted for the largest proportion in the G1 phase (Figure 1B-C). Figure 1D showed the top 5 marker genes for 11 cell types. UMAP and violin plots were utilized to visualize nCount-RNA, nFeature-RNA, G2M.Score, and S.Score across all cells, demonstrating that proliferating cells exhibited the highest proliferative activity and vigorous division (Figure 1E-G). At the same time, the distribution of marker genes on UMAP for each cell type was presented (Figure 1H). Among all cell types, MCs drew our attention. MCs play a crucial role in allergic reactions, pathogen immune responses during infections, angiogenesis, and the regulation of both innate and adaptive immunity. In addition to all these roles, MCs were increasingly recognized as regulators of the tumor microenvironment. Despite the accumulating evidence for MCs in tumors, their exact role in the tumor microenvironment remained incompletely understood(40). Therefore, we next performed a further analysis of mast cells. 3.2 Visualization of MCs subtypes in EC Next, we analyzed the scRNA-seq data from tumor and pericarcinoma tissues, identified MCs, and performed further sub-clustering. This analysis resulted in eight distinct cell subtypes from a total of 5,001 mast cells: C0 EGR1+ MCs, C1 SRSF7+ MCs, C2 TXNIP+ MCs, C3 DUSP1+ MCs, C4 S100A8+ MCs, C5 HSPA6+ MCs, C6 IL32+ MCs, C7 RPL35A+ MCs (Figure 2A), and showed the distribution of phases and tissue sources in the subtypes, while faceting gived a clearer picture of the distribution of MCs from different sample sources (Figure 2B-D). The bar graphs illustrated that the C1 SRSF7+ MCs had the highest proportion of tumor tissues of P1 and P3 origin and was enhanced over the pericarcinoma tissues share, and similarly, the Ro/e preference graph corroborated this, suggesting that the C1 SRSF7+ MCs was more preferred to tumor tissues (Figure 2E-F). In order to better explore the characteristics of different MCs subtypes, we visualized their typical genes. As shown in Figure 2G, C1 SRSF7+ MC highly express DDX5, which had been shown to be associated with a variety of key tumor promoting molecular interactions and was involved in tumorigenesis and tumor progression signaling pathways(41). This suggested that C1 SRSF7+ MCs in EC might be involved in tumor promoting effect. Several related features (CNVscore, ncount-RNA, S.Score and G2M.Score) of eight MCs subtypes were visualized (Figure 2H-I). From the Figures, we could learn that C7 RPL35A+ MCs had the highest expression level of CNVscore and G2M.Score, while C1 SRSF7+ MCs and C6 IL32+ MCs had the highest nCount-RNA expression level, and all subtypes had basically the same expression level of S.Score. In the end, bar plots showed the expression level of gene makers in each MCs subtype, validating the basis for delineating subtypes (Figure 2J). 3.3 Slingshot analysis of proposed temporal trajectories of MCs subtypes To infer the lineage trajectory and pseudotime sequence of MCs, we employed slingshot analysis to assess the distribution of MCs differentiation trajectories across all MCs, visually represented through UMAP plots (Figure 3A). Then, we found 3 cell lineage trajectories of the MCs subtypes (Figure 3B-E). Including: lineage 1: C4 → C2 → C3 → C0; lineage 2: C4 → C2 → C3 → C1; lineage 3: C4 → C2 → C3 → C6. Slingshot analysis revealed that the differences among the three trajectories mainly reside in the late stages. Combined with Figure 3C-E to determine, lineage 1's endpoint was located in C0, which showed no preference for tumor tissue, lineage 3's endpoint was located in C6, which had a very small number of cells and a low percentage of tumor tissue, while lineage 2's endpoint was located in C1, which not only showed a preference for tumor tissue, but also had a high percentage of tumor tissue. Therefore, we concluded that lineage 2 represented the differentiation line of MCs associated with the tumor. In addition, we also noted that MCs are influenced by some cytokines or tumor cell-secreted proteins during development in TME, resulting in a possible transformation of the MCs phenotype to a tumor-associated or pro-tumorigenic phenotype(18), whereas C1 belonged to the terminal end and consisted predominantly of MCs originating from tumor tissues, and based on this observation, we hypothesized that C1 may play a crucial role in the differentiation of tumor-associated MCs(TAMCs) process. Subsequently, we confirmed the biological processes corresponding to the three cell lineage trajectories of MCs subtypes using GO-BP enrichment analysis (Figure 3F). It was found that C1 in lineage 2 was associated with biological processes such as endopeptidase and cysteinetype, C2 was linked to processes such as protein folding, C3 was related to leukocyte functions, and C4 was involved in processes such as lamellipodium formation, contraction, and production. Finally, the dynamic trends plot demonstrated the expression variation and distribution of marker genes for MCs subtypes along the three differentiation trajectories in pseudotime (Figure 3G). 3.4 Expression of stemness gene sets in subtypes of MCs To investigate the expression of stem cell genes in MCs subtypes and to understand their differentiation potential, we used bubble plots to illustrate the different expression of stem cell genes in MCs subtypes. The results showed the corresponding expression of stem cell genes KDM5B, EPAS1, CTNNB1, EZH2, KLF4, CD44, BMI1, and HIF1A in MCs subtypes and different tissue types (Figure 4A). Subsequently, we visualized the cell stemness AUC scores of different MCs subtypes using a UMAP plot (Figure 4B). We then combined this with other analyses to assess the expression levels of stemness-related genes in different subtypes of MCs, and violin plots showed the different expression levels of stemness genes in different sample sources, tissues, subtypes of MCs, and phases, respectively (Figure 4C-F). The results showed that C1 SRSF7+ MCs exhibited a higher level of cell stemness, indicating a lower degree of differentiation and higher differentiated potential; and it also showed that pericarcinoma tissues had the higher level of cell stemness. In addition, there was no significant difference in the expression levels of stemness genes in different cell phases. By CytoTRACE analysis, C1 SRSF7+ MCs showed the lowest degree of differentiation and the highest cell stemness among all subtypes, which we hypothesized might be related to the transformation of MCs to TAMCs (Figure 4G-H). Afterwards, the expression profiles of stemness genes with relatively elevated expression levels in Figure 4A were demonstrated in all MCs by UMAP plots (Figure 4I). 3.5 Enrichment analysis of MCs subtypes in EC First, we utilized volcano plots to represent the DEGs profiles between subtypes of MCs (Figure 5A). The results showed that the up-regulated DEGs in C1 SRSF7+ MCs were mainly DDX5, EEF1A1, TPSB2, TPSAB1, and CPA3. In addition, we performed GO-BP enrichment analysis of the DEGs in the subtypes of MCs to reveal their enrichment in biological processes. The heatmap showed the results of the top four enrichment entries in the MCs subtypes (Figure 5B). The C0 EGR1+ MCs subtype was mainly associated with pathways such as response to unfolded protein, response to topologically incorrect protein and regulation of neuron death; The C1 SRSF7+ MCs subtype was enriched in pathways such as protein folding, protein refolding, chaperone-mediated protein folding and ‘de novo’ protein folding; The C2 TXNIP+ MCs subtype revealed their close association with cytoplasmic translation, oxidative phosphorylation, ribosome biogenesis and rRNA processing; The C3 DUSP1+ MCs subtype showed enrichment in pathways such as negative regulation of transferase activity, response to muscle stretch, response to mechanical stimulus and negative regulation of phosphorylation; The C4 S100A8+ MCs subtype was enriched in pathways related to leukocyte migration, myeloid leukocyte migration, response to molecule of bacterial origin and leukocyte chemotaxis; The C5 HSPA6+ MCs subtype mainly exhibited enrichment in pathways such as protein refolding, response to temperature stimulus, myeloid cell differentiation and regulation of hemopoiesis; The C6 IL32+ MCs subtype revealed pathways related to leukocyte mediated cytotoxicity, lymphocyte mediated immunity, natural killer cell mediated immunity and positive regulation of leukocyte cell-cell adhesion; The enrichment analysis conducted on the C7 RPL35A+ MCs subtype revealed their association with cytoplasmic translation, ribosomal small subunit biogenesis, rRNA processing and rRNA metabolic process. The word cloud plots illustrated the enrichment results of DEGs across various pathways in the eight MC subtypes (Figure 5C). The results showed that the C1 SRSF7+ MCs subtype was mainly enriched in leukocyte, immune and activation, and it was hypothesized that C1 SRSF7+ MCs subtype might be related to MCs activation and participation in immune regulation. In addition, the results of GSEA enrichment analysis were also shown in the form of bubble plots (Figure 5D). It showed that C1 SRSF7+ MCs subtype was significantly expressed in regulation of immune system process, cell motility and migration, protein folding, response to immune and external stimulus pathways. All of the above pathways would suggest that MCs in the C1 SRSF7+ MCs subtype had likely transformed into TAMCs. Finally, we performed GSEA on the DEGs of the C1 SRSF7+ MCs subtype according to GO-BP terminology. The results were shown in Figure 5E. We observed that pathways associated with protein refolding, skeletal muscle cell differentiation, chaperone cofactor-dependent protein refolding and ‘de novo’ protein folding were upregulated in the C1 subtype. In contrast, pathways associated with ATP synthesis coupled electron transport, mitochondrial ATP synthesis coupled electron transport, aerobic electron transport chain and cytoplasmic translation were downregulated in the C1 subtype. Combining the above up-regulated genes and enriched pathways with previous studies, we believed that the C1 subtype was affected by the endoplasmic reticulum stress state, which disrupts the homeostasis of the original proteins and generates aberrant protein folding, and that this stress state could control a variety of pro-tumorigenic attributes of cancer cells, dynamically re-programming the function of immune cells, transforming MCs into TAMCs, thus exerting pro-tumorigenic effects, and conferring a greater tumorigenic, metastatic, and drug-resistant capacity to the malignant cells(42). 3.6 Transcription factors regulate the carcinogenic mechanism of C1 SRSF7+ MCs Transcription factors can directly act on the genome and regulate gene transcription and affect the biological function of cells by combining specific nucleotide sequences in the upstream of the gene. Therefore, we used scenic to analyze the gene regulatory network of C1 SRSF7+ MCs. First of all, we carried out cluster analysis of MCs according to regulator activity (Figure 6A). It was obvious that the discretization of UMAP diagram based on regulator activity was smaller, the interference factors were better excluded, and all MCs were clustered and distributed. Among them, C1 SRSF7+ MCs were mainly distributed on the right side of UMAP plot without significant discretization. By further analyzing the key regulators of different MCs subtypes, the five major regulators of C1 SRSF7+ MCs, ATF4, JUNB, NFκB2, MAFK and JUN, were identified (Figure 6B-C). After analyzing these five key regulators in depth in conjunction with previous studies and Figure 6D, ATF4 and JUNB caught our attention. ATF4, which was expressed at higher levels in C1 SRSF7+ MCs than in other subtypes, was a major transcriptional regulator of the unfolded protein response to hypoxia, activated genes that promoted recovery of normal endoplasmic reticulum function and hypoxic survival(43), regulated mast cells through endoplasmic reticulum stress(44), and had been associated with programmed cell death in a variety of tumors, particularly ER stress-induced iron death(45,46). As for JUNB, its expression level was high in C1 SRSF7+ MCs, C4 S100A8+ MCs and C5 HSPA6+ MCs subtypes, and it is a potent inhibitor of endoplasmic reticulum stress and apoptosis, and, in particular, its modulation of endoplasmic reticulum stress is associated with ATF4 alterations(47). 3.7 CellChat analysis among all cells In order to systematically elucidate complex cellular responses, we aimed to investigate cell-to-cell relationships and ligand-receptor communication networks to better understand interactions between cells. Using CellChat analysis, we initially established intercellular communication networks involving various cells such as tumor cells, fibroblasts, T-NK cells, and different subtypes of MCs, etc (Figure 7A). After establishing the intercellular communication networks using CellChat analysis, we calculated both the number of interactions (represented by the thickness of the connecting lines between two cell types) and the strength of interactions (indicated by the weight of the lines, where thicker lines denote stronger interaction strengths). This approach helped quantify the complexity and intensity of communication pathways between different cell types in the network. We utilized gene expression pattern analysis methods available through CellChat to investigate how cells and signaling pathways interact. Initially, we assessed the relationship between inferred potential communication patterns and groups of cells that secrete signaling molecules to decipher outgoing communication patterns. Three distinct signaling patterns were identified through our analysis: pattern 1 (subtypes of MCs), pattern 2 (Neurons cells, fibroblasts, SMCs, tumor-cells and ECs) and pattern 3 (myeloid-cells, B-Plasma cells, neutrophils, proliferating-cells and T-NK cells) (Figure 7B). To identify the key incoming and outgoing signals associated with the eight MCs subtypes, we quantitatively analyzed the ligand-receptor network using CellChat. This approach allowed us to predict the primary incoming signals from secreting cells (signal senders) releasing various cytokines or ligands. Additionally, we assessed which cell types acted as targeting cells (signal receivers), and how ligand-receptor-mediated communications between different cell types contributed to the progression of EC. This analysis helped illustrate how receptors on these cells were targeted by ligands released either from the same type of cell or from other cell types (Figure 7C). In addition to examining detailed communication within individual pathways, an important aspect was understanding how multiple cell populations and signaling pathways coordinate their functions. To address this, CellChat employed a pattern recognition method based on nonnegative matrix decomposition. This method identified global communication patterns and key signals across different cell groups, shedding light on how various cells and pathways collaborate in their functions. The application of this analysis revealed three distinct incoming signaling patterns and three outgoing signaling patterns. For instance, this output indicated that the majority of outgoing MCs signaling was characterized by pattern 1, which represented multiple pathways, including but not limited to CD99, ANNFXIN, EGF, PARs, ICAM, CSF, etc. All output tumor-cells, fibroblasts, ECs, SMCs, neurons signalings were characterized by pattern 2, which represented pathways such as COLLAGEN, LAMININ, FN1, APP, PTN and so on. On the other hand, the analysis of communication patterns in target cells indicated that incoming signalings to tumor-cells, SMCs, and neurons were predominantly characterized by pattern 1. This pattern included signaling pathways such as EGF, TENASCIN, JAM, MPZ, CADM, and TWEAK. In contrast, the majority of incoming signalings to subtypes of MCs, B-plasma cells, T-NK cells, proliferating-cells, myeloid-cells, and neutrophils were characterized by pattern 2, which was driven by pathways such as CXCL and ANNEXIN (Figure 7D). Combining the above analysis and the demonstration of all incoming and outgoing signal intensities in Figure 7E, the signaling molecule EGF caught our attention. EGF was present in the incoming pathway of tumor-cells, i.e., tumor-cells were the target cells, and EGF is again present in the outgoing pathway of C1 SRSF7+ MCs subtype, i.e. C1 SRSF7+ MCs subtype was the secreting cell, which links C1 SRSF7+ MCs subtype and tumor-cells, we speculated that this signaling pathway might be related to tumor progression, so we next focused on EGF. 3.8 Analysis of AREG-EGFR/AREG-(EGFR+ERBB2) signal pathway The circular displayed the inferred cell-cell communication network between MCs and other cells (Figure 8A-B). The results showed that there was a strong crosstalk between C1 SRSF7+ MCs and tumor cells. We considered all identified cell types in ECEC as source cells for the EGF signaling pathway, and the results indicated that all subtypes of MCs could target tumor cells with released EGF. In addition to the senders and receivers of EGF signaling, based on the relative importance of each cell type in EGF signaling-mediated intercellular communication, we identified the cell types that act as mediators and influencers in this process, which is referred to as the "centrality measurement" algorithm. As can be seen from the Figure, C1 SRSF7+ MCs subtype had higher expression as a ‘sender’ in the EGF signaling pathway, whereas tumor-cells were acting as ‘receiver’, ‘mediator’ and ‘influencer’ in this signaling pathway (Figure 8C). Similarly, the heatmap corroborated this conclusion (Figure 8D). The violin plot showed the cell-cell interactions while giving the different ligands and receptors in the EGF signaling pathway, and the results showed that C1 SRSF7+ MCs subtype and tumor-cells were mainly contacted with AREG as a ligand and EGFR or ERBB2 as receptors (Figure 8E). Bubble and circle plot as well as hierarchical plots likewise corroborated this conclusion (Figure 8F-H). Combined with the results of previous results in this paper, it can be concluded that the C1 SRSF7+ MCs and tumor cells crosstalk through the AREG-EGFR/AREG-(EGFR+ERBB2) signal pathway, thereby exerting a tumor-promoting effect. 3.9 In vitro experimental validation of EGFR To further investigate the role of EGFR in EC, we conducted in vitro experiments using the TE-10 and KYSE-30 cell lines. Initially, we knocked down EGFR and measured the mRNA and protein expression levels before and after knockdown. We observed a significant reduction in both mRNA and protein expression levels in both cell lines compared to the control group (Figure 9A). Subsequently, the CCK-8 assay revealed a marked decrease in EC cell viability post-EGFR knockdown (Figure 9B). Colony formation assays and EDU experiments confirmed that EGFR knockdown inhibited EC cell proliferation (Figure 9C, E-F). Additionally, scratch and transwell assays were employed to assess the migration and invasion capabilities of EC cells post-EGFR knockdown, demonstrating a significant reduction in migration and invasion levels (Figure 9D, F-H). These results collectively indicate that EGFR knockdown suppresses the activity, proliferation, migration, and invasion of EC cells, thereby inhibiting tumor growth. 3.10 Enrichment Analysis and Construction of Predictive Models To further investigate the impact of MCs with high SRSF7 expression on EC patients, we divided the TCGA cohort patients into high and low SMRS (SMRS: SRSF7+MCs risk score) groups according to the gene expression levels of the SRSF7+ MCs subtype. A heatmap illustrated the expression profiles of the top 30 DEGs (Figure 10A), and a volcano plot depicted the up-regulation and down-regulation of DEGs (Figure 10B). Subsequently, we employed various enrichment methods to gain insights into the associated biological processes. KEGG enrichment analysis revealed that DEGs were primarily enriched in pathways such as cholesterol metabolism, PPAR signaling pathway, and fat digestion and absorption (Figure 10C). Whereas a large number of studies have shown that cholesterol metabolism has functional significance in tumorigenesis, cancer progression and metastasis through its regulatory role in immune response, iron death, autophagy, cell stemness and DNA damage response(48). In GO-BP analysis, enrichment was observed in the triglyceride metabolic process, acylglycerol metabolic process, and neutral lipid metabolic process (Figure 10D). In GO-CC analysis, enrichment included chylomicron and high- and low-density lipoproteins (Figure 10E). And in GO-MF analysis, glycosaminoglycan binding and lipoprotein particle receptor binding were highlighted, insterestingly, dysregulated expression of the former was present in cancer and has been reported to correlate with clinical prognosis in several malignancies(49) (Figure 10F). We then visualized the primary enrichment terms for each gene set and used t-SNE plots to graphically represent the risk score distribution of these enrichment terms (Figure 10G). GSEA results showed that the up-regulated genes were mainly enriched in processes such as intestinal absorption, peptidyl methionine modification, intestinal lipid absorption and protein lipid complex assembly, while down-regulated genes were enriched in processes like regulation of release of sequestered calcium on into cytosol, external encapsulating structure organization, B cell receptor signaling pathway and collagen fibril organization (Figure 10H). Additionally, we constructed a prognostic model to explore the clinical significance of MCs with high SRSF7 expression. Univariate Cox regression analysis identified 11 genes significantly associated with prognosis (Figure 10I), with AHR as a protective factor (HR < 1) and the others as risk factors. To address the issue of multicollinearity among these genes, we further screened them using LASSO regression analysis, ultimately identifying eight prognostic-related genes (Figure 10J). Cox regression analysis was then used to calculate the coefficient values of these genes (Figure 10K-L). Curve and scatter plots demonstrated the differences in risk scores and survival outcomes between the two groups, indicating that the high SMRS group was associated with poorer prognosis (Figure 10M). Moreover, a heatmap displayed the differential expression patterns of genes used in model construction (Figure 10N). Kaplan-Meier survival curves further confirmed the conclusion that the high SMRS group had a worse survival outcome (Figure 10O). ROC curves and AUC values for 1-year, 3-year, and 5-year outcomes indicated that the model had good predictive value (Figure 10P). 4. Discussion In recent years, the rapid development and application of scRNA-seq in cancer research has revolutionized our understanding of the biological features and dynamics within cancer lesions, greatly facilitating the diagnosis, treatment, and prognosis prediction of a range of tumors(50-52). Overall, the present study focused on mast cells in esophageal cancer, and we validated the pro-carcinogenic role of this pathway by launching a comprehensive profiling of mast cell subtypes with an eye on the C1 SRSF7+ MCs and obtaining its reciprocal receptor, EGFR, using cellular communication analysis, and subsequently verifying the pro-carcinogenic role of this pathway through cellular knockdown experiments. In this study, we comprehensively characterized the cellular heterogeneity of EC using scRNA-seq technology. We identified immune cells including T-NK cells, MCs, and myeloid cells and so on, as well as non-immune cells such as smooth muscle cells and neuronal cells. In addition, we carefully analyzed the sample origin of these cell types and the distribution characteristics during the phase. Among them, MCs caught our attention. Until more than a hundred years ago, MCs were regarded as effectors of allergy, and it is only in the last two decades that MCs have gained recognition for their involvement in other physiological and pathological processes. MCs maturation, phenotype and function as a direct result of the local microenvironment(53), and by releasing a range of bioactive mediators has a significant effect on their ability to specifically recognize and respond to a variety of strategies(54-56). Therefore, depicting and analyzing the TME is important for MCs. And in previous studies, MCs have been shown to correlate with pro-tumorigenic effects(57-59). Despite the accumulating evidence for MCs in tumors, their exact role in the TME remains incompletely understood(40). We therefore focused our attention on the study of MCs. By further dimensionality reduction clustering, we obtained eight MCs subtypes, i.e., C0 EGR1+ MCs, C1 SRSF7+ MCs, C3 TXNIP+ MCs, C4 S100A8+ MCs, C5 HSPA6+ MCs, C6 IL32+ MCs, and C7 RPL35A+ MCs. By integrating the proportions of MCs subtypes in sample sources and cell phases, Ro/e preference analyses, cell stemness analyses, and slingshot proposed pseudotime analyses, we identified the target subtype in this study: the C1 SRSF7+ MCs. C1 SRSF7+ MCs were significantly more abundant in tumor tissues than in pericancer tissues in P1 and P3 samples, and this was confirmed by Ro/e preference analysis. In slingshot proposed pseudotime analysis, Lineage 2 was considered to be representative of the differentiation trajectory of MCs associated with tumors. And the endpoint of Lineage 2 was a subtype of C1 SRSF7+ MCs, this result may prove that MCs are affected by some cytokines or tumor cell-secreted proteins during their development in the TME, leading to the transformation of MCs into a tumor-associated or pro-tumor phenotype, which is in line with the previous study(60). Meanwhile, cell stemness analysis by AUC value and CytoTRACE showed that the C1 SRSF7+ MCs subtype had the strongest cell stemness among all subtypes, with high differentiation potential, which did not contradict slingshot's results, and it is understandable that the transformation from normal phenotype to TAMCs phenotype would result in an increase in cell stemness. It can be seen that the C1 SRSF7+ MCs subtype is intricately linked to tumor progression. To further investigate the tumor-promoting related roles of the C1 SRSF7+ MCs subtype, we performed enrichment analysis and obtained the upregulated genes DDX5, TPSB2, and CPA3, of which DDX5 interacts with a variety of key pro-tumorigenic molecules and participates in tumorigenic and tumor progression signaling pathways, and when DDX5 is expressed or its post-translational modifications are deregulated, the normal cellular signaling network collapses, leading to many pathological states, including tumorigenesis and tumor progression(61,62). Moreover, the enriched pathways obtained by GO-BP and GSEA on the C1 SRSF7+ MCs subtype showed that the C1 SRSF7+ MCs subtype was extensively involved in protein folding and refolding, regulation of immune system processes, and response to external stimuli. All these pathways suggest that the C1 SRSF7+ MCs subtype has probably been transformed into TAMCs. Finally, combining the above up-regulated genes and enriched pathways, we suggest that the C1 SRSF7+ MCs subtype is affected by the endoplasmic reticulum stress state(63), which disrupts the original protein equilibrium(64) and produces aberrant protein folding(65,66), and this stress state dynamically reprograms the function of MCs, transforming MCs into TAMCs, which exerts pro-tumorigenic effects(67) and confers cancer cells with enhanced tumorigenic, metastatic, and drug-resistant capabilities. In this regard, we can treat patients with esophageal cancer by targeting the abnormal protein folding to prevent MCs from entering the endoplasmic reticulum stress state in patients, thus preventing the conversion of MCs into TAMCs, and thus controlling the progression of the cancer. In addition, gene regulatory network of C1 SRSF7+ MCs was revealed by scenic analysis, in which the most valuable key regulators were ATF4 and JUNB. ATF4 showed a dual role in iron death and cancer under endoplasmic reticulum stress(68), and under sustained stress conditions, ATF4 promotes apoptotic cell death induction. Characterizing the mechanisms that regulate ATF4-mediated transcription and its effects on cellular metabolism may identify novel targets for cancer therapy(45). As for JUNB, more and more studies have shown that it is involved in tumorigenesis by regulating cell proliferation, differentiation, senescence, and metastasis, and in particular, it affects the TME by transcriptionally promoting or repressing oncogenes in tumor cells or immune cells(69). Furthermore, previous mechanistic studies have shown that JUNB overexpression regulates the mitochondrial apoptosis pathway, mediating resistance to FasL and TRAIL-induced cell death, and thus tumor resistance to immunotherapy(70). This study of ours provides a theoretical basis for subsequent analysis of drug sensitivity and provides new insights into the development of innovative targeted therapeutics. To explore the interactions involving the C1 SRSF7+ mast cell subtype and other cell types, we employed CellChat communication pattern analysis. This approach helped reveal coordinated responses and interactions between different cell types in the context of their communication pathways. Different cell types can simultaneously activate common cell type-independent signaling pathways or different cell type-specific signaling transduction transduction pathways(39). Through CellChat analysis, we established the intercellular communication network between most cells, including tumor cells, fibroblasts, T-NK cells, and various subtypes of MCs, etc., as a way to characterize the relationship between the subtype of C1 SRSF7+ MCs and other cell types, and at the same time, we identified the three modes of outgoing, incoming and their corresponding signaling pathway expression. The C1 SRSF7+ MCs subtype belongs to mode 1 in the outgoing pathway, and its communication molecules, i.e., ligands, include ANNEXIN, PARs, CSF, ICAM, etc.; and it belongs to mode 2 in the incoming pathway, and its communication molecules, i.e., receptors, include BAFF, CLEC, ALCAM, SELPG, etc. It is also worth noting that tumor cells, which can be learned after our careful observation, belong to mode 2 on the outgoing and mode 1 on the incoming, echoing the subtype of C1 SRSF7+ MCs, which drew our attention. By targeting tumor cells and the C1 SRSF7+ MCs subtype for interactions analysis, we have identified the secretion of AREG ligands by a subtype of C1 SRSF7+ MCs in the EGF signaling pathway that act on the protein receptor EGFR on the membrane of the tumor cells. In previous studies, the EGFR family has been validated to play a key role in EGFR signaling through the activation of many important cellular processes, including cell division, growth, and differentiation. Playing a key role in mediating cell growth factor signaling(71), overexpression of EGFR signaling widely promotes tumor progression and leads to promotion of proliferation and inhibition of apoptosis(72). And cancer immunotherapies, particularly immune checkpoint blockade (ICB), have transformed oncology care over the past decade and significantly improved survival in a wide range of metastatic tumors. Based on significant treatment benefits, ICB therapy is approved by the FDA as monotherapy or in combination with other cancer therapies for cancers such as melanoma, breast cancer, renal cell carcinoma, head and neck squamous cell carcinoma, and lung cancer(73-77). However, the MCs-mediated pro-tumor axis AREG-EGFR in EC has not yet been mentioned. Therefore, our study provides new EC target therapeutic approaches and provides a scientific basis for the treatment and prognosis of EC. Meanwhile, to further investigate the role of EGFR in EC, we performed in vitro experiments using TE-10 and KYSE-30 cell lines. We observed that EGFR knockdown inhibited tumor cell activity, migration and proliferation, thereby suppressing tumor growth. However, previous studies have shown that epidermal growth factor receptor inhibitors (EGFRIs) produce a variety of dermatologic side effects in the majority of patients, and this targeted therapeutic regimen needs to be further refined(78). Given their role in promoting tumor growth and immune evasion, mast cells are considered potential therapeutic targets. Contemporary therapeutic strategies may include the use of mast cell stabilizers, mast cell mediator inhibitors, or blocking mast cell recruitment to tumor tissues and organs. Finally, we constructed a prognostic model to indicate that the higher the SMRS score, the worse the prognosis. Our study will direct attention to MCs in the progression of esophageal cancer, trigger attention to them, and promote researchers' understanding of the tumor microenvironment in esophageal cancer. Moreover, we discovered the existence of a specific subtype of MCs, C1 SRSF7+ MCs, which will facilitate the precise treatment of esophageal cancer. At the same time, we discovered the communication pathway between the tumor and our target MCs subtype. Although EFGR antagonists are still proved to have certain side effects, we believe that the development of targeted therapy will be further advanced in the future. However, this study still has some limitations. The relatively small sample size chosen is one aspect, and secondly, we only performed transcriptomics studies and in vitro experiments. The analysis of mast cell dynamics in EC using SCENIC and AUCell in our article is well-founded though and provides a detailed understanding of the regulatory networks that drive mast cell behavior. However, to draw more reliable conclusions, these findings must be validated by further experiments and compared across different cancer types. Next, we will integrate in vivo and in vitro experiments to provide a more comprehensive validation. In conclusion, the innovative features of our study lie in the use of high-resolution single-cell analysis technology, the construction of cell-cell interaction networks, the analysis of dynamic evolutionary trajectories, the identification of regulatory networks, and experimental verification, which provide new ideas for the targeted treatment of MCs in EC and new cell carriers for the development of EGFR targeted drugs. These will help to promote the in-depth development of the research on EC and provide new strategies for the disease.