Research Paper Volume 16, Issue 12 pp 10252—10270
Characterization of tumor endothelial cells (TEC) in gastric cancer and development of a TEC-based risk signature using single-cell RNA-seq and bulk RNA-seq data
- 1 Department of Gastrointestinal Surgery, Zhu Cheng People’s Hospital, Weifang, China
Received: December 7, 2023 Accepted: April 22, 2024 Published: June 12, 2024
https://doi.org/10.18632/aging.205928How to Cite
Copyright: © 2024 Fan et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Abstract
Background: Tumor endothelial cells (TECs) are essential participants in tumorigenesis. This study is focused on elucidating the TEC traits in gastric cancer (GC) and constructing a prognostic risk model to predict the clinical outcome of GC patients.
Methods: Single-cell RNA sequencing (scRNA-seq) data were obtained from the GEO database. Using specific markers, the Seurat R package aided in processing scRNA-seq data and identifying TEC clusters. Based on TEC cluster-associated genes identified by Pearson correlation analysis, TEC-related prognostic genes were screened by lasso-Cox regression analysis, thereby constructing a risk signature. A nomogram was created by combining the risk signature with clinicopathological features.
Results: Based on the scRNA-seq data, 5 TEC clusters were discovered in GC, with 3 of them showing prognostic associations in GC. A total of 163 genes were pinpointed among 3302 DEGs as significantly linked to TEC clusters, leading to the formulation of a risk signature comprising 8 genes. Furthermore, there was a notable correlation between the risk signature and the immune cell infiltration. Multivariate analysis findings indicated that the risk signature served as an independent prognostic factor for GC. Moreover, its efficacy in forecasting immune response was validated.
Conclusion: TEC-based risk model is highly effective in predicting the survival outcomes of GC patients and can forecast the immune response. Targeting TECs may significantly inhibit tumor progression and enhance the efficacy of immunotherapy.
Introduction
Gastric cancer (GC) is a common worldwide cancer that presents a major risk to health [1, 2]. Its intricate evolution involves the interplay of dietary factors, host genes, Helicobacter pylori infection, and environmental elements [3–5]. Because of atypical symptoms, gastric cancer at an early stage is frequently not identified until it has progressed to a more advanced stage. This leads to a poor prognosis characterized by a high likelihood of local recurrence and distant metastasis [6, 7].
Tumor endothelial cells (TECs) are crucial components of the tumor microenvironment (TME) and play a pivotal role in promoting tumorigenesis and metastasis [8–13]. Throughout tumor development, TECs contribute not only to the formation of new blood vessels but also impact the biological behavior of tumor cells by secreting a variety of bioactive molecules, including growth factors, cytokines, and enzymes [12, 14, 15]. Additionally, TECs possess the capability to modulate immune responses. They can produce immunosuppressive molecules such as TGF-β and IL-10, which dampen the activity of immune cells, thereby facilitating tumor cells to evade immune detection and destruction [14].
Additionally, TECs are also able to influence the invasive and metastatic ability of tumor cells through direct interaction with them [16]. For instance, by secreting enzymes like matrix metalloproteinases (MMPs) that aid in tumor cell extravasation and basement membrane penetration, TECs assist tumor cells in traversing the vessel wall, entering the bloodstream, and establishing distant metastases [14, 17–19]. Of note, TEC has the capacity to interact with a variety of cells and molecules present in the TME [14]. Other cell types within the tumor microenvironment (TME), such as cancer-associated fibroblasts (CAFs) and immune cells, synergistically cooperate with TEC in establishing a favorable microenvironment for tumor progression and metastasis by secreting signaling molecules or engaging in direct cell-to-cell interactions [20, 21].
It has been reported that distinguished from normal endothelial cells (NECs), TECs exhibit chromosomal instability, altered gene expression, enhanced proliferative capacities, and resistance to anti-angiogenic drugs [22, 23]. Targeting TECs can inhibit tumor progression and prolong patient survival by blocking the formation of tumor blood vessels [16, 22, 24]. Based on the scRNA-seq analysis, the study by Yin et al. revealed a TEC cluster in GC tissues that exclusively expressed IGFBP5 and displayed malignant characteristics [23]. Patients with TECs overexpressing IGFBP5 and IGFBP3 demonstrated significantly lower overall survival (OS). While numerous studies have been undertaken, a comprehensive understanding of the systematic characteristics of TECs and their correlation with GC prognosis and response to immunotherapy remains elusive.
In this research, we acquired scRNA-seq and transcriptome data, identifying distinct TEC clusters and developing a risk signature based on TECs. We assessed the clinical significance of the risk signature and examined the relationship between the risk signature and TME, immune response. Subsequently, we created a new nomogram that combines the risk signature with clinicopathological characteristics, to enhance its performance in predicting the clinical outcome of GC. Our study provides new insights into the molecular mechanisms underlying the occurrence and progression of GC, allowing for more individualized therapies.
Methods
Data acquisition and processing
From the public online website (https://dna-discovery.stanford.edu/research/datasets/) [25], we acquired scRNA-seq data of GC tissues, which included 9 GC samples, 10 normal samples, 2 samples of peripheral blood mononuclear cells (PBMC), and 1 sample of metaplasia. Our focus was specifically on the 9 GC tissue samples and their corresponding normal samples.
We used strict standards for scRNA-seq data, ensuring that each gene was expressed in a minimum of 3 cells and that each cell expressed at least 250 genes. To ensure quality control (QC) with the Seurat R package [26], cells with mitochondrial gene expression percentages higher than 15% were eliminated. Additional filtering based on QC criteria (nFeature RNA ≥ 200 and nFeature RNA ≤ 5000 and nCount_RNA ≥ 500) was applied. Normalization was executed using Seurat R package’s “NormalizeData” function, resulting in 37,440 cells. We adjusted the number of principal components (PCs) to 30 for the generation of cell clusters. The identification of cell clusters was performed using the FindClusters function with a resolution set to 1, and the results were visualized through UMAP method.
The Cancer Genome Atlas (TCGA) provided RNA-seq data, single-nucleotide variant (SNV), copy number variant (CNV) data, and clinical information for GC. After additional curation, the transcriptome data yielded a total of 345 tumor samples and 32 para-cancerous samples. Furthermore, for further validation, GSE62254 (300 GC samples) and GSE15459 (192 GC samples) cohorts were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Moreover, ten pathways related to cancer were acquired from the literature [27].
Definition of TEC
The Seurat R package enabled a thorough re-examination of scRNA-seq data to characterize the signature of TECs. The function ‘FindIntegrationAnchors’ was used to handle batch effects in 18 samples. The uniform manifold approximation and projection method with 30 principal components (PCs) and a resolution of 0.25 was employed for non-linear dimensional reduction. The process of clustering was executed by utilizing the functions ‘FindNeighbors’ and ‘FindClusters’. The function ‘RunUMAP’ was utilized to apply UMAP dimensionality reduction. Endothelial cells were characterized by 8 marker genes, namely VWF, ENG, PLVAP, MCAM, PECAM1, CLDN5, SELE, and SELP. Furthermore, UMAP dimensionality reduction was applied to the endothelial cell clusters. Marker genes for each TEC cluster were determined using the FindAllMarkers function. The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis on TEC cluster marker genes was conducted using the clusterProfiler R package [28]. Additionally, the CNV characteristics within the TEC clusters were analyzed using the Copykat R package to discern between malignant and non-malignant cells in each sample [29].
TEC-related gene identification
The limma R package was used to identify DEGs between tumor and normal tissue, with an adjusted p-value of less than 0.05 and an absolute log2 (fold change) greater than 1 [30]. Correlations between DEGs and TEC clusters were evaluated, with a focus on the identification of crucial TEC-associated genes. The genes associated with prognosis were further determined through univariate Cox regression analysis in the survival package with a significance level of p < 0.05. In order to reduce the number of genes, a lasso-Cox regression analysis was conducted, followed by a multivariate Cox regression analysis using a stepwise regression approach. Based on the outcomes of the multivariate Cox model, a risk signature was developed utilizing the formula: risk score = Σβi × Expi, where i represents the gene in the risk signature, expi denotes the expression of gene i, and βi indicates the coefficients of gene i in the multivariate Cox model. Subsequent to zero-mean normalization, the patients were stratified into high- and low-risk groups. The predictive performance of the risk signature was evaluated through receiver operating characteristic curve (ROC) analysis using the timeROC R package.
Immune landscape analysis
The CIBERSORT algorithm evaluated the proportions of 22 distinct immune cell types within the context of GC [31]. By computing StromalScores, ImmuneScores, and ESTIMATEScores, the ESTIMATE algorithm facilitated the investigation of the connection between the risk genes and the TME [32].
The ability to respond to immune checkpoint inhibitors
We acquired transcriptomic and clinical information from GC patients who received treatment with immunotherapy from the IMvigor210 and GSE78220 cohorts [33, 34]. The assessment was conducted to determine the potential significance of the risk signature in predicting the effectiveness of immune checkpoint inhibitors (ICIs).
Statistical analysis
R software (Version 4.2.1) was utilized for all the analyses. A p-value less than 0.05 was deemed to be statistically significant.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article.
Results
Screening of TECs in scRNA-seq samples
Figure 1 depicts the study’s flow chart. At first, a total of 37,440 cells were obtained from scRNA-seq data after an initial QC. After log-normalization and dimensionality reduction, a total of 36 subgroups were observed, uncovering a single TEC population distinguished by 8 marker genes: VWF, ENG, PLAVP, MCAM, PECAM1, CLDN5, SELE, and SELP (Supplementary Figure 1A, 1B). Afterward, TEC group cells went through further clustering and dimensionality reduction. The TEC populations were divided into five clusters using the identical clustering algorithm (Supplementary Figure 1C, 1D). In Figure 2A, the UAMP plot shows the distribution of 18 samples. As a result, five TEC clusters were finally generated and used for subsequent analysis (Figure 2B). 601 DEGs were identified among the 5 TEC clusters. Figure 2C displays the expression of the top 5 DEGs, which act as marker genes for the TEC clusters. Figure 2D indicates the proportion of the five TEC clusters in each sample. The analysis of KEGG showed enrichment in different pathways, including mineral absorption, ferroptosis, necroptosis, ribosome, focal adhesion, rap 1 signaling pathway, regulation of actin cytoskeleton, platelet activation, ECM-receptor interaction, adherens junction, and cell adhesion molecules (Figure 2E). Moreover, five TEC clusters were identified to comprise 839 tumor cells and 1137 normal cells according to CNV traits (Figure 2F).
Figure 1. The flow chart of this study.
Figure 2. The identification of TEC clusters based on scRNA seq data of GC patients. (A) UMAP plot of the distribution of 18 samples. (B) UMAP plot of the distribution of five TEC clusters after clustering. (C) Dot plot of the top 5 marker gene expression of TEC clusters. (D) The proportion of the five TEC clusters in tumor samples and normal samples. (E) KEGG enrichment analysis of 5 TEC clusters. (F) UMAP distribution map of malignant and non-malignant cells predicted by Copykat package.
Expression of cancer-related pathways in TEC
To explore the relationship between TECs and tumorigenesis, we calculated GSVA scores for ten oncogenic signaling pathways across distinct TEC subtypes, utilizing the GSVA R package. The heatmap indicated the significant activation of multiple oncogenic signaling pathways in TEC_1 and TEC_4 (Figure 3A). Notably, the TEC_0 and TEC_3 clusters exhibit a substantially higher proportion of malignant cells compared to the other three clusters (Figure 3B). Furthermore, we conducted a comprehensive comparison of GSVA scores for oncogenic signaling pathways among five TEC clusters, revealing distinctions between malignant and non-malignant cells. The findings suggest that the cell cycle, HIPPO, NOTCH, and RAS pathways are highly activated in non-malignant cells (Figure 3C–3G).
Figure 3. The characteristics of tumor-related pathways in TEC clusters. (A) Heatmap of 10 tumor-related pathway scores enriched in TEC cells. (B) Comparison of TEC clusters in malignant and non-malignant cells. (C–G) Comparison of GSVA score of each pathway between malignant and non-malignant cells in TEC_0 (C), TEC_1 (D), TEC_2 (E), TEC_3 cluster (F), and TEC_4 cluster (G). *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, Abbreviation: ns: not significant.
In order to evaluate the correlation between TEC clusters and prognosis, the ssGSEA scores of the marker genes (the top 5 DEGs of TEC clusters) were calculated for each TEC cluster using the TCGA cohort. The findings showed elevated TEC scores for the TEC_2 cluster in tumor tissues compared to normal tissues, whereas the TEC_3 and TEC_4 clusters displayed a contrasting pattern (Figure 4A–4E). Kaplan-Meier (KM) survival analysis showed that, in the high TEC score group, GC patients exhibited a worse prognosis across TEC_1, TEC_2, TEC_3, and TEC_4 clusters compared to the low-TEC score group. Conversely, TEC_0 showed no significant association with GC prognosis (Figure 4F–4J). Supplementary Figure 2 revealed that the difference in TEC scores between different clinical variables such as T stage, N stage, M stage, and pathological stage.
Figure 4. The associations between the five TEC cluster and prognosis of GC patients. (A–E) Comparison of five TEC scores in cancer and normal tissues, **P < 0.01, ***P < 0.001, ****P < 0.0001, Abbreviation: ns: not significant. (F–J) K-M curves of the high and low TEC score groups in the TEC_0 cluster (F), TEC_1 cluster (G), TEC_2 cluster (H), TEC_3 cluster (I), and TEC_4 (J).
Discovering hub genes related to TEC
To construct a risk signature, we conducted a differential analysis between tumor and normal tissues and identified a total of 2717 DEGs, comprising 2259 genes with up-regulation and 458 genes with down-regulation, as illustrated in Figure 5A. Regarding biological processes (BP), these genes exhibited enrichment in organelle fission, nuclear division, chromosome segregation, DNA replication, mitotic nuclear division, and mitotic cell cycle phase transition (Figure 5B). Furthermore, KEGG analysis indicated that these genes are predicted to positively regulate the cell cycle, Fanconi anemia pathway, DNA replication, motor proteins, and cellular senescence (Figure 5C).
Figure 5. Identification of the hub genes to construct a risk signature. (A) Volcano plot of differentially expressed genes of cancer and normal tissues in TCGA cohort. (B) GO analysis. (C) KEGG analysis. (D) Volcano plot of prognosis-related genes identified from univariate Cox regression analysis. (E) The trajectory of each independent variable with lambda. (F) Plots of the produced coefficient distributions for the logarithmic (lambda) series for parameter selection (lambda). (G) The multivariate Cox coefficients for each gene in the risk signature. (H) K-M curves of risk model constructed by 8 genes in TCGA cohort. (I) ROC curves of risk model constructed by 8 genes in TCGA cohort. (J) K-M curves of risk model constructed by 8 genes in GSE62254 cohort. (K) ROC curves of risk model constructed by 8 genes in GSE62254 cohort. (L) K-M curves of risk model constructed by 8 genes in GSE15459 cohort. (M) ROC curves of risk model constructed by 8 genes in GSE15459 cohort.
A total of 163 genes with predictive performance were identified through univariate Cox regression analysis (Figure 5D). The utilization of lasso-Cox regression analysis led to an additional decrease in the number of genes to 16, as illustrated in Figure 5E, 5F. Following a multivariate Cox regression analysis employing a stepwise regression approach, 8 genes were selected for incorporation into the risk signature (Figure 5G). These genes are ATP citrate lyase (ACLY), Sushi Domain Containing 1 (SUSD1), Angiopoietin 2 (ANGPT2), Ribonuclease A Family Member 1 (RNASE1), Neuronal Pentraxin 1 (NPTX1), RAB19, Kinesin Family Member 24 (KIF24), and Transcription Factor 20 (TCF20). The formula for calculating the risk score is as follows: RiskScore = 0.276 × ANGPT2 − 0.349 × KIF24 + 0.175 × RNASE1 + 0.459 × ACLY + 0.143 × NPTX1 − 0.229 × RAB19 + 0.369 × SUSD1 − 0.384 × TCF20. GC samples were categorized into high- and low-risk groups based on the median risk scores. Figure 5I, 5K, 5M show that the AUC values of the risk model varied between 0.7, 0.7, and 0.74 in the TCGA cohort, 0.67, 0.69, and 0.68 in the GSE62254 cohort, and 0.63, 0.61, and 0.59 in the GSE15459 cohort for 1-, 3-, and 5-year survival. In the TCGA and GEO cohorts, the KM survival analyses revealed that high-risk patients had worse clinical outcomes than low-risk patients, as shown in Figure 5H, 5J, 5L.
Identification of independent risk factors and nomogram development
To enhance the predictive accuracy of the risk model, we integrated clinical pathological features and risk scores and conducted univariate and multivariate Cox regression analyses. The results confirmed that the risk score was the most significant independent prognostic factor for GC (Figure 6A, 6B). Additionally, based on the risk score and clinical variables such as T stage and N stage, we constructed a nomogram to assess the clinical outcomes of GC patients (Figure 6C). The calibration curve demonstrated that the nomogram can effectively forecast the actual clinical outcomes (Figure 6D). As shown in Figure 6E, decision curve analysis (DCA) confirmed that the nomogram and risk score had better predictive efficacy in identifying high-risk patients compared to the N stage and T stage. Time-dependent ROC analysis demonstrated that, in the TCGA cohort, the area under the curve (AUC) of the risk score and nomogram was higher than other indicators (Figure 6F).
Figure 6. The development of a nomogram for predicting the prognosis of GC. (A, B) Univariate and multivariate Cox analysis of risk score and clinicopathological characteristics. (C) Nomogram model integrating the risk score and T stage, N stage was constructed. (D) Calibration curves for 1, 3 years of nomogram. (E) Decision curve for nomogram. (F) Comparison of predictive capacity of clinicopathological features and the nomogram using time-ROC analysis. ***P < 0.001.
Mutation and pathway analysis of the hub genes
Next, we analyzed the SNV mutation status of 8 genes in the risk model. The results showed that TCF20, ANGPT2, SUSD1, KIF24, ACLY, and NPTX1 exhibited SNV mutations in GC samples, while no SNV mutations were found in RAB19 and RNASE1 (Figure 7A). The probability of co-occurrence among the identified key genes and the top 10 most mutated genes was examined. Figure 7B shows no significant co-occurrence probability between mutations in RAB19 and NPTX1. However, ACLY, SUSD1, ANGPT2, and TCF20 exhibit a notable probability of co-occurrence with mutations in TTN, ARID1A, FAT4, and PCLO. Among the 8 genes, that only a minimal number of GC samples experienced gain/loss of copy number variation (CNV) (Figure 7C). To further elucidate the associations between the risk genes and GC, we analyzed the correlations between these genes and several molecular signatures of GC. The results demonstrated that RNASE1 had significantly negative correlations with aneuploidy score, fraction altered, and number of segments, whereas ANGPT2, KIF24, ACLY, and SUSD1 showed significantly positive correlations with homologous recombination defects, fraction altered, and number of segments (Figure 7D). Furthermore, the analysis of pathways indicated a strong correlation between eight genes and six pathways, including myogenesis, UV-response-DN, mitotic spindle, spermatogenesis, E2F targets, and G2M checkpoints (Figure 8A, 8B).
Figure 7. The characteristics of mutations of the genes included in the risk signature. (A) Waterfall diagram of SNV mutations of 8 key genes. (B) Colinearity and mutual exclusion analysis of 8 key genes and the 10 most mutated genes in tumors. (C) CNV mutations (gain, loss, none) of 8 key genes. (D) Correlation heatmap of 8 key genes with Aneuploidy Score, Homologous Recombination Defects, Fraction Altered, Number of Segments, and Nonsilent Mutation Rate.
Figure 8. Identification of the pathways in which the risk genes are implicated. (A) Heatmap showing the correlation between genes and pathways. (B) Heatmap displaying the enrichment scores for key pathways. *P < 0.05, **P < 0.01, ***P < 0.001.
Associations between risk genes and tumor immunity
Our data revealed significant positive correlations between RNASE1 and StromalScore, ImmuneScore, and ESTIMATEScore, whereas KIF24 exhibited significantly negative correlations with the three scores (Figure 9A, 9B). Based on the median of gene expression, we further evaluated the ImmuneScore difference between patients with high- and low-risk genes. The high-expression group for RNASE1 and NPTX1 genes exhibited a higher ImmunScore (Figure 9C). Correlation analysis reveals a significant negative correlation between ANGPT2, SUSD1, and ACLY with CD8 T cells, regulatory T cells, and activated NK cells, while ANGPT2, SUSD1, and ACLY show a significant positive correlation with resting NK cells, M0 macrophages, activated mast cells, and neutrophils (Figure 9D). Moreover, the MCPcounter algorithm revealed an association between ANGPT2 and neutrophils, endothelial cells, and fibroblasts (Figure 9E).
Figure 9. The relationship between the risk genes and immune landscape. (A, B) The correlation matrix of the risk genes and StromalScore, ImmuneScore, and ESTIMATEScore. (C) Comparison of high and low expression of 8 key genes and ImmuneScore. (D) Correlation between 8 key genes and immune cell score predicted by CIBERSORT analysis. (E) Correlation between 8 key genes and 10 immune cell types predicted by MCPcounter analysis. *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001.
Assessment of the predictive efficacy of risk models for immunotherapy
Immunotherapy, represented by ICIs, has become a crucial therapeutic approach for extending the survival of patients with advanced tumors. The response of patients to immunotherapy determines the ultimate treatment benefits. Therefore, we evaluated the predictive efficacy of risk features for immune response based on the IMvigor210 and GSE78220 cohorts. In the IMvigor210 cohort, patients in the low-risk group showed a favorable prognosis and a prolonged OS compared to the high-risk group (Figure 10A, p = 0.00011). Similar patterns were observed in the GSE78220 cohort (Figure 10F, p = 0.002). Furthermore, in the IMvigor210 and GSE78220 datasets, the proportion of stable disease (SD) and progressive disease (PD) was higher in the high-risk group compared to the low-risk group (Figure 10B, 10C, 10G, 10H). Of note, there is a significant difference in survival between patients in Stages I+II (Figure 10D, p = 0.0068) and Stages III+IV (Figure 10E, p = 0.0081) in the high- and low-risk groups.
Figure 10. The response of risk score to immune checkpoint inhibitors in IMvigor210 cohort and GSE78220 cohort. (A) Prognostic differences among risk score groups in the IMvigor210 cohort. (B) Differences in risk scores among immunotherapy responses in the IMvigor210 cohort. (C) Distribution of immunotherapy responses among risk score groups in the IMvigor210 cohort. (D) Prognostic differences between risk score groups in early-stage patients in the IMvigor210 cohort. (E) Prognostic differences between risk score groups in advanced patients in the IMvigor210 cohort. (F) Prognostic differences in risk score groups in the GSE78220 cohort. (G) Differences in risk scores among immunotherapy responses in the IMvigor210 cohort. (H) Distribution of immunotherapy responses among risk score groups in the GSE78220 cohort. *P < 0.05, ***P < 0.001.
Discussion
The dynamic interplay among various components within the TME, including cancer cells, stromal cells, immune cells, and the extracellular matrix (ECM), plays a crucial role in fostering tumor heterogeneity, clonal evolution, and multidrug resistance mechanisms [35–38]. This process ultimately contributes to tumor progression and metastasis. It has been reported that proangiogenic factors secreted by TECs and their interaction with tumor and immune cells are vital for tumor proliferation, angiogenesis, metastasis, and chemoresistance [9, 16, 22, 39]. The main objective of this research is to examine the variety of TEC clusters to analyze and categorize TECs in GC. We discovered five TEC clusters exhibiting unique characteristics. Previous studies have confirmed that a specific TEC cluster is essential for GC development, implying worse clinical outcomes [40]. Our study identified five TEC clusters, four of which were significantly associated with the prognosis of GC. Notably, differences in cell cycle, HIPPO, NOTCH, and RAS among TEC clusters may participate in GC growth and metastasis.
The Hippo pathway is a critical tumor suppressor signaling pathway that maintains tissue homeostasis by regulating cell differentiation and proliferation [41–43]. The activation of YAP/TAZ in the Hippo signaling pathway greatly enhances cell proliferation, migration, invasion, and anti-apoptotic processes in different solid tumors, including GC [43–45]. It has been shown that DUB1, a deubiquitinating enzyme, shows significantly elevated expression in GC tissues and is correlated with the activated TAZ protein and patient prognosis [46]. Mechanistically, DUB1 inhibits TAZ K48-linked polyubiquitination, reducing TAZ degradation and increasing its stability, thereby promoting GC stemness and progression. Furthermore, the Hippo pathway interacts with other important oncogenic signaling pathways, such as Wnt, Notch, EGFR, PI3K/AKT, MAPK, etc., [47]. This interaction affects the key components of the Hippo pathway, thereby determining cell fate.
Based on the predictive values of three TEC clusters, we developed a risk signature consisting of 8 genes. It consisted of three protective genes (RAB19, KIF24, and TCF20) and five risk genes (ACLY, SUSD1, ANGPT2, RNASE1, and NPTX1). In our study, SNV mutations were observed in RAB19 and RNASE1 without significant co-occurrence probability. SNV mutations affect protein activity or function, leading to GC development. Notably, CDC27 and FLG genes were mutated at the single-cell level but not detected in the corresponding tumor tissues, and they could promote cell growth by modulating the inflammatory response. KLF4 is a tumor suppressor, and KIF4 SNV may affect its DNA binding ability and apoptotic function [48]. Research has discovered a prevalent SNV in the zinc finger 2 region of the KLF4 gene in the foveolar-type gastric adenoma (FGA) tissues of Helicobacter pylori-negative patients. Compared to the wild-type KLF4 gene, the mutated KLF4 significantly inhibits cell proliferation and induces early apoptosis [48].
Additionally, we discovered that the 8 genes exhibited a significant correlation with 6 pathways. The protective genes showed a strong positive correlation with MITOTIC SPINDLE, E2F targets, G2M checkpoints, and spermatogenesis, while the risk genes (such as KIF24 and NPTX1) were significantly linked to myogenesis and UV-response-DN. E2F-1 participates in the G2/M checkpoint process, and as a transcription factor, it regulates the cell cycle, proliferation, and apoptosis. Overexpression of E2F-1 can decrease the expression of c-Myc, Skp2, Bcl-2, cyclin D1, and survivin, while increasing the expression of Bax, ultimately inhibiting the growth of GC cells [49, 50]. Furthermore, G2/M checkpoint arrest significantly inhibits the occurrence and development of GC [51, 52]. Hence, this data offers guidance for further investigating the regulatory mechanism of these risk genes in GC.
Recent evidence suggests that the interplay between TECs and the tumor immune microenvironment (TIME) plays a crucial role in tumor progression [9, 22]. Tumor-associated lymphatic endothelial cells (LECs) promote the growth and remodeling of lymphatic vessels by responding to factors in the TME, increasing lymph node metastasis [53]. Our research discovered a positive correlation between RNASE1 and ImmuneScore and a negative correlation between KIF24, ACLY, and ImmuneScore. These results imply potential crosstalk between risk genes and the TME in GC, highlighting their potential as therapeutic targets.
In the TME, various immune cell populations collectively determine the anti-tumor immune status of GC patients [35, 54, 55]. TECs can interact with these immune cells, establishing an immunosuppressive TME that facilitates tumor cell evasion from immune surveillance [8, 14]. In the risk model, memory B cells, plasma cells, regulatory T cells, activated NK cells, and activated dendritic cells show negative correlation with risk genes. Furthermore, our data suggest that a risk signature based on TECs can predict responses to immunotherapy. These findings provide new insights into the biological role of TECs in tumor immunity. However, the molecular mechanisms of TECs in GC and their potential value in immunotherapy still require further exploration.
Despite our progress in understanding the role of TECs in the development of GC, there remain some limitations. Firstly, the TEC clusters and the TEC-based risk signature we constructed are based on retrospective data from public databases, meaning our conclusions could be limited by the data collection and analysis methodologies. To ensure the accuracy and reliability of these findings, future research should be validated in a prospective, multicentric cohort of GC patients to better simulate the situation in clinical practice and provide stronger evidence to support our conclusions. Secondly, our study mainly focused on exploring the potential prognostic value of the TEC-based risk signature, that is, predicting the prognosis of GC patients by analyzing TEC characteristics. However, we have not fully revealed the specific molecular mechanisms involved by TECs in the occurrence and development of GC. To deeply understand how TECs affect the biological process of GC, future research needs to adopt a variety of experimental methods, including gene expression analysis, proteomics studies, and experiments using cell and animal models, to reveal the molecular action pathways of TECs in GC and potential therapeutic targets. Moreover, we also need to consider the heterogeneity of GC, that is, different types of GC patients may have different TEC characteristics and risk patterns. Therefore, future research should consider a more detailed classification of GC subtypes and explore the role of TECs in different subtypes, in order to provide more personalized treatment plans for patients.
Finally, although our research provides a new perspective for the prognosis assessment and treatment of GC, we must recognize that any single biomarker or risk signature cannot completely predict the progression of the disease. Future studies should explore integrating multiple biomarkers with clinical parameters to develop a more comprehensive prognostic model for risk assessment, enhancing the accuracy and effectiveness of GC treatment.
Conclusion
Our study delineated five distinct TEC clusters within GC tissues, four of which exhibit significant associations with the prognosis of GC patients. Utilizing lasso-Cox analysis, we identified 8 risk genes to formulate a TEC-based risk signature, which holds high clinical application value by accurately predicting the prognosis of GC patients and their response to immunotherapy.
Supplementary Materials
Supplementary Figures
Author Contributions
XFX designed, revised, and supervised the study. MF analyzed and organized the data. YH generated the figure and tables. MF and YH wrote this manuscript. All authors reviewed and approved the final manuscript.
Acknowledgments
The original data in our study were retrieved from the TCGA (https://portal.gdc.cancer.gov) and GEO (https://www.ncbi.nlm.nih.gov/geo) databases and we gratefully acknowledge the contributions from public available databases. In addition, we would also like to thank the Shengxin Bean Sprout Platform (http://www.sxdyc.com/index) for its help in data statistical analysis.
Conflicts of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Funding
No funding was used for this paper.
References
- 1. Thrift AP, Wenker TN, El-Serag HB. Global burden of gastric cancer: epidemiological trends, risk factors, screening and prevention. Nat Rev Clin Oncol. 2023; 20:338–49. https://doi.org/10.1038/s41571-023-00747-0 [PubMed]
- 2. Lordick F, Carneiro F, Cascinu S, Fleitas T, Haustermans K, Piessen G, Vogel A, Smyth EC, and ESMO Guidelines Committee. Gastric cancer: ESMO Clinical Practice Guideline for diagnosis, treatment and follow-up. Ann Oncol. 2022; 33:1005–20. https://doi.org/10.1016/j.annonc.2022.07.004 [PubMed]
- 3. Wang Z, Han W, Xue F, Zhao Y, Wu P, Chen Y, Yang C, Gu W, Jiang J. Nationwide gastric cancer prevention in China, 2021-2035: a decision analysis on effect, affordability and cost-effectiveness optimisation. Gut. 2022; 71:2391–400. https://doi.org/10.1136/gutjnl-2021-325948 [PubMed]
- 4. Malfertheiner P, Camargo MC, El-Omar E, Liou JM, Peek R, Schulz C, Smith SI, Suerbaum S. Helicobacter pylori infection. Nat Rev Dis Primers. 2023; 9:19. https://doi.org/10.1038/s41572-023-00431-8 [PubMed]
- 5. Jung YS, Xuan Tran MT, Park B, Moon CM. Association Between Family History of Gastric Cancer and the Risk of Gastric Cancer and Adenoma: A Nationwide Population-Based Study. Am J Gastroenterol. 2022; 117:1255–63. https://doi.org/10.14309/ajg.0000000000001837 [PubMed]
- 6. Hirata Y, Noorani A, Song S, Wang L, Ajani JA. Early stage gastric adenocarcinoma: clinical and molecular landscapes. Nat Rev Clin Oncol. 2023; 20:453–69. https://doi.org/10.1038/s41571-023-00767-w [PubMed]
- 7. Roy S, Kanda M, Nomura S, Zhu Z, Toiyama Y, Taketomi A, Goldenring J, Baba H, Kodera Y, Goel A. Diagnostic efficacy of circular RNAs as noninvasive, liquid biopsy biomarkers for early detection of gastric cancer. Mol Cancer. 2022; 21:42. https://doi.org/10.1186/s12943-022-01527-7 [PubMed]
- 8. Yang D, Kuang T, Zhou Y, Su Y, Shen J, Yu B, Zhao K, Ding Y. Tumor-associated endothelial cell prognostic risk model and tumor immune environment modulation in liver cancer based on single-cell and bulk RNA sequencing: Experimental verification. Int Immunopharmacol. 2023; 124:110870. https://doi.org/10.1016/j.intimp.2023.110870 [PubMed]
- 9. Maishi N, Annan DA, Kikuchi H, Hida Y, Hida K. Tumor Endothelial Heterogeneity in Cancer Progression. Cancers (Basel). 2019; 11:1511. https://doi.org/10.3390/cancers11101511 [PubMed]
- 10. Gu Y, Becker MA, Müller L, Reuss K, Umlauf F, Tang T, Menger MD, Laschke MW. MicroRNAs in Tumor Endothelial Cells: Regulation, Function and Therapeutic Applications. Cells. 2023; 12:1692. https://doi.org/10.3390/cells12131692 [PubMed]
- 11. Xie P, Tan SY, Li HF, Tang HD, Zhou JH. Transcriptomic landscape of endothelial cells: Key tumor microenvironment components indicating variable clinical outcomes in pancreatic ductal adenocarcinoma. Environ Toxicol. 2024; 39:572–82. https://doi.org/10.1002/tox.23881 [PubMed]
- 12. Yin Z, Wang L. Endothelial-to-mesenchymal transition in tumour progression and its potential roles in tumour therapy. Ann Med. 2023; 55:1058–69. https://doi.org/10.1080/07853890.2023.2180155 [PubMed]
- 13. Dzobo K, Senthebane DA, Dandara C. The Tumor Microenvironment in Tumorigenesis and Therapy Resistance Revisited. Cancers (Basel). 2023; 15:376. https://doi.org/10.3390/cancers15020376 [PubMed]
- 14. Nagl L, Horvath L, Pircher A, Wolf D. Tumor Endothelial Cells (TECs) as Potential Immune Directors of the Tumor Microenvironment - New Findings and Future Perspectives. Front Cell Dev Biol. 2020; 8:766. https://doi.org/10.3389/fcell.2020.00766 [PubMed]
- 15. Hida K, Maishi N, Annan DA, Hida Y. Contribution of Tumor Endothelial Cells in Cancer Progression. Int J Mol Sci. 2018; 19:1272. https://doi.org/10.3390/ijms19051272 [PubMed]
- 16. Maishi N, Ohba Y, Akiyama K, Ohga N, Hamada J, Nagao-Kitamoto H, Alam MT, Yamamoto K, Kawamoto T, Inoue N, Taketomi A, Shindoh M, Hida Y, Hida K. Tumour endothelial cells in high metastatic tumours promote metastasis via epigenetic dysregulation of biglycan. Sci Rep. 2016; 6:28039. https://doi.org/10.1038/srep28039 [PubMed]
- 17. Mustafa S, Koran S, AlOmair L. Insights Into the Role of Matrix Metalloproteinases in Cancer and its Various Therapeutic Aspects: A Review. Front Mol Biosci. 2022; 9:896099. https://doi.org/10.3389/fmolb.2022.896099 [PubMed]
- 18. Abdel-Hamid NM, Abass SA. Matrix metalloproteinase contribution in management of cancer proliferation, metastasis and drug targeting. Mol Biol Rep. 2021; 48:6525–38. https://doi.org/10.1007/s11033-021-06635-z [PubMed]
- 19. Gonzalez-Avila G, Sommer B, Mendoza-Posada DA, Ramos C, Garcia-Hernandez AA, Falfan-Valencia R. Corrigendum to "Matrix metalloproteinases participation in the metastatic process and their diagnostic and therapeutic applications in cancer" [Crit. Rev. Oncol. Hematol. 137, May (2019) 57-83]. Crit Rev Oncol Hematol. 2019; 138:172. https://doi.org/10.1016/j.critrevonc.2019.04.017 [PubMed]
- 20. Arima Y, Matsueda S, Saya H. Significance of Cancer-Associated Fibroblasts in the Interactions of Cancer Cells with the Tumor Microenvironment of Heterogeneous Tumor Tissue. Cancers (Basel). 2023; 15:2536. https://doi.org/10.3390/cancers15092536 [PubMed]
- 21. Raudenska M, Balvan J, Hanelova K, Bugajova M, Masarik M. Cancer-associated fibroblasts: Mediators of head and neck tumor microenvironment remodeling. Biochim Biophys Acta Rev Cancer. 2023; 1878:188940. https://doi.org/10.1016/j.bbcan.2023.188940 [PubMed]
- 22. Zeng Q, Mousa M, Nadukkandy AS, Franssens L, Alnaqbi H, Alshamsi FY, Safar HA, Carmeliet P. Understanding tumour endothelial cell heterogeneity and function from single-cell omics. Nat Rev Cancer. 2023; 23:544–64. https://doi.org/10.1038/s41568-023-00591-5 [PubMed]
- 23. Yin H, Guo R, Zhang H, Liu S, Gong Y, Yuan Y. A Dynamic Transcriptome Map of Different Tissue Microenvironment Cells Identified During Gastric Cancer Development Using Single-Cell RNA Sequencing. Front Immunol. 2021; 12:728169. https://doi.org/10.3389/fimmu.2021.728169 [PubMed]
- 24. Duan W, Xia S, Tang M, Lin M, Liu W, Wang Q. Targeting of endothelial cells in brain tumours. Clin Transl Med. 2023; 13:e1433. https://doi.org/10.1002/ctm2.1433 [PubMed]
- 25. Sathe A, Grimes SM, Lau BT, Chen J, Suarez C, Huang RJ, Poultsides G, Ji HP. Single-Cell Genomic Characterization Reveals the Cellular Reprogramming of the Gastric Tumor Microenvironment. Clin Cancer Res. 2020; 26:2640–53. https://doi.org/10.1158/1078-0432.CCR-19-3231 [PubMed]
- 26. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018; 36:411–20. https://doi.org/10.1038/nbt.4096 [PubMed]
- 27. Sanchez-Vega F, Mina M, Armenia J, Chatila WK, Luna A, La KC, Dimitriadoy S, Liu DL, Kantheti HS, Saghafinia S, Chakravarty D, Daian F, Gao Q, et al, and Cancer Genome Atlas Research Network. Oncogenic Signaling Pathways in The Cancer Genome Atlas. Cell. 2018; 173:321–37.e10. https://doi.org/10.1016/j.cell.2018.03.035 [PubMed]
- 28. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–7. https://doi.org/10.1089/omi.2011.0118 [PubMed]
- 29. Gao R, Bai S, Henderson YC, Lin Y, Schalck A, Yan Y, Kumar T, Hu M, Sei E, Davis A, Wang F, Shaitelman SF, Wang JR, et al. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat Biotechnol. 2021; 39:599–608. https://doi.org/10.1038/s41587-020-00795-2 [PubMed]
- 30. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]
- 31. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol. 2018; 1711:243–59. https://doi.org/10.1007/978-1-4939-7493-1_12 [PubMed]
- 32. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
- 33. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE III, Koeppen H, Astarita JL, Cubas R, Jhunjhunwala S, Banchereau R, Yang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018; 554:544–8. https://doi.org/10.1038/nature25501 [PubMed]
- 34. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, Berent-Maoz B, Pang J, Chmielowski B, Cherry G, Seja E, Lomeli S, Kong X, et al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell. 2016; 165:35–44. https://doi.org/10.1016/j.cell.2016.02.065 [PubMed]
- 35. Pitt JM, Marabelle A, Eggermont A, Soria JC, Kroemer G, Zitvogel L. Targeting the tumor microenvironment: removing obstruction to anticancer immune responses and immunotherapy. Ann Oncol. 2016; 27:1482–92. https://doi.org/10.1093/annonc/mdw168 [PubMed]
- 36. Furukawa N, Popel AS. Peptides that immunoactivate the tumor microenvironment. Biochim Biophys Acta Rev Cancer. 2021; 1875:188486. https://doi.org/10.1016/j.bbcan.2020.188486 [PubMed]
- 37. Ma C, Xi S, Sun H, Zhang M, Pei Y. Identifying the oncogenic roles of FAP in human cancers based on systematic analysis. Aging (Albany NY). 2023; 15:7056–83. https://doi.org/10.18632/aging.204892 [PubMed]
- 38. Ma C, Ma RJ, Hu K, Zheng QM, Wang YP, Zhang N, Sun ZG. The molecular mechanism of METTL3 promoting the malignant progression of lung cancer. Cancer Cell Int. 2022; 22:133. https://doi.org/10.1186/s12935-022-02539-5 [PubMed]
- 39. Zhou S, Ou H, Wu Y, Qi D, Pei X, Yu X, Hu X, Wu E. Targeting tumor endothelial cells with methyltransferase inhibitors: Mechanisms of action and the potential of combination therapy. Pharmacol Ther. 2023; 247:108434. https://doi.org/10.1016/j.pharmthera.2023.108434 [PubMed]
- 40. Li Y, Hu X, Lin R, Zhou G, Zhao L, Zhao D, Zhang Y, Li W, Zhang Y, Ma P, Ren H, Liao X, Niu P, et al. Single-cell landscape reveals active cell subtypes and their interaction in the tumor microenvironment of gastric cancer. Theranostics. 2022; 12:3818–33. https://doi.org/10.7150/thno.71833 [PubMed]
- 41. Qiao Y, Li T, Zheng S, Wang H. The Hippo pathway as a drug target in gastric cancer. Cancer Lett. 2018; 420:14–25. https://doi.org/10.1016/j.canlet.2018.01.062 [PubMed]
- 42. Messina B, Lo Sardo F, Scalera S, Memeo L, Colarossi C, Mare M, Blandino G, Ciliberto G, Maugeri-Saccà M, Bon G. Hippo pathway dysregulation in gastric cancer: from Helicobacter pylori infection to tumor promotion and progression. Cell Death Dis. 2023; 14:21. https://doi.org/10.1038/s41419-023-05568-8 [PubMed]
- 43. Kwon H, Kim J, Jho EH. Role of the Hippo pathway and mechanisms for controlling cellular localization of YAP/TAZ. FEBS J. 2022; 289:5798–818. https://doi.org/10.1111/febs.16091 [PubMed]
- 44. Kofler M, Kapus A. Nuclear Import and Export of YAP and TAZ. Cancers (Basel). 2023; 15:4956. https://doi.org/10.3390/cancers15204956 [PubMed]
- 45. Tiffon C, Giraud J, Molina-Castro SE, Peru S, Seeneevassen L, Sifré E, Staedel C, Bessède E, Dubus P, Mégraud F, Lehours P, Martin OCB, Varon C. TAZ Controls Helicobacter pylori-Induced Epithelial-Mesenchymal Transition and Cancer Stem Cell-Like Invasive and Tumorigenic Properties. Cells. 2020; 9:1462. https://doi.org/10.3390/cells9061462 [PubMed]
- 46. Wang D, Li Z, Li X, Yan C, Yang H, Zhuang T, Wang X, Zang Y, Liu Z, Wang T, Jiang R, Su P, Zhu J, Ding Y. DUB1 suppresses Hippo signaling by modulating TAZ protein expression in gastric cancer. J Exp Clin Cancer Res. 2022; 41:219. https://doi.org/10.1186/s13046-022-02410-5 [PubMed]
- 47. Kim W, Khan SK, Gvozdenovic-Jeremic J, Kim Y, Dahlman J, Kim H, Park O, Ishitani T, Jho EH, Gao B, Yang Y. Hippo signaling interactions with Wnt/β-catenin and Notch signaling repress liver tumorigenesis. J Clin Invest. 2017; 127:137–52. https://doi.org/10.1172/JCI88486 [PubMed]
- 48. Mishiro T, Shibagaki K, Fukuyama C, Kataoka M, Notsu T, Yamashita N, Oka A, Nagase M, Araki A, Kawashima K, Ishimura N, Maruyama R, Kinoshita Y, Ishihara S. KLF4 Mutation Shapes Pathologic Characteristics of Foveolar-Type Gastric Adenoma in Helicobacter pylori-Naive Patients. Am J Pathol. 2022; 192:1250–8. https://doi.org/10.1016/j.ajpath.2022.06.005 [PubMed]
- 49. Wang XT, Xie YB, Xiao Q. Lentivirus-mediated RNA interference targeting E2F-1 inhibits human gastric cancer MGC-803 cell growth in vivo. Exp Mol Med. 2011; 43:638–45. https://doi.org/10.3858/emm.2011.43.11.072 [PubMed]
- 50. Wei WY, Yan LH, Wang XT, Li L, Cao WL, Zhang XS, Zhan ZX, Yu H, Xie YB, Xiao Q. E2F-1 overexpression inhibits human gastric cancer MGC-803 cell growth in vivo. World J Gastroenterol. 2015; 21:491–501. https://doi.org/10.3748/wjg.v21.i2.491 [PubMed]
- 51. Liao X, Qian X, Zhang Z, Tao Y, Li Z, Zhang Q, Liang H, Li X, Xie Y, Zhuo R, Chen Y, Jiang Y, Cao H, et al. ARV-825 Demonstrates Antitumor Activity in Gastric Cancer via MYC-Targets and G2M-Checkpoint Signaling Pathways. Front Oncol. 2021; 11:753119. https://doi.org/10.3389/fonc.2021.753119 [PubMed]
- 52. Li Y, Zhang D, Yu K, Hu Y, Wu Q, Qian F, Wang Z. CMPD1 inhibited human gastric cancer cell proliferation by inducing apoptosis and G2/M cell cycle arrest. Biol Res. 2018; 51:11. https://doi.org/10.1186/s40659-018-0159-6 [PubMed]
- 53. Garnier L, Pick R, Montorfani J, Sun M, Brighouse D, Liaudet N, Kammertoens T, Blankenstein T, Page N, Bernier-Latamani J, Tran NL, Petrova TV, Merkler D, et al. IFN-γ-dependent tumor-antigen cross-presentation by lymphatic endothelial cells promotes their killing by T cells and inhibits metastasis. Sci Adv. 2022; 8:eabl5162. https://doi.org/10.1126/sciadv.abl5162 [PubMed]
- 54. Bagaev A, Kotlov N, Nomie K, Svekolkin V, Gafurov A, Isaeva O, Osokin N, Kozlov I, Frenkel F, Gancharova O, Almog N, Tsiper M, Ataullakhanov R, Fowler N. Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. Cancer Cell. 2021; 39:845–65.e7. https://doi.org/10.1016/j.ccell.2021.04.014 [PubMed]
- 55. Xiang X, Wang J, Lu D, Xu X. Targeting tumor-associated macrophages to synergize tumor immunotherapy. Signal Transduct Target Ther. 2021; 6:75. https://doi.org/10.1038/s41392-021-00484-9 [PubMed]