A signature based on efferocytosis-related genes for the evaluation of prognosis and the tumour microenvironment in gastric cancer.

Introduction
Gastric cancer (GC) represents a significant malignancy that poses a serious threat to human health, with its mortality rate remaining alarmingly high and inflicting profound distress on countless families1. China is recognized as one of the regions with the highest incidence of GC, contributing approximately 40% of new cases and deaths globally each year2. Due to the insidious nature of early symptoms associated with GC, many patients are diagnosed at advanced stages. The predominant treatment strategy involves a combination of immunotherapy and chemotherapy; however, numerous patients fail to derive substantial benefits from these interventions, resulting in a persistently low 5-year survival rate among GC patients3,4. Consequently, there is an imperative need for comprehensive research into GC aimed at identifying effective therapeutic targets and developing robust prognostic assessment tools and treatment prediction methodologies to mitigate the mortality rate associated with this disease and provide renewed hope for affected individuals.
During the progression of cancer, various forms of cell death can occur in tumor cells due to mutations, hypoxia, and treatment with radiation therapy and chemotherapy, such as apoptosis, necrosis, ferroptosis, and pyroptosis5. Efferocytosis is an important process for removing dead cells6,7. In simple terms, the main functions of efferocytosis include two steps: “find me” and “eat or don’t eat me”. Apoptotic cells release “find me” signals, such as CX3CL1, which attract macrophages to the location of the dying cells. Apoptotic cells undergo changes in their cell membrane properties during the process of cell death, exposing phosphatidylserine (PS) to the cell surface as a “eat-me” signal. After migrating to the vicinity of apoptotic cells, macrophages can recognize and bind to “eat-me” signal molecules through their surface receptors, such as PS receptors8.
Research indicates that efferocytosis pro within tumors can play a unique regulatory role in promoting the formation of an immunosuppressive microenvironment9. For example, the cytokines produced by apoptotic cells during efferocytosis, such as sphingosine-1-phosphate (S1P) and transforming growth factor-beta (TGF-beta), can stimulate macrophages to polarize into anti-inflammatory and tumor-promoting M2 type10. Previous studies have shown that efferocytosis can facilitate metastasis, inhibiting efferocytosis or specifically eliminating the protein granulocyte colony-stimulating factor from macrophages can improve the function of CD8+T cells and reduce liver metastasis of pancreatic cancer11. At present, the study of efferocytosis in GC is still relatively limited, given the prominent role that efferocytosis plays, relevant research is highly necessary.
Single-cell sequencing (scRNA-seq) can evaluate the role of specific genes in specific cells, and more and more researchers are using scRNA-seq in combination with bulk RNA-seq to identify new biomarkers. The purpose of this study is to use a combined analysis approach to identify the role of efferocytosis-related genes (ERGs) in patient prognosis and immune infiltration, and to construct corresponding prediction models.

Materials and methods

Data acquisition
In this study, we obtained single-cell mRNA sequencing (scRNA-seq) data from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo), which includes 10 GC samples (GSE167297)12. The Bulk RNA sequencing dataset and clinical information of gastric adenocarcinoma (STAD) were downloaded from the TCGA database using UCSC Xena. A total of 343 GC samples and 31 para-cancerous samples were included. In addition, GSE26901 dataset was acquired from the GEO database and utilized as the validation cohort13. The data utilized in our study are publicly available or have been published in the scientific literature.

scRNA-seq data analysis and cell annotation
The processing of scRNA-seq data was conducted using the “Seurat” R package, version 4.4.0. This study implemented the following criteria to filter out low-quality cells1: Each gene must be expressed in a minimum of three cellular entities2; Cells with fewer than 200 detected genes were excluded3; Cells exhibiting more than 10% expression of mitochondrial genes were also excluded. The “scDblFinder” function was employed to refrain from doublets. The scRNA-seq data was subsequently standardized using the “NormalizeData” function. The “ElbowPlot” function is employed with variable genes as input to perform principal component analysis (PCA) for the identification of significant principal components (PCs). Additionally, the top PCs are selected as the statistically significant input of t-Distributed Stochastic Neighbor Embedding (t-SNE). Cell clustering and identification of cell types within the clusters were analyzed using the “FindClusters” function. We used the R package “SingleR” in conjunction with classical cell markers to annotate our scRNA-seq data14.

Screening ERGs in scRNA-seq data
We used “FindAllMarkers” function to analyze differential expression to get marker genes in each cell type. 137 ERGs were sourced from previously published literature15. Then, we obtained single-cell ERGs (SCERGs) by intersecting cell markers and ERGs. The expression of SCERGs was scored using the AUCell package and we categorized these cells into high-efferocytosis and low-efferocytosis groups based on the quartiles of SCERGs activity.

