Single-cell data revealed CD14-type and FCGR3A-type macrophages and relevant prognostic factors for predicting immunotherapy and prognosis in stomach adenocarcinoma

Introduction

Stomach adenocarcinoma (STAD) threatens the lives and health of at least one million individuals every year and is an urgent public health problem in the world (Smyth et al., 2020). Similarly, China had a huge number of STAD patients (Abbas et al., 2018; Yang et al., 2021). Poor lifestyle habits such as high oil and salt, virus infection or smoking and alcohol abuse in some regions were the leading causative factors of STAD (Jafari-Sales et al., 2022; Yang et al., 2021). With limited cancer symptoms and unpopular awareness of early screening, many STAD patients progress to advanced stages before being detected (Hoft, Noto & DiPaolo, 2021; Li et al., 2022). Endoscopy is recognized as the best early screening approach until now and has improved the survival advantage of STAD patients by at least 30% (Machlowska et al., 2020). However, lack of awareness of early screening, existing malignant metastasis and treatment resistance (Chen et al., 2022; Choi et al., 2022), STAD remained the most lethal malignancy with a 5-year survival rate of less than 20% (Ilic & Ilic, 2022). Therefore, innovative mining of prognostic and therapeutic biomarkers in STAD was urgent.

Up to data, massive genes models such as oxidative stress and metabolism-related (Dong et al., 2023), memory B-cell-corrected miRNA (Liu et al., 2023), complement system-related (Tong et al., 2022) as well as ferroptosis-relevant (Zhu et al., 2022) signature were established hoping to predict the potential risks that patients may face, and help doctors develop more personalized treatment plans, for improving patients’ treatment effectiveness and survival rate. However, the heterogeneity of STAD was a huge challenge for treatment, and an in-depth understanding of the mechanisms involved was the key to overcome the treatment tolerance of STAD (Katzenelenbogen et al., 2020). The emerging single-cell sequencing technology (scRNA-seq) sequenced individual cells to precisely obtain the genetic features and clonal diversity in cells, which allowed the discrimination of tumor heterogeneity by cell clustering (Yu et al., 2021). For example, the heterogeneous immune landscape between LUAD and squamous carcinoma were discovered by scRNA seq approach (Wang et al., 2022), indicating different immunotherapy methods can be adopted for these two kinds of lung cancers. In LUAD patients with EGFR mutations, the heterogeneous tumors and immune cell populations in early period were also discovered as a promising diagnostic method for LUAD (He et al., 2021). Shen, Chen & Gao (2023) constructed a mesenchymal stem cell-relevant model for assessing prognosis and drug sensitivity in STAD patients via scRNA-seq approach. Similarly, through integrating scRNA-seq with bulk data, CXCR4 on tumor-infiltrating B cells showed inferior survival outcomes and was considered as a therapeutic targets for STAD (Su et al., 2023). However, macrophage relevant study based on scRNA-seq method in LUAD is rare, only one M2 macrophage-related prognostic and immunotherapeutic signature in ovarian metastasis of STAD was published (Gao et al., 2023). Therefore, precise identification of macrophage cell types and deeply explore its actions in STAD are meaningful.

Hereby, we portrayed the cell types in STAD as well as the landscape through scRNA-seq data in an attempt to reveal the cell types that threaten STAD prognosis. We found that macrophage subgroup score was higher than other immune cells. Subsequently, using bulk sequencing data we attempted to establish an effective macrophage related prognostic assessment system to precisely assess prognostic and therapeutic features for STAD patients.

Materials and Methods

Results

Cell types in STAD

After screening and processing, sequencing data of 22,258 cells in 14 samples from GSE167297 were included in this study (Figs. S1A and S1B). We selected the top 3,000 HVGs for PCA and the top 50% of PCs for subsequent clustering analysis (Fig. S1C). The batch effect between samples was eliminated before annotation (Figs. S1D and S1E). According to UMAP method clustering, cells in all samples were clearly divided into 26 cell clusters at resolution = 0.8 (Fig. 1A). The 26 cell clusters were defined as 11 cell types according to marker genes expression characteristics (Fig. 1B and Table 2). The marker genes specifically expressed in the 11 cell types were shown by different visualization forms (Figs. 1C, S2 and S3). The number of cells in the samples after quality check and the proportion of annotated cell types are shown in Figs. S1F and S1G. The proportions of B cells, T cells, and Macrophages cell types in the samples and the distribution of cells in clusters were illustrated in Figs. 1D and 1E, with a smaller proportion of Macrophages. Macrophages showed higher AUCell scores, indicating higher cellular activity in Macrophages (Figs. 1F and 1G). These results indicated a greater impact on STAD despite the lower percentage of Macrophages.

Table 2:

