This article is mentioned in:

Introduction

Globally, gastric cancer (GC) is one of the most common tumours of the alimentary system, ranking fifth in terms of incidence and third in terms of mortality among all cancer types in 2018 (1). Chemotherapy is one of the main treatments for GC, with oxaliplatin or cisplatin plus fluoropyrimidine being used as the first-line treatment (2). Nevertheless, the duration of the response to chemotherapy is limited, and the 5-year overall survival (OS) rate is still low in advanced GC. Therefore, there is an urgent requirement to find a novel and effective treatment for patients with GC.

The tumour microenvironment (TME) is composed of a variety of different components, including endothelial cells, fibroblasts, inflammatory mediators, mesenchymal cells, immune cells and stromal cells (3). The function of stromal cells and immune cells in the TME is crucial. However, our understanding of the TME is still far from sufficient. ESTIMATE, an algorithm designed by Yoshihara et al (4), can be used to evaluate stromal cells and immune cells according to the expression of specific genes. Recently, the ESTIMATE algorithm has been used in the study of various malignancies, such as acute myeloid leukaemia (AML) (5), clear cell renal cell carcinoma (ccRCC) (6) and glioma (7). However, to the best of our knowledge, the ESTIMATE algorithm has rarely used in the study of GC.

In the present study, ESTIMATE was utilized to perform TME evaluation to calculate stromal and immune scores. Hub genes associated with the TME were identified using Cytoscape software (version 3.7.1; http://cytoscape.org/). More importantly, the impact of these hub genes on survival was evaluated and validated, and the association between the expression level of RGS1 and the clinicopathological characteristics of patients with GC was analysed.

Materials and methods

The cancer genome atlas (TCGA) data collection and processing

The gene transcriptome and clinical profiles of 343 patients with stomach adenocarcinoma from TCGA database (http://tcga-data.nci.nih.gov/tcga/) were downloaded and processed (March 13, 2020). Immune, stromal and ESTIMATE scores were calculated by the ESTIMATE algorithm to predict the level of infiltrating immune cells, stromal cells and the tumour purity in tumour tissue, respectively (4).

Survival and association analysis

Kaplan-Meier survival curves were plotted, and the log-rank test was employed to compare the differences in OS for patients with GC based on the immune, stromal and ESTIMATE scores. Moreover, the association between score levels and multiple subgroups of clinical variables, including grade, American Joint Committee on Cancer (AJCC) stage and TNM stage (8), was also analysed.

Differentially expressed genes (DEGs), heatmaps and clustering analysis

The patients with GC were divided into high and low scoring groups according to the median immune, stromal and ESITIMATE scores (1088.7, 92.0 and 1106.1, respectively). DEGs were screened out with |log (FC)|>1 and false discovery rate (FDR)<0.05. Moreover, clustering analysis was applied to identify significant upregulated and downregulated gene sets between the immune and stromal subgroups. Corresponding heatmaps were plotted to illustrate the DEGs. Furthermore, the ‘VennDiagram’ package was utilized to identify and visualize the intersecting genes of upregulated and downregulated gene sets from the differential analysis of immune scores and stromal scores (9).

Functional enrichment analysis and protein-protein interaction (PPI) network construction

Intersect genes were subjected to comprehensive Gene Ontology (GO; http://geneontology.org/) analysis, including biological processes (BPs), cellular components (CCs) and molecular functions (MFs), and they were also subjected to Kyoto Encyclopaedia of Genes and Genomes (KEGG; http://www.kegg.jp/) pathway analysis. At the same time, the STRING database (https://string-db.org/) was used to construct the PPI network with the highest confidence (0.900).

Batch survival analysis of hub genes

The hub genes were identified from the PPI network using cytoHubba (10) from Cytoscape software (version 3.7.1; http://cytoscape.org/), and used for further survival analysis. Kaplan-Meier curves were plotted and the log-rank test was applied to evaluate the difference between the high and low groups based on the gene expression level with regard to the 5-year OS rate.

Collection of tissue samples and clinicopathological characteristics

Between August 2016 and December 2018, a total of 50 tumour tissue specimens were collected from patients with primary GC without metastasis who underwent curative gastrectomy in the Department of General Surgery of Tianjin Medical University General Hospital (Tianjin, China). All tissues were stored at −80°C. The patients were followed up, and corresponding clinicopathological characteristics, including age, sex, BMI, tumour location, differentiation degree, T stage, N stage and AJCC stage, were collected. The age range of all patients was 52–85 years, and 34 patients were male and 16 were females. All patients provided written informed consent before tissue collection.

Reverse transcription quantitative PCR (RT-qPCR) assay

Total RNA was extracted from tissue specimens with an RNAprep Pure Tissue kit (Tiangen Biotech, Co., Ltd). After quantification with a NanoDrop-2000 spectrophotometer (Thermo Fisher Scientific, Inc.), a total of 500 ng RNA was reverse transcribed into cDNA using FastQuant RT Super mix kit (Tiangen Biotech, Co., Ltd.) using the following temperature protocol: 15 min at 42°C and 3 min at 95°C. The cDNA was stored at −20°C. Subsequently, mRNA quantitation in tumour tissues was detected using a RT-qPCR assay with SYBR-Green qPCR Master Mix (Bimake Biotechnology) on a DNA Engine Opticon 2 (Bio-Rad Laboratories, Inc.). The thermocycling conditions were as follows: 95°C for 30 sec, then 40 cycles of 95°C for 15 sec, 60°C for 30 sec and 72°C for 30 sec, followed by dissociation at 95°C for 15 sec, 60°C for 1 min and 95°C for 15 sec. mRNA expression was quantified using the 2−∆∆Cq method (11) and normalized to the internal reference gene (β-actin). The primers are listed in Table I.

Table I. Summary of primer sequences used for reverse transcription-quantitative PCR.

Table I.

Summary of primer sequences used for reverse transcription-quantitative PCR.

Gene	Forward primer (5′-3′)	Reverse primer (5′-3′)
CXCR4	AATGACTTGTGGGTGGTTG	AAGAAAGCCAGGATGAGGA
PTGFR	GTTCCTGTCCACTCTGAAGG	TCACCATGCTCTGTAGTTTGA
RGS1	CAGCCAAGAAGATTAAAGCAC	TTGAGGAACCTGGGATAAGA
β-actin	CTCCTCCACCTTTGACGCTG	TCCTCTTGTGCTCTTGCTGG

Statistical analysis

R software (version 3.6.0; http://www.r-project.org/) was applied for the statistical analysis of TCGA data. SPSS (version 23.0; IBM, Corp.) was applied for the statistical analysis of the data from the Department of General Surgery of Tianjin Medical University General Hospital. The χ2 test and Fisher's exact test were applied for categorical variables, and Student's t-test was applied for continuous variables. The Kruskal-Wallis test and Wilcoxon rank sum test were used to compare ESTIMATE, immune and stromal scores between 2 or more subgroups of clinical variables, respectively. P<0.05 was considered to indicate a statistically significant difference.

Results

Stromal scores and immune scores are associated with the histological grade, AJCC stage, T stage and prognosis of GC patients

A total of 343 patients with GC from the TCGA database were included. The patients were divided into high and low scoring groups based on the median stromal, immune and ESTIMATE scores. There was no significant association between the ESTIMATE score and OS (P=0.4073; Fig. 1A). However, higher stromal scores were associated with worse outcomes in patients with GC (P=0.0178; Fig. 1B). There was no significant association between immune scores and OS (P=0.6624; Fig. 1C).

Figure 1. Stromal scores are associated with the survival outcome of GC, and immune scores and stromal scores are associated with GC tumour grade, AJCC stage and T staging. Kaplan-Meier survival curves suggested that (A) higher ESTIMATE scores (P=0.4073) were not significantly related to survival outcome, (B) stromal scores (P=0.0178) were significantly related to worse prognosis and (C) immune scores (P=0.6624) were not significantly related to survival outcome. (D-F) The distribution of ESTIMATE scores, stromal scores and immune scores for histological grade; (G-I) the distribution of ESTIMATE scores, stromal scores and immune scores for AJCC stage; and (J-L) the distribution of ESTIMATE scores, stromal scores and immune scores in T staging. *P<0.05, **P<0.01 and ***P<0.001. AJCC, American Joint Committee on Cancer; low, low score group; high, high score group; G, histological grade; stage, AJCC stage; T, AJCC T staging; GC, gastric cancer.

In addition, the associations between scores and clinical variables were also analysed. There were significant differences in ESTIMATE, stromal and immune scores (Fig. 1D-F) among different histological grades. Similar to those results observed for histological grade, there were significant differences in ESTIMATE, stromal and immune scores (Fig. 1G-I) for AJCC stage. For TNM staging, significant differences in ESTIMATE, stromal and immune scores (Fig. 1J-L) were found only between the subgroups of T staging. For ESTIMATE scores, no significant difference was observed in N and M staging (Fig. S1A and B). Similar results were also observed in stromal scores (Fig. S1C and D) and immune scores (Fig. S1E and F).

DEGs of stromal scores and immune scores for GC

The patients were divided into high and low groups based on the median scores. The results of the clustering analysis in the immune score group are shown in Fig. S2. The heatmap based on stromal scores is shown in Fig. S3. A total of 1,174 DEGs were identified based on immune scores, including 853 upregulated genes and 321 downregulated genes. Moreover, differential gene analysis was conducted in the same way based on the stromal score. In addition, 1,513 upregulated genes and 218 downregulated genes were identified (Fig. 2A and B). A Venn diagram was applied to distinguish 770 intersect genes, consisting of 644 upregulated and 126 downregulated intersect genes (Fig. 2A and B). The details of the DEGs of immune and stromal scores are provided in Tables SI and SII, respectively.

Figure 2. Identification of intersecting genes and functional enrichment analysis. Venn plots were generated to visualize the number of (A) upregulated and (B) downregulated intersecting genes. (C) The top 10 Gene Ontology terms in biological processes, cellular components and molecular functions were subjected to functional enrichment clustering analysis and visualized as a bubble chart. (D) The top 30 Kyoto Encyclopaedia of Genes and Genomes pathways were identified and visualized as a bubble chart.

Functional enrichment analysis and PPI network

The functional enrichment analysis results indicated that the intersecting genes were significantly associated with the immune response. The top 10 GO terms in each category of BP, CC and MF (Fig. 2C) were screened out. ‘Regulation of leukocyte activation’, ‘side of membrane’ and ‘carbohydrate binding’ were the top GO terms in the GO analysis. For the KEGG pathway enrichment analysis, the top 30 pathways are shown in Fig. 2D. Among them, ‘cytokine-cytokine receptor interaction’, ‘chemokine signalling pathway’ and ‘haematopoietic cell lineage’ were the top KEGG pathways, in which 770 intersecting genes might be involved. Intersecting genes were used to construct the PPI network, consisting of 294 nodes and 1,984 edges. As shown in Fig. 3A, a total of 71 TME-associated hub genes were identified by 12 algorithms from cytoHubba.

Figure 3. Hub genes of algorithms from cytoHubba and further batch survival analysis. (A) The PPI network data from STRING were further analysed by Cytoscape, and hub gene identification was performed by cytoHubba based on 12 algorithms. Three hub genes, (B) CXCR4 (P=0.0146), (C) PTGFR (P=0.0482) and (D) RGS1 (P=0.0150), were identified by batch survival analysis. The colours of the nodes represent the ranks of the nodes, and the darker the colour of the nodes, the higher the rank. Low, low expression group; high, high expression group; PPI, protein-protein interaction; DMNC, density of maximum neighbourhood component; MCC, maximal clique centrality; MNC, maximum neighbourhood component; EPC, edge percolated component.

Batch processing of hub gene survival analysis

To further elaborate the relationship between hub genes and outcome, the 71 hub genes from 12 algorithms were used for survival analysis. Among them, only three hub genes, CXCR4 (P=0.0146), PTGFR (P=0.0482) and RGS1 (P=0.0150), were identified with significant differences in 5-year OS rate between the high and low level groups (Fig. 3B-D).

Validation of survival of CXCR4, PTGFR and RGS1

To further verify the impact of the three hub genes on prognosis in the study population, 50 patients who were diagnosed with GC and underwent gastrectomy were included and followed up. These patients were classified into a high and low group based on the median 2−∆∆Cq value of each gene (Fig. S4A-C; CXCR4, 0.8353; PTGFR, 1.0189; RGS1, 1.000). However, the results revealed that the high and low CXCR4 and PTGFR groups exhibited no significant differences in OS rate in the population (Fig. 4A and B). However, the prognosis of the high RGS1 group (P=0.0172) was significantly worse than that of the low RGS1 group (Fig. 4C).

Figure 4. Survival analysis of three hub genes in the validation cohort. Kaplan-Meier curves revealed that higher expression of (A) CXCR4 and (B) PTGFR was not significantly associated with overall survival, but that higher expression of (C) RGS1 was significantly related to worse outcome.

Clinicopathological characteristics of the validation cohort

The clinicopathological characteristics of 50 patients were retrospectively collected. These patients were divided into low and high level groups according to the median 2−∆∆Cq value of RGS1. The mean age of all patients was 68.1 years old, and 34 patients (68%) were male. The baseline characteristics of the patients are shown in Table II. There was no significant difference in age between the low and high level groups (67.6±7.6 vs. 68.7±9.0 years). In addition, there were also no significant differences between the low and high level groups in terms of sex, BMI index, tumour location, T stage, N stage or AJCC stage. With regard to the degree of differentiation, however, the number of poorly differentiated tumours in the low level group was significantly less than that in the high level group (P<0.01; Table II).

Table II. Clinicopathological characteristics of validation cohort based on RGS1 expression.

Table II.

Clinicopathological characteristics of validation cohort based on RGS1 expression.

Characteristic	Low	High	P-value
Mean age ± SD, years	67.6±7.6	68.7±9.0	0.625
Male, n (%)	18 (72)	16 (64)	0.544
Mean BMI ± SD, kg/m2	22.71±3.42	22.16±2.23	0.509
Location, n			0.758
Gastric body	8	7
Antrum	17	18
Degree of differentiation, n			0.005
Poorly differentiated	4	14
Moderately differentiated	15	5
Highly differentiated	6	6
T stage, n			0.663
1	5	4
2	3	3
3	11	8
4	6	10
N stage, n			0.844
0	8	7
1	7	5
2	5	6
3	5	7
AJCC stage, n			0.693
I	6	5
II	9	7
III	10	13

[i] AJCC, American Joint Committee on Cancer.

Discussion

GC is a global health threat with high incidence and mortality. Recently, immunotherapy has been used to treat various solid tumours, playing an anti-tumour role by inhibiting immune checkpoints, such as programmed cell death protein 1 and programmed cell death ligand 1, to enhance the activity of T cells. For example, nivolumab, an emerging immune checkpoint inhibitor, has shown promising therapeutic prospects in phase III clinical trials of metastatic GC (12). Thus, the TME may play an indispensable role in the treatment of tumours. Due to the development of next-generation sequencing technology, we have a better understanding of the molecular biological level of tumours. To the best of our knowledge, TME-related genes have been identified in other malignancies, such as AML (5) and ccRCC (6), but in GC, they are limited. Therefore, the present study aimed to identify novel genetic targets for GC in the TME.

Immune, stromal and ESTIMATE scores were calculated in the present study based on the ESTIMATE algorithm. The results showed that higher stromal scores were significantly associated with worse outcomes in GC patients. In addition, higher scores were significantly related to higher tumour grades, advanced AJCC stage and higher AJCC T-stage with regard to immune, stromal and ESTIMATE scores. Thus, these findings suggested that stromal cells may play a more important role in the survival and prognosis of GC than immune cells. However, not only stromal cells but also immune cells may be involved in the development of GC. According to the results of GO analysis, the present study also provided some evidence for the biological basis of immunotherapy. In addition, functional enrichment analysis was performed and a PPI network was constructed. ‘Regulation of leukocyte activation’, ‘T-cell activation’ and ‘regulation of lymphocytes’ were the top 3 BP terms of GO analysis. The results revealed that the regulation of different immune cells might be involved in the progression of GC. Potential KEGG pathways, such as ‘cytokine-cytokine receptor interaction’, ‘chemokine signalling pathway’ and ‘haematopoietic cell lineage’, were identified. Moreover, a total of 71 hub genes of the PPI network were identified. Subsequently, higher expression of CXCR4 (P=0.0146), PTGFR (P=0.0482) and RGS1 (P=0.0150) was found to be significantly related to worse outcome.

CXCR4 is a member of the C-X-C chemokine receptor family. The binding of CXCL12, also known as stromal cell-derived factor-1, to CXCR4 triggers a series of downstream signalling pathways that result in several aspects of tumour progression, including angiogenesis, metastasis and survival (13,14). CXCR4 is highly expressed in >20 different types of cancer, such as breast (15), colorectal (16) and pancreatic ductal (17) cancer. A meta-analysis by Jiang et al (18) revealed that high CXCR4 expression is related to a worse prognosis in GC and that CXCR4 might be a prognostic biomarker in gastrointestinal cancer. In addition, Xiang et al (19) also found that CXCR4 is upregulated in GC tissues and is associated with more advanced tumour stage and poor outcomes. This study also reported that CXCR4 and CXCR2 can activate each other to promote the metastasis of GC. Together, these findings suggest that CXCR4 is closely related to the prognosis of GC. However, in the population from the Department of General Surgery of Tianjin Medical University General Hospital, there was no significant difference in survival rate between the low and high level groups. Although the effect of CXCR4 on patient survival was not verified, we hypothesize this may be due the small sample size that was included.

The PTGFR protein belongs to the G-protein coupled receptor family, which is the receptor of prostaglandin F2-α. Romanuik et al (20) found that PTGFR might have a positive effect on the proliferation of castration-recurrent prostate cancer based on transcriptome data. Akiyama et al (21) found that the expression of PTGFR is upregulated in human tumour endothelial cells and that PTGFR is expressed in human tumour blood vessels in vivo based on immunostaining. The results indicated that PTGFR may be applied for antiangiogenic therapy for RCC. However, in colorectal cancer, PTGFR is expressed at low levels in tumours but is hypermethylated in >40% of tumours (22). Accordingly, PTGFR may serve as a biomarker associated with GC, and the epigenetic modification of PTGFR may play a crucial role in the progression of GC. However, similarly to the results for CXCR4, there was no significant difference in OS of the validation population according to PTGFR.

RGS1, known as regulator of G-protein signalling 1, is a member of the regulator of the G-protein family that can activate G-protein signalling (23). RGS1 is highly expressed in a variety of immune cells, such as B lymphocytes (24,25), T lymphocytes (26), natural killer cells (27) and dendritic cells (28). The present results showed that high expression of RGS1 was related to a low degree of differentiation and worse outcomes. In metastatic head and neck squamous cell carcinoma (HNSCC), RGS1 expression is significantly upregulated, and knockdown of RGS1 was shown to inhibit the anchoring growth of HNSCC, which is not conducive to cell transformation (29). A previous study has shown that RGS1 expression is negatively correlated with the migration ability of regulatory T (Treg) cells (26). Moreover, the expression of RGS1 in Treg cells in peripheral blood from patients with metastatic castration-resistant prostate cancer was significantly higher than that of healthy donors. The overexpression of RGS1 may inhibit the migration of Treg cells and lead to Treg accumulation (30). The immune escape TME formed by the accumulation of Treg cells is conducive to the distant metastasis of prostate cancer. Therefore, we hypothesize that RGS1 may inhibit the migration of Treg cells in the TME, leading to their accumulation, which is conducive to the immune escape and distant metastasis of GC cells. The results of the present study have shown that RGS1 can also inhibit the differentiation of GC and increase the degree of malignancy, leading to the poor prognosis of GC patients.

In summary, the present study identified three TME-related hub genes from the PPI network according to the ESTIMATE algorithm. In addition, the effect of RGS1, one of the three hub genes, on survival was further verified in patients from the Department of General Surgery of Tianjin Medical University General Hospital, and RGS1 might play an important role in the differentiation of GC cells. However, there were also some limits to the present study. Firstly, the sample size of GC specimens from the Department of General Surgery of Tianjin Medical University General Hospital was relatively small, which may cause some bias to the results. Furthermore, only the gene expression level of RGS1 was verified, and the influence of RGS1 protein expression level on the prognosis of patients with GC was not explored. In the future further expansion of the sample size of GC specimens and revalidation of the relationship between the expression of CXCR4, PTGFR and RGS1 and the prognosis of patients with GC will be performed. Further mechanism studies will also be conducted to gain an in-depth understanding of the functions of RGS1 in the TME for GC.

In conclusion, the present study analysed transcriptome data from TCGA public database based on bioinformatics algorithms and identified genes associated with immune cells and stromal cells in the TME of GC. In addition, the validity of the three hub genes for predicting survival was verified. The present study lays a foundation for future research on the mechanism of the TME in GC.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

Funding information is not applicable.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author upon reasonable request.

Authors' contributions

ShiL and HY performed the data processing and data analysis, and these authors contributed equally to the manuscript. ShuL collected and interpreted the patient data regarding gastric cancer. DW and ZZ participated in the study design. WF led the design of the study and is the corresponding author. All authors have read and approved the manuscript, and confirm the authenticity of all the raw data.

Ethics approval and consent to participate

This study was approved by the Ethical Committee of Tianjin Medical University General Hospital (Tianjin, China) and written informed consent was obtained from all patients.

Patient consent for publication

Consent for publication was obtained from all patients.

Competing interests

The authors declare that they have no competing interests.

References