Intercellular communication and trajectory analyses
The “Monocle2” R package was utilized to unveil pseudo-time trajectories. We employed the DDRTree method for dimensionality reduction, projecting individual cells onto a spatial plane and organizing them into a trajectory with branch points. The dynamic expression heatmap was generated using the “plot-pseudo-time-heatmap” function.
Cell-cell communication plays a crucial role in various biological processes, achieved through the interaction between ligands and cell surface receptors. To investigate the interactions between different cell types, we have developed a tool called “CellChat”, which is an R package specifically designed for quantitative inference and analysis of cell communication networks based on single-cell transcriptome data16. By constructing networks of ligands, receptors, and their associated factors, CellChat simulates the complex and diverse communication relationships between cells. It utilizes gene expression profiles of ligands and receptors in different cell types to infer their interactions and reveals a rich diversity of receptor-ligand signaling mechanisms between two distinct cell types.

KEGG and GO enrichment analysis
We used “clusterProfile” package to conduct Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses17–19. The visualization of the results was conducted by using the “ggplot2” package in R. In addition, we also conducted enrichment analysis using the “BEST” online tool (https://rookieutopia.hiplot.com.cn/app_direct/BEST/)20. These analyses aim to reveal the key pathways that drive GC progression driven by ERGs.

Construction of a prognostic signature
Using the “FindAllMarkers” function to identify differential genes between the high-efferocytosis and low-efferocytosis groups. The intersection of specific genes in the high efferocytosis group and bulk RNA-seq differentially expressed genes (DEGs) was used to obtain differentially expressed single-cell ERGs (SCERDEGs). The initial step involved the utilization of univariate Cox analysis to identify prognostic risk genes associated with overall survival (OS) in GC patients. Subsequently, a LASSO analysis was employed to select significant genes, followed by multivariate Cox analysis for the construction of a prognostic signature. An optimal model was constructed with 4 SCERDEGs, output as gene risk score.

riskScore = ∑ni Coef(genei) * Expression(genei).

coef (genei), expression (genei) and n represent the coefficient of multivariate Cox regression, the expression of each gene and the number of genes, respectively.

Establishment and validation of the nomogram model
The risk scores and clinicopathological features were incorporated through the “RMS” package and “regplot” package to generate a nomogram. Time-dependent receiver operating characteristic curve (tdROC) and calibration curve were used to evaluate the accuracy and discriminative ability of the nomogram in the training and validation sets, respectively.

Immune infiltration analysis
Patients were divided into high-risk group and low-risk group according to the best cut-off value of gene risk score. The “ESTIMATE” package was utilized in the calculation StromalScore, ImmuneScore and EstimateScore based on gene expression data from tumor patients. The proportional infiltration levels of immune cells in two risk types were quantified and compared using “CIBERSORT” package.

Immune checkpoint analysis, immunotherapy response analysis and chemotherapy sensitivity analysis
The correlation between risk score and immune checkpoints was also examined to identify potential populations sensitive to immunotherapy. Tumor immune dysfunction and exclusion (TIDE) can predict the response of tumors to ICB, TIDE score of the sample were obtained from the TIDE website (http://tide.dfci.harvard.edu/) and compared between the two risk groups in the TCGA dataset. In addition, we used the “oncoPredict” package to analyze different drug sensitivity between high- and low- risk groups.

Statistics analysis
R software (version 4.4.1) was used for statistical analysis. The Kaplan–Meier method was used to draw the OS survival curve, and the log-rank test was used for statistical analysis. Two-tailed P value was less than 0.05 was considered as statistically significant.

Results

Identifying cell types through scRNA-seq analysis
As shown in Fig. 1, the flowchart illustrates the main research steps in this studay. We performed single-cell analysis on 10 tumor samples from GSE167297 and divided all cells into 16 clusters (0–15) (Fig. 2A). All the clusters of cells were divided into 7 cell subtypes, including T/NK cells, B cells, macrophages, plasma cells, endothelial cells, epithelial cells and stromal cells. (Fig. 2B). T/NK cells account for the highest proportion in GC tissues, followed by B cells. In Fig. 2C and D, we show the number of cells of each type and the proportion of each cell type in each sample. The violin plot shows the most representative marker genes for each cell type (Fig. 2E). Additionally, the top 5 marker genes for each cell type were also presented (Fig. 2F).

Screening of differentially expressed scergs in GC
Table S1 shows 137 ERGs obtained from published literature. We used “Seurat” package in R to identify differentially expressed genes for the 7 cell types relative to the other cell types (Table S1). In the end, we identified 20 SCERGs by intersecting the differential genes of the co-cultured cells and ERGs (Fig. 3A). Additionally, the heatmap shows the expression of 20 genes in each cell type (Fig. 3B). We can find that multiple SCERGs are highly expressed in macrophages.

To analyze its biological functions and associated pathways, we performed GO and KEGG analyses. In the GO analysis, we found that the gene set functions were mainly enriched in the functions of myeloid cell migration and activation, as well as antigen presentation, which are all related to anti-tumor immunity. In the KEGG analysis, the gene set was mainly enriched in chemokine pathway, immune cell receptor signaling pathway, JAK/STAT signaling pathway, and some infectious diseases, which are closely related to human immune function (Fig. 3C and D).

Scoring the cell type based on the expression of 20 scergs
Given the important function of efferocytosis in cancer development, we employed the “AUCell” method to evaluate each cell according to the expression levels of 20 SCERGs. The t-SNE plot illustrates the efferocytosis scores assigned to each individual cell (Fig. 4A). We compared the SCERG activity of different cell types and found that macrophages had higher SCERG activity, while B cells had the lowest activity (Fig. 4B). Then, we divided the cells into two groups based on the median of the cell burial score: high-efferocytosis and low-efferocytosis. It can be seen from Fig. 4C that high-efferocytosis group are mainly concentrated in the macrophages. The volcano plot primarily shows the expression differences between cancerous and adjacent tissues in GC patients in the TCGA database (Fig. 4D). Meanwhile, we compared the high efferocytosis group with the low efferocytosis group and identified 1432 genes. And then we intersected the differential genes with those from TCGA, ultimately obtaining 133 SCERDEGs for further analysis (Fig. 4E). The KEGG and GO enrichment analysis results showed that 133 SCERDEGs were enriched in molecular signaling pathways related to the immune microenvironment, such as the IL-17 signaling pathway and chemokine signaling pathway (Fig. 4F).

Analysis of the developmental trajectory and cellular communication of macrophages
Given the highest expression of SCERGs in macrophages, we performed dimensionality reduction clustering to extract them, ultimately dividing them into five cell subpopulations (MAC1-5) (Fig. 5A). Figure 5B corresponds to the expression of classical macrophage markers for each cell. Through further analysis, we believe that MAC1 and MAC4 have higher M2 polarization characteristics. We conducted pseudotime analysis to elucidate the differentiation transformation of macrophage subtypes. These cells formed a continuum with 7 cell states (Fig. 5C and D). During the pseudo-time, the number of MAC1, MAC3, MAC4, and MAC5 cells gradually increases, while the number of MAC2 cells gradually decreases (Fig. 5E). The bubble chart shows the expression of the 20 SCERGs in different macrophage subtypes (Fig. 5F). Among them, MAC2 has visible expression differences from other subtype cells. The heatmap shows the dynamic changes of 20 SCERGs expressions along pseudo-time (Fig. 5G). Genetic markers associated with efferocytosis may mediate the differentiation transformation between different subtypes of macrophages. Conducting cell chat analysis to explore the various interactions between macrophages and other cells in the TME. The number and strength of intracellular communication are shown in Fig. 5H and I. Overall, macrophages communicate frequently with endothelial cells and stromal cells, and the interaction between macrophages and B cells is the most intense.

Establishment and validation of a nomogram model
Perform prognostic analysis using 133 SCERDEGs, univariate Cox regression identified 13 genes that were associated with prognosis in GC patients, including CXCL3, VCAN, TSC22D3, CKLF, APOE, C4orf48, ZNF331, TUBA1A, ZFP36, MKNK2, DUSP1, GADD45B, and PER1 (Fig. 6A). Then, we used LASSO and multivariate Cox analysis to select 4 genes (CXCL3, VCAN, ZFP36 and MKNK2) and calculate gene risk scores (Fig. 6B/C/D). The Kaplan-Meier curve illustrated the survival curves for the 4 genes based on the optimal cutoff values, demonstrating significant survival differences between the high-expression and low-expression groups (Fig. 6E/F/G/H). Based on the optimal cut-off value of genetic risk score, patients were divided into high-risk group and low-risk group. The Kaplan-Meier curve results showed survival differences, which were statistically significant in both the TCGA database and the GSE26901 dataset (Fig. 6I/J). Meanwhile, the analysis of clinical features showed that the patients in the high-risk group had a higher TNM stage compared to those in the low-risk group, and the analysis of tumor differentiation grade also showed that the proportion of undifferentiated or poorly differentiated tumors was higher in the high-risk group (Fig. 6K/L). In the end, we successfully constructed a nomogram (Fig. 6M) for predicting survival time of GC patients by combining gene risk scores and various clinical pathological features, which showed good predictive performance in both external and internal validation sets, as evidenced by calibration curves(Fig. 6N/O) and tdROC (Fig. 6P/Q), validating the accuracy of the model.

Association of gene risk scores and immune infiltration
The significant difference in survival times between the two risk groups suggests a significant genetic heterogeneity between them. In the comparison of immune cell infiltration analysis, we found that patients in the high-risk group had higher abundance of M2 macrophages, and M2 macrophages are often closely related to the progression of cancer (Fig. 7A/B/C). Additionally, we found that there were fewer activated immune cells and more resting immune cells in the high-risk group. ESTIMATE calculates three TME scores, and the high-risk group exhibited a significantly higher ImmuneScore and StromalScore compared to the low-risk group; however, the corresponding tumor purity was lower (Fig. 7D).

The relationship between gene risk scores and immunotherapy prediction and drug sensitivity
We analyzed the expression differences of some immune checkpoints between the two risk groups, with most immune checkpoints expressing higher levels in the low-risk group, especially CD274 (Fig. 8A). Immune checkpoint inhibitors are an important treatment option for GC, and in the STAD sample from TCGA, the TIDE score was significantly higher in the high-risk scoring group, suggesting that patients with high-risk scores may benefit less from immunotherapy. (Fig. 8B). Additionally, we used the “oncoPredict” package to analyze drug sensitivity, obtaining the maximum half-maximal inhibitory concentration (IC50) for each sample corresponding to all chemotherapy drugs. Here, we mainly show the differences in sensitivity to several common chemotherapy drugs for GC, and the results show that the high-risk group has poorer sensitivity to drugs (Fig. 8C).

Discussion
GC is a highly aggressive and metastatic cancer that is a major threat to global health21. And, due to the heterogeneity of the tumor microenvironment, many GC patients do not benefit from chemotherapy or immunotherapy22, so identifying prognostic markers and biomarkers is crucial for accurately assessing the response to treatment. Recent studies have elucidated the role of efferocytosis in tumorigenesis and progression, particularly in the remodeling of the tumor microenvironment. For instance, it has been demonstrated that tumors can exploit efferocytosis to polarize macrophages towards an M2 phenotype, thereby impairing immune cell clearance functions and facilitating tumor-associated angiogenesis23. Regarding GC, research on efferocytosis remains limited; thus, our study aims to elucidate the role of efferocytosis in the development and immune infiltration of GC through a multi-omics approach.
ScRNA-seq is a cutting-edge technology that can reveal changes in transcriptomes at the single-cell level. It has unique advantages in exploring cellular heterogeneity and the cell types related to tumor development, progression, and metastasis24. In our study, we identified 7 cell types, including T/NK cells, B cells, macrophages, plasma cells, endothelial cells, epithelial cells and stromal cells, through analyzing 10 GC samples from the GSE167297 dataset. By conducting an intersection analysis between the selected marker genes of individual cells and EGRs, we identified 20 SCERGs. The AUCell analysis showed that the SCERGs score of macrophages were the highest. Meanwhile, the results of cellular communication analysis show that macrophages have close communication with T/NK cells and stromal cells. Macrophages are an important component of the tumor microenvironment and play a unique paradoxical role in the progression of tumors, while also serving as the main actor in the process of efferocytosis25. We subdivided the macrophages identified in this study into different subpopulations, and a total of five distinct macrophage types were identified. After analyzing the data, we believe that the MAC2 and MAC4 types exhibit significant M2 polarization characteristics, with high expression of CD163. During development, these two types of cells gradually increase. At the same time, the expression of some SCERGs also increases during development, which may be closely related to M2 polarization.
Subsequently, we identified a total of 133 SCERDEGs for analysis by intersecting the differential genes from the high-efferocytosis group with those from TCGA. Through a combined analysis utilizing LASSO and Cox regression, we pinpointed 4 independent prognostic genes: CXCL3, VCAN, ZFP36, and MKNK2. CXLC3 is a member of the chemokine family, and studies have shown that it is mainly released by inflammatory cells such as macrophages. CXLC3 is expressed at higher levels in tumor tissue and can affect the tumor microenvironment, such as tumor-associated fibroblasts, and promote tumor metastasis26,27. However, the relationship between the expression level of CXCL3 and tumor prognosis has been contradictory in previous studies. According to the results of this study, high expression of CXCL3 suggests a better prognosis, but the specific mechanism of prognostic influence deserves further exploration27,28. Several studies have already demonstrated that VCAN and ZFP36 are associated with poor prognosis in patients. The results of this study are consistent with those of previous studies29–32. Additionally, research has shown that MKNK2 exhibits an inverse relationship with M2 macrophage markers CD163 and CD206. These data, combined with findings from previous studies, demonstrate the role of MKNK2 in regulating macrophage pre-tumorigenic phenotypes33. Utilizing the optimal cut-off value of the gene risk score across all patients, GC patients were stratified into two groups, effectively distinguishing survival differences between them and addressing some limitations of TNM staging in clinical practice. Furthermore, clinically relevant studies have demonstrated that the high-risk group is significantly associated with advanced TNM staging and poor differentiation grade among patients. Ultimately, we constructed a prognostic nomogram for GC patients based on genetic risk scores and clinical characteristics. In the training and validation cohorts, the calibration curve and ROC analysis confirmed the accuracy and reliability of the model. We believe that utilizing this model can enable clinicians to more accurately assess patient outcomes and make better clinical judgments.
This study found that the ERGs-related risk score is closely related to the immune microenvironment of patients, among which the most closely related is macrophages. We found that the ERGs-related risk score is positively correlated with M2 macrophage infiltration. Multiple studies have confirmed that M2 macrophages can impair anti-tumor immunity and promote tumor progression, which may be the reason why the risk score we have constructed is significantly correlated with prognosis34,35. Meanwhile, KEGG analysis revealed that chemokine signaling pathway and IL-17 signaling pathway were highly enriched, and these two are considered to be closely related to the remodeling of the tumor microenvironment36,37. And relevant studies have confirmed that IL-17 is associated with M2 polarization of macrophages38. In addition, the study also found that the number of CD8 positive T cells and activated CD4 memory T cells decreased as the risk score increased, with these cells being considered the mainstay of anti-tumor immunity39,40. In conclusion, we speculate that efferocytosis may influence the polarization of macrophages in GC patients through the IL-17 and chemokine pathways, mainly promoting M2 polarization, thereby affecting the patients’ anti-tumor immune response.
Moreover, this study also found that the ERGs-related risk score was associated with the treatment response and prognosis of patients. Immunotherapy is one of the principal treatment modalities for GC patients. Nevertheless, numerous patients fail to reap benefits from immunotherapy41. Hence, exploring predictive markers related to the therapeutic efficacy of immunotherapy is of great significance. In this study, we predicted the response of tumors to immunotherapy and chemotherapy based on ERGs-related risk score, and we can find that the TIDE score is higher in the high-risk group, which may benefit less from immunotherapy. Poor response to immunotherapy may be related to the immunosuppressive microenvironment caused by efferocytosis. Meanwhile, the results of analyzing the sensitivity of common GC chemotherapy drugs show that the high-risk group has lower sensitivity to drugs. Therefore, by evaluating the ERGs-related scores, it will be beneficial to predict the patient’s treatment response. In addition, we found that the ERGs-related risk score was associated with the clinicopathological characteristics of patients. Patients in the high-risk group had a higher stage and poorer differentiation. Ultimately, we constructed a prognostic model for GC patients based on ERGs-related risk score and clinical features. This model can accurately predict the prognosis of GC patients and is conducive to early identification of high-risk patients by doctors.
Currently, the process of efferocytosis in GC is still poorly understood, and our study aims to provide new insights from this angle. However, the study also has some limitations. Firstly, the study mainly analyzes data that has already been published, and it is necessary to validate our findings in real-world cohorts to ensure the robustness of ERGs as a prognostic marker for GC. Secondly, the mechanism by which efferocytosis affects tumors needs to be verified through more experiments. Despite these limitations, the study still provides some compelling evidence for the study of efferocytosis in GC.

Conclusion
In summary, we constructed a 4 gene signature related to efferocytosis by integrating scRNA-seq and bulk RNA-seq data, which can effectively predict GC survival time and treatment response. At the same time, we found that ERGs are closely related to the infiltration of immune cells such as M2 macrophages. This study provides our own insights into the efferocytosis function in GC and targeting ERGs such as VCAN may be a potential treatment for GC.

Electronic supplementary material
Below is the link to the electronic supplementary material.