Annotation information of cell types.

Cell type	Marker gene	Seurat clusters
B cells	CD79A, MS4A1	1,3,20
T cells	CD3D, CD3E	0,2,4,8,14,21,24
Macrophages	CD14, FCGR3A	5
Plasma cells	MZB1, SDC1	6,10,11,24
Epithelial cells	EPCAM, KRT19	7,17,23,25
Endothelial cells	TM4SF1	9,11,12,17,22,25
Fibroblasts	COL1A1, FGF7	12
Mast cells	TPSAB1, CPA3	16
Multipotent progenitor cells	UBE2S, H2AFX	12,17,18,20,21,25
Multilymphoid progenitor cells	COTL1, TMSB10	0,1,2,3,4,5,7,8,11,12,13,16,17,18,19,20,21,24,25
Mesenchymal cells	ACTA2, AKAP12	12,22

DOI: 10.7717/peerj.16776/table-2

Macrophage-related gene modules

Based on six immune cell (B cell, T cell CD4+, T cell CD8+, neutrophil, macrophage, myeloid dendritic cell) immune scores obtained in TIMER from STAD patients, they were divided into high and low immune score groups based on median values. We found that survival differences were exhibited only in the macrophage groups (Fig. 2A). Macrophage affected survival in STAD patients, and higher macrophage activity was also shown in single-cell data. Macrophage might be a critical factor in regulating the prognosis of STAD. Further, we identified gene modules highly associated with macrophage by WGCNA. The network was biologically significant when the soft threshold β = 8 and the correlation coefficient = 0.85 (Fig. 2B). According to the dynamic shear number procedure, six gene modules were identified (Figs. 2C and 2D). A Pearson correlation analysis was performed with the six gene modules using age, gender, T stage, N stage, M stage, stage, grade, and macrophage as trait vectors. There was a remarkable positive correlation between the pink module and macrophage (r = 0.69, p < 1e−5), thus, we defined pink as a macrophage-related gene module (Fig. 2E). Pathway enrichment analysis for genes within the pink module exhibited that these genes were mainly enriched in pathways in cancer, PI3K-Akt signaling pathway, proteoglycans in cancer, regulation of multicellular organismal process as well as regulation of localization. These genes had molecular functions in signaling receptor binding and extracellular matrix structural constituent, and participated in cellular component such as extracellular region part and extracellular space (Fig. 2F). It indicated that genes within the pink module actively regulated the intercellular communication exchange in STAD.

RiskScore for predicting STAD prognosis

In single cell analysis, 670 marker genes were significantly differentially expressed in Macrophage cell types versus other cells. It was extracted with 86 duplicated genes within the pink module to construct the LASSO COX model (Fig. 3A). The model was observed to be optimal at λ = 0.0389 under 5-fold cross-validation (Figs. 3B and 3C). RGS2, GNAI2, ANXA5, MARCKS, CD36, NRP1, and PDE4A were the prognostic factors of STAD, among which, GNAI2 and PDE4A were protective indicators, while RGS2, ANXA5, MARCKS, CD36 and NRP1 were risk indicators (Fig. 3D). The prognostic assessment value of RiskScore was judged based on K-M survival analysis and ROC curve in the training set (TCGA_GTEx-STAD) and validation set (GSE26901, GSE15459). We found that the AUC values in 1 year among three dataset were 0.68, 0.67 and 0.72 respectively. The AUC values in 3 years among three dataset were 0.70, 0.71 and 0.71 respectively. The AUC values in 5 years in two validation sets were 0.71 and 0.68 respectively. In short, all AUC values were higher than 0.6 in all three datasets (Figs. 3E, 3G and 3I). We also observed that low RiskScore showed prognostic advantages in all three datasets (Figs. 3F, 3H and 3J). As a result, RiskScore was a trustworthy prognostic system.

Expression landscape and mutation landscape of prognostic factors

Expression analysis was conducted in order to analyze the expression levels of prognostic factors and protein products in the RiskScore system. First, the protein product expression levels of RGS2, GNAI2, ANXA5, MARCKS, CD36, NRP1, and PDE4A were queried in the HPA database. The expression of the protein encoded by MARCKS gene was lower in STAD tissues than in normal gastric tissues (Fig. 4A). The mRNA levels of RGS2, PDE4A, and CD36 were elevated in STAD tissues, and the mRNA levels of GNAI2, MARCKS, and NRP1 were elevated in normal tissues (Fig. 4B). The expression of seven prognostic factors was also expressed with tissue specificity. RGS2 was specifically expressed in Cytosol, ANXA5 was specifically expressed in Nuclear membrane, and CD36 was specifically expressed in Golgi apparatus (Fig. 4C). Besides, we found a significant copy number deletion of GNAI2 (Fig. 4D).

Nomogram based on the composition of RiskScore, stage, and age

We discussed the proportion of patients with different clinical factors in the RiskScore groups, with significant differences between patients in the T stage and grade subgroups (Fig. 5A). As the clinical factors increased, the RiskScore of patients also showed a tendency to increase (Fig. 5B). Univariate and multivariate COX analyses based on clinical factors and RiskScore showed that RiskScore, stage, and age were independent prognostic factors for STAD, and all were risk factors (Hazard Ratio > 1, p < 0.05) (Figs. 5C and 5D). Therefore, based on RiskScore, stage, and age, we plotted the nomogram used to predict 1-year, 3-year, and 5-year survival of STAD patients (Fig. 5E). The calibration curve showed a good fit between the prediction curve and the actual observed 1-year, 3-year, and 5-year survival curves; therefore, the nomogram was a reliable prognostic prediction tool for STAD (Fig. 5F). Decision curve also further confirmed the strong robustness of nomogram and RiskScore, both of which had higher benefit than the extreme curves (Fig. 5G). Compared to the prognostic value of other clinical factors, the AUC values of nomogram and RiskScore were remarkably superior, showing strong prognostic predictive value (Fig. 5H).

Biological pathways and immune landscape of differences in RiskScore groups

This study also discussed the biological pathways regulated by DEGs in the RiskScore groups. We found 506 up-regulated DEGs in the high RiskScore group (Fig. 6A), while the number of down-regulated DEGs was minor, so the Over-Representation Analysis was selected. The ORA results showed that the upregulated DEGs were mainly involved in the processes of formation of vascular morphology, extracellular matrix components, and glycosaminoglycan binding (Fig. 6B). The presence of 50 HALLMARK pathways in the RiskScore subgroup was identified in ssGSEA, and we assessed the Spearman correlation between them and RiskScore (Fig. 6C). The GSEA results showed that HALLMARK_TNFA_SIGNALING_VIA_NFKB, HALLMARK_HYPOXIA, HALLMARK_IL6_JAK_STAT3_SIGNALING, HALLMARK_APOPTOSIS, HALLMARK_APOPTOSIS, and HALLMARK_ANDROGEN_RESPONSE, HALLMARK_MYOGENESIS, HALLMARK_INTERFERON_GAMMA_RESPONSE, HALLMARK_APICAL_JUNCTION, HALLMARK_COMPLEMENT, HALLMARK_EPITHELEMENT HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION were enabled in high risk group (Fig. 6D). In addition, we assessed the immune landscape of patients in the RiskScore subgroup. ImmuneScore, StromalScore were higher and exhibited higher T cell CD8+, macrophage, myeloid dendritic cell scores in the high RiskScore group (Figs. 6E and 6F). The RiskScore showed a positive correlation with the TIDE score, indicating that the high risk group was more prone to immune escape and less prone to benefit from immunotherapy (Fig. 6G).

Experimental validation of the reliability of RiskScore models

To validate the RiskScore models, we examined the expression of RGS2, RNAI2, ANXA5, MARCKS, PDE4A, CD36, and NRP1 in MKN7, MKN28, and GES1 cell lines using qRT-PCR (Figs. 7A–7G). The results of qRT-PCR matched the expression levels in RNA-seq data analysis. WB was conducted to detect the expression of MARCKS, which contributed most to the RiskScore models. The results showed that MARCKS expression was elevated in gastric cancer cell lines MKN7 and MKN28 (Fig. 7H). After inhibition of MARCKS in MKN7 and MKN28, decreased cellular activity in MKN7 and MKN28 cells were observed by cell cycle (Figs. 7I and 7J) and CCK8 assays (Figs. 7K and 7L).

Discussion

STAD continues to be one of the most common causes of cancer death worldwide (Sung et al., 2021). Proven targeted therapeutic markers still lack therapeutic specificity; therefore, there is an urgent need to find novel and effective biomarkers (Choi et al., 2022).

In this study, macrophages were found to be the most active cell type in STAD and specifically express CD14 and FCGR3A. Tumor-associated macrophages (TAMs) were important anti-cancer factors in the anti-tumor immune response. There were numerous biological functions in anti-cancer or pro-cancer responses. On the one hand, TAMs were major regulators of angiogenesis in diffuse STAD, where a high density of M2-type macrophages was the main trigger of the immunosuppressive phenomenon in STAD. On the other hand, TAMs also positively influenced the treatment resistance in STAD patients (Balkwill & Mantovani, 2012; Chung & Lim, 2014; Ge et al., 2018; Joyce & Pollard, 2009). Macrophages were indeed factors affecting the prognosis or treatment of STAD (Wang et al., 2022). Macrophage marker gene CD14, could promote M2-type macrophage abundance (Shi et al., 2021), and FCGR3A+ macrophages might attenuate cytokinesis or the ability to engulf apoptotic cells (Lescoat, Lecureur & Varga, 2021). M2-type macrophages stimulated tumor growth by promoting tumor immunosuppression (Pan et al., 2020). Here, we revealed a specific type of macrophages in STAD, revealing the heterogeneity of macrophages. And high level of macrophage score affected the poor prognosis of STAD patients.

Macrophage-related genes were mainly enriched in ECM-receptor interaction. Another remarkable feature of cancer patients was the high level of ECM accumulation, and ECM accumulation also worsened the prognosis (Izzi, Davis & Naba, 2020). Based on macrophage-related genes, we established a RiskScore prognostic system for LUAD patients. RGS2, GNAI2, ANXA5, MARCKS, CD36, NRP1, and PDE4A were prognostic factors for STAD. Up-regulation of RGS2 expression resulted in a biomarker of delayed value-added and poor prognosis in non-small cell lung cancer, relieving the dormancy of slow-cycling/dormant cancer cells after surgery/chemotherapy and causing tumor recurrence (Cho et al., 2021). Coherently, high level of RGS2 was also found in gastric cancer cells, and was closely connected with TIDE, and CD8+ T-cell infiltration in other cancer types (Yang et al., 2022), indicating that RGS2 could studied further as a new target for immunotherapy on gastric cancer. Yu et al. (2022) showed that GNAI2 facilitated STAD cell growth and migration by promoting PI3K/AKT pathway bioactivity. ANXA5 served as a bridge between the innate and adaptive immune systems, and the ANXA5-phosphatidylserine axis targeted the tumor environment (TME) and delivered chemotherapeutic agents to the TME to kill tumor cells therein (Woodward, Faria & Harrison, 2022). Data from scRNA in the ascites of STAD patients indicated that MARCKS marked high plasticity and severe survival in STAD (Huang et al., 2023) MARCKS could trigger NF-κB pathway in smoke-relevant lung cancer (Liu et al., 2021), and MARCKS-targeting drug (BIO-11006) has displayed preliminary success in clinical trials of lung diseases (Yadav et al., 2023). Here, our data further showed that inhibition of MARCKS activity markedly reduced the viability of gastric cancer cell lines. Excessive accumulation of fatty acid molecules in metastatic STAD induced CD63 expression amplifying O-GlcNAcylation signaling levels and stimulating the value-added potential of STAD cells (Jiang et al., 2019). NRP1+ lung cancer cells has been ascertained to have tumor-initiating properties (Jimenez-Hernandez et al., 2018), and the blockade of NRP1 synergises with anti-PD-1 could strengthen CD8+T-cell proliferation, cytotoxicity as well as tumor inhibition (Leclerc et al., 2019). Thus, it was recognized as a novel target for CAR-T immunotherapy in STAD (Bebnowska et al., 2020), demonstrating that it is a risk indicator for LUAD. PDE4A was a trigger for metastasis in hepatocellular carcinoma, and inhibition of PDE4A function had the potential to suppress tumor metastasis in patients with hepatocellular carcinoma (Peng et al., 2018). Prognostic factors in STAD were mainly associated with tumor cell metastasis, immune cell response, and post-treatment recurrence, and they comprised a RiskScore that took into account important factors in tumor progression. Furthermore, we plotted nomogram, ROC curve and Decision curve also showed that nomogram and RiskScore exhibited excellent robustness and prognostic advantages, and predicting STAD prognosis by our system is of potential clinical value.

Finally, we carried out qRT-PCR analysis of the expression levels of RGS2, GNAI2, ANXA5, MARCKS, CD36, NRP1, and PDE4A in gastric cancer cells. The results were shown to be compatible with the expression levels of 7 prognostic factors in RNA-seq. The results of qRT-PCR indicated that the STAD prognostic factors in the model were reliable. Besides, MARCKS was a prognostic factor that inhibited the viability of gastric cancer cells. Studies demonstrated that MARCKS was a marker of poor prognosis in highly plastic gastric cancer (Huang et al., 2023). However, the functional studies of MARCKS remained poorly characterized, and subsequent focus on its function is still needed.

Our study aimed to identify the specific cell types that affect the prognosis of STAD patients and attempted to find novel cellular targets. An attempt was also made to establish a novel prognostic system, validating their prognostic potential in multiple data. The limitations of present research were as below: although we finally disclosed a single cell atlas and established a RiskScore signature based on microphage-relevant genes, more database and clinical samples are needed to test the accuracy of present findings. Moreover, the mechanism of how model feature genes regulate cancer needs to be further verified in in vivo/in vitro assays to enhance their clinical values.

Supplemental Information