We selected the fragments per kilobase of transcript per million mapped read format of pan-cancer from The Cancer Genome Atlas (TCGA) database (https://www.cancer.gov/) (tumour, n = 10363; normal, n = 730) using the Sento Academic website (https://www.xiantao.love/) to analyse the differential expression of CHSY3 in pan-cancer (p-value < 0.05; Wilcoxon rank-sum test). Then, STAD expression matrix files (tumour, n = 375; normal, n = 32) were downloaded from TCGA database for paired and unpaired difference analyses by using the “limma” R package (p-value < 0.05; t-test). Table 1 exhibits the clinical information of high and low CHSY3 expression groups in STAD in TGCA database.

Samples were collected from patients who underwent gastric cancer surgery in the Department of Gastroenterology, The First Affiliated Hospital of Bengbu Medical College, from January 2017 to December 2018. A total of eight gastric cancer tissue samples and their paraneoplastic tissues were collected for Western blot experiments and 197 gastric cancer tissue samples and 30 normal paraneoplastic tissues for immunohistochemical staining. No patients had received chemotherapy, radiotherapy or biological treatment before surgery and were diagnosed with gastric cancer before and after surgery. The collected tissue samples were stored in a −80°C refrigerator immediately after surgery until protein extraction.

Rabbit anti-human antibody CHSY3 (100 μg) was obtained from OriGene China. The primary antibody to β-actin and CD3+/CD4+/CD8+ T cells alpha rabbit monoclonal antibody were provided by Cell Signaling Technology, Inc. (Danvers, MA, United States). Horseradish peroxidase (HRP)-conjugated anti-rabbit antibody was provided by Jackson ImmunoResearch Inc. (West Grove, PA, United States). Bovine serum albumin was obtained from Sigma-Aldrich (St. Louis, MO, United States). Skimmed milk and Tween-20 were purchased from Sangon Biotech Co., Ltd. (Shanghai, China).

All samples were fixed in 4% paraformaldehyde, embedded in paraffin, sectioned to 4 μm and adhered to slides. After de-affinity under different density gradients of xylene, the slides were rehydrated, and antigens were retrieved with citric acid buffer (pH 7.8, 0.1 M) at approximately 82°C for 24 min. The sections were evenly covered with endogenous peroxidase blocking solution for 15 min at room temperature to block the activity of endogenous peroxidase. After incubation with gene primary antibody or overnight at 4°C, the slides were washed gently with phosphate-buffered saline and incubated with biotin-conjugated secondary antibody for 10 min at room temperature and incubated with streptavidin peroxidase for 5 min. All sections were stained with hematoxylin and then cleaned. After sections were dried and cleared, immunohistochemical evaluation was performed.

Fresh gastric cancer tissue was obtained, and total protein was extracted. PowerPac HV high-voltage power supply (Bio-Rad Laboratories, Inc., CA, United States) was used for protein electrophoresis. Total protein was electrophoretically transferred to polyvinylidene fluoride membrane after sodium dodecyl sulphate polyacrylamide gel electrophoresis. After closure with fresh 5% skimmed milk, the membranes were incubated overnight at 4°C with CHSY3 primary antibody diluted at 1:800 and β-actin diluted at 1:3000. After washing with TBST, the membranes were incubated with horseradish peroxidase-coupled anti-rabbit immunoglobulin G diluted at 1:3000 for 2 h at room temperature. Finally, specific protein bands were detected for visualisation using the Bio-Rad Chemical XRS Imaging System (Bio-Rad Laboratories).

To quantify the differential expression of CHSY3 in STAD and normal tissues, we extracted parameters from the results of Western blot experiments by ImageJ (version: 1.8.0) software and visualized them by GraphPad Prism (version: 9) software for bar chart.

Correlation between CHSY3 expression (Affymetrix ID: 242100_at) and overall survival (OS; n = 881), first progression (FP) (n = 645) and post-progression survival (PPS) (n = 503) were analysed in the Kaplan–Meier Plotter website (https://kmplot.com/analysis/) by selecting the ‘Gastric Cancer’ module. Database samples were divided into high and low expression groups based on the mean value of CHSY3 expression (OS: cut-off value = 124, high = 234, low = 397; FP: cut-off value = 100, high = 258, low = 264; PPS: cut-off value = 125, high = 138, low = 246). The log-rank test was used to compare the survival differences between the two groups. We also used TCGA database to combine CHSY3 with clinical factors (gender, age, grade, T, N, M and stage) in the univariate and multivariate Cox analyses to look for independent indicators of STAD prognosis. Finally, we used the “rms” R package to construct a nomogram survival prediction system to integrate CHSY3 expression with clinical factors in each patient with STAD to predict 1-, 3- and 5-year survival rates. This system uses scoring criteria based on the magnitude of the regression coefficients of all independent variables and gives each independent variable a score for each value considered. A total score can be calculated for each patient. Subsequently, the probability of outcome occurrence for each patient was calculated using a conversion function between the scores and probability of outcome occurrence (Iasonos et al., 2008; Balachandran et al., 2015). In addition, the accuracy of the nomogram prediction system in predicting 1-, 3- and 5-year survival rates was assessed using calibration curves. If the model prediction curve matched, located above or found below the reference line, the predicted value was considered equal, greater or lower than the actual value, respectively.

CHSY3 expression-associated methylation sites were analysed using MEXPRESS (https://mexpress.be) by selecting the “STAD” section (p-value < 0.05). Then, SurvivalMeth (http://bio-bigdata.hrbmu. edu.cn/survivalmeth/) was used to analyse survival differences between the CHSY3 expression-related methylation site and STAD prognosis (log-rank test, p-value < 0.05).

To further explore the oncogenic correlation between CHSY3 and STAD, we used GSEA for oncogenic pathway enrichment. GSEA ranks genes according to the degree of differential expression in the two types of samples using a predefined set of genes, usually from functional annotations or results of previous experiments, and then tests whether the predefined set of genes is enriched at the top or bottom of this ranking table (Subramanian et al., 2005). Gene sets under pathways with |normalised enrichment score |>1, nominal p-value < 0.05 and false discovery rate q-value < 0.25 are generally considered significant.

The Hallmark consensus pathway gene set downloaded from the Molecular Signatures Database (http://www.gsea-msigdb.org/gsea/login.jsp) was subjected to gene set variation analysis (GSVA) using the “GSVA” R language package to assess the enrichment of the Hallmark pathway in the high and low CHSY3 expression groups (p-value < 0.05). GSVA is a non-parametric, unsupervised algorithm that assesses the enrichment of different metabolic pathways between samples by converting the sample-to-sample gene expression matrix into a sample-to-sample genomic expression matrix (Hänzelmann et al., 2013). Single-sample GSEA (ssGSEA) enrichment analysis was also performed using the “GSVA” R language package to assess the differences in immune cells and functions between the high and low CHSY3 expression groups (p-value < 0.05).

The estimation of stromal and immune cells in malignant tumour tissues using expression data (ESTIMATE) algorithm was used to assess the immune cell score, stromal score and tumour purity score of the high and low CHSY3 expression groups in the TIME of STAD. ESTIMATE analysis was performed using transcriptional profiles of cancer samples to evaluate the number of tumour cells, infiltrating immune cells and stromal cells, and the “estimateScore” function was used to calculate tumour purity, immune cell score and stromal cell score of all samples (p-value < 0.05).

The “limma” R language package was used to analyse CHSY3-related marker genes in high and low CHSY3 expression groups (p-value < 0.05). Sixty-six immune cell surface marker genes in STAD were analysed by the “reshape2” R package to determine correlation with CHSY3 expression (p-value < 0.05). The prognosis of CHSY3-related marker genes in STAD was further analysed using prognosis-related genes for the stepwise multivariate Cox proportional hazard regression analysis to obtain the optimal candidates and construct an immune-related risk model. The formula for calculating the risk score was as follows: R i s k s   c o r e = ∑ i = 1 n c o e f i × X i

The “coefi” and “Xi” represent the coefficient and expression level of each CHSY3 prognosis-related marker gene, respectively. According to the risk score of the model, TCGA samples can be divided into the high- and low-risk groups. The log-rank test was used to compare the survival differences between the two groups (p-value < 0.05). The “survivalROC” package was used to performreceiver operating characteristic (ROC) curve, and the Area Under Curve (AUC) values were obtained to evaluate the prognostic model’s reliability. To further analyse the prognostic risk of model scores, we performed univariate and multivariate Cox analyses by combining model scores with clinical factors (sex, age, grade, T, N, M, and stages) to assess whether model scores could be used as independent prognostic factors. Finally, the Tumor Immune Dysfunction and Exclusion (TIDE) (http://tide.dfci.harvard.edu/login/) was used to evaluate the risk of immune evasion and the effect of immunotherapy between the high- and low-risk groups of the model.

Statistical analyses were performed using the R software version 4.1.1 (R Foundation for Statistical Computing, Vienna, Austria). For quantitative data in the article data analysis, the significance of normally distributed variables was estimated using Student’s t-test, and non-normally distributed variables were analysed using the Wilcoxon rank-sum test. The log-rank test was used to compare data between two groups, and the Kruskal–Wallis test was performed to compare data between more than two groups. A p-value < 0.05 was considered significant.

Pan-cancer differential expression analysis revealed that CHSY3 was significantly differentially expressed in various cancers, with significant down-regulation in glioblastoma multiforme, lung adenocarcinoma, lung squamous carcinoma, etc., and significant up-regulation in STAD, breast-invasive carcinoma, kidney renal clear cell carcinoma, etc. (p-value < 0.05, Figure 1A). Paired difference analysis demonstrated significantly higher CHSY3 expression in STAD tissues than in paracancerous tissues (p-value < 0.05, Figure 1B) and exhibited consistent results in the unpaired difference analysis (p-value < 0.05, Figure 1C). At the tissue level, CHSY3 expression was higher in STAD tissues than in paraneoplastic tissues by immunohistochemical staining analysis (Figure 1D). This result was also observed at the protein level by Western blot experiments (Figure 1E). In addition, quantitative analysis of Western blot experiments demonstrated that CHSY3 expression was significantly higher in tumour tissues than in normal tissues in 7 out of 8 pairs of tissues (Figure 1F). Altogether, CHSY3 expression was significantly higher in STAD than in normal tissues.

The survival analysis revealed that the high CHSY3 expression group had poorer OS (Figure 2A), FP (Figure 2B) and PPS (Figure 2C) than the low CHSY3 expression group in the Kaplan–Meier Plotter database (p-value < 0.05). In addition, the univariate Cox analysis presented that T3, T4, N1, N3, M1, stage III, stage IV, age >65 years and CHSY3 expression were all risk factors significantly associated with prognosis (p-value < 0.05). The multivariate Cox analysis indicated that age and CHSY3 could be used as risk factors for STAD independent of the clinical factors in Table 2 (p-value < 0.05). This implies that CHSY3 can be used as an indicator of STAD prognosis.

To further evaluate the practical benefits of CHSY3 in clinical applications, we integrated CHSY3 expression data in STAD with clinical factors to construct a nomogram to predict patient survival at 1, 3 and 5 years (Figure 2D). The nomogram system survival prediction curves have a high coincidence with the calibration curve, which indicates the high accuracy of our nomogram (Figure 2E). Therefore, the proposed nomogram survival prediction system has a good clinical utility value.

MEXPRESS identified a total of 12 methylation sites associated with CHSY3 expression in STAD. Among these CHSY3 expression-associated methylation sites, cg10678749, cg11572844, cg02589568, and cg06610705 presented a significant positive correlation with CHSY3 expression (r > 0, p-value < 0.05), while cg04729562, cg18829263, cg02458929, cg26226142, cg09608073, cg24642372, cg20694933, and cg02571738 (r < 0, p-value < 0.05) had a significant negative correlation with CHSY3 expression (Figure 3A). SurvivalMeth revealed that seven CHSY3 expression-associated methylation sites had significant survival value. Cg06610705 and cg11572844 were positively correlated with CHSY3 expression, and high CHSY3 expression groups of these sites had a poorer prognosis than the low CHSY3 expression groups (Figure 3B, p-value < 0.05). Cg24642372, cg20694933, cg09608073, cg02571738 and cg26226142 were negatively correlated with CHSY3 expression, of which the low CHSY3 expression groups had a poorer prognosis than the high CHSY3 expression groups (Figure 3C, p-value < 0.05). Consistent with previous findings, the group with high CHSY3 expression-associated methylation sites had a poor prognosis of STAD.

GSEA showed that CHSY3 is involved in multiple signalling pathways in STAD. Table 3 exhibits the top 10 up-regulated signalling pathways associated with CHSY3 according to the GSEA score in STAD. We found that the transforming growth factor-beta (TGF-β), Wnt and Hedgehog signalling pathways were associated with the tumour mesenchymal pathway. The chemokine, Toll-like receptor and NOD-like receptor signalling pathways are involved in immune and inflammatory activation (Figure 4). Interestingly, high expression of CHSY3 in STAD did not perform a better prognosis.

To explore the effect of CHSY3 expression on the TIME in STAD, the Hallmark pathway gene set was enriched in the high and low CHSY3 expression groups using GSVA enrichment analysis. The results presented that the high CHSY3 expression group was up-regulated in immune activation and inflammatory signalling pathways, such as the TNFA signalling via NF-KB, allograft rejection, complement, IL6–JAK–STAT3 signalling, IL2–STAT5 signalling and inflammatory response. Moreover, the high CHSY3 expression group was up-regulated in epithelial–mesenchymal transition (EMT), TGF-β signalling, angiogenesis, Wnt/beta-catenin signalling, Notch signalling and other mesenchymal signalling pathways, while the low CHSY3 expression group demonstrated the opposite phenomenon (Figure 5A). The heatmap and differential analysis of immune cells and functions in STAD using the ssGSEA method demonstrated higher enrichment in the high CHSY3 expression group than in the low CHSY3 expression group (Figures 5A,B). The ESTIMATE assessment of the TIME in STAD also revealed that the immune cell score (Figures 5A,C) and tumour stroma score (Figures 5A,D) were significantly higher in the high CHSY3 expression group than in the low CHSY3 expression group, while the tumour purity score (Figures 5A,E) demonstrated the opposite phenomenon. These analyses consistently revealed that the high CHSY3 expression group was associated with severe immune cell infiltration. However, this immune advantage did not exhibit a survival advantage, which became the focus of our attention. Studies have demonstrated that the immune-exclusive tumour phenotype is characterised by numerous immune cells that are retained in the stroma surrounding the nest of tumour cells without penetrating these cells (Chen and Mellman, 2017). Immune-evading tumours are generally characterised by high TGF-β expression, myeloid inflammation and tumour neovascularisation as microenvironmental features (Fukumura et al., 2018; Metelli et al., 2018; Hegde and Chen, 2020). Figure 5F displays that signalling pathways such as EMT, angiogenesis, TGF-β and Wnt were significantly higher in the high CHSY3 expression group than in the low CHSY3 expression group. Therefore, a high CHSY3 expression mediates TIME immune evasion in STAD.

Currently, programmed death-ligand 1 (PD-L1) is a remarkable discovery in immunotherapy, and its expression level is an important predictor of the response to anti-PD-1/L1 therapy (Fukumura et al., 2018). We further selected PD-L1 to assess the effect of immunotherapy, and the results revealed that PD-L1 expression was higher in the high CHSY3 expression group than in the low CHSY3 expression group, implying that the high CHSY3 expression group may have immune-evasion ability and better immunotherapy effects (Figure 5G). Studies have also reported that the newly discovered TIDE score is an effective predictor of anti-PD1 and anti-CTLA4 therapy among all efficacy characteristics of ICI suppression therapy (Jiang et al., 2018). The predictive function of the TIDE score for efficacy was stable regardless of the degree of tumour-infiltrating cytotoxic T-cells (Jiang et al., 2018). The TIDE score analysis indicated that the high CHSY3 expression group had significantly higher score than the low CHSY3 expression group; this finding also demonstrated that high CHSY3 expression mediates the immune-evasion status of the TIME in STAD (Figure 5H).

Using the PubMed database, we retrieved 66 marker genes localised on the surface of immune cells. Proteins encoded by these genes, also known as immunomodulators, are classified as immunostimulators and immunoinhibitors, and studies have demonstrated that immunomodulators have significant effects on prognosis (Hadden, 1993; Hengge et al., 2001). In the correlation analysis, 46 of 66 immunomodulators in STAD were associated with CHSY3 expression (Table 4). In the differential expression analysis, 41 of 46 CHSY3-related immunomodulators were significantly different between the high and low CHSY3 expression groups (Figure 6A, p-value < 0.05). In the Cox survival analyses, CHSY3 expression-related genes, namely, CSF1R, TGFB1, CXCR4, and TNFSF18, were prognosis-related risk factors, while CTLA4 was a prognosis-related favourable factor (Figure 6B, p-value < 0.05). Moreover, we constructed a CSHY3-related immune risk model based on prognostic immunomodulators. The model risk scoring formula is as follows: risk core = (coefficient × CXCR4 expression) + (coefficient × CTLA4 expression) (Figure 6C). The immune risk model was constructed to classify TCGA samples according to risk and prognosis (Figure 6D). According to the model, the survival of the high-risk group was significantly worse than that of the low-risk group (Figure 6E). The accuracy assessment of the ROC curve for the risk model revealed that the risk scoring of the model (AUC = 0.706) and the model risk scoring combined with clinical factors (AUC = 0.738) had a high accuracy. In addition, the univariate Cox analysis bared that age, stage, T, N and risk score were prognostically relevant risk factors (Figure 6G), and the multivariate Cox analysis revealed that age and risk score could be used as independent prognostic risk factors (Figure 6H). These analyses suggest that the CHSY3-mediated TIME of STAD is associated with a poor prognosis. The TIDE score indicated that the high-risk group demonstrated a significantly higher risk of immune evasion than the low-risk group (Figure 6I). This is consistent with the previous conclusion that CHSY3 mediates the immune-evasion status in the TIME of STAD.

To further verify that CHSY3 mediates TIME immune evasion in STAD, we selected tissue samples expressing “+” and “+++” of CHSY3 in STAD in Figure 1D for immunohistochemical staining analysis. The purpose of staining was to observe the infiltration of CD3⁺ T cells, CD4⁺ T cells and CD8⁺ T cells in STAD tumour tissues and surrounding tissues to determine whether the high expression of CHSY3 in TIME is consistent with the characteristics of immune evasion. The results showed that CD3⁺ T cells, CD4⁺ T cells and CD8⁺ T cells in CHSY3 expressing “+++” tissues were mainly clustered in the tumour peripheral stroma, with few immune cells penetrating the stroma into the tumor parenchyma (Figure 7A). In contrast, in CHSY3 expressing “+”, tumour peripheral CD3⁺ T cells, CD4⁺ T cells and CD8⁺ T cells clustered less and more immune cells penetrated the stroma into the tumour parenchyma (Figure 7B). In summary, this phenomenon further confirms that the CHSY3 high expression group mediates immune evasion in the TIME of STAD.