THBS2 is Closely Related to the Poor Prognosis and Immune Cell Infiltration of Gastric Cancer

Background: The potential functions of Thrombospondin 2 (THBS2) in the progression and immune infiltration of gastric cancer (GC) remain unclear. The purpose of this study was to clarify the role of THBS2 in GC prognosis and the relationship between THBS2 and GC immune cell infiltration. Material and Methods: The differential expression levels of THBS2 in the GC and cancer-adjacent tissues were identified using the TCGA databases and verified using real-time polymerase chain reaction (PCR), immunohistochemical staining and two datasets from Gene Expression Omnibus (GEO). THBS2 related differential expressed genes (DEGs) were identified and used for further functional enrichment analysis and Gene Set Enrichment Analysis (GSEA). Furthermore, a THBS2-related immune infiltration analysis was also performed. Kaplan-Meier and Cox regression analyses were utilized to illustrate the effects of THBS2 on the prognosis and clinical variables of GC. Finally, a nomogram was constructed to predict the survival probability of patients with GC. Results: The THBS2 expression in GC was significantly higher than that in cancer-adjacent tissues (p < 0.001), which was verified using real-time PCR, immunohistochemical staining and datasets from GEO. The 599 identified DEGs were primarily enriched in pathways related to tumorigenesis and tumor progression, including the focal adhesion pathway, signaling by vascular endothelial growth factor, and Wnt signaling. THBS2 expression was positively correlated with the enrichment of the macrophages (r = 0.590, p < 0.001), which was also confirmed by immunohistochemistry; however, negatively correlated with the enrichment of Th17 cells (r = 0.260, p < 0.001). The high expression of THBS2 was significantly correlated with the pathological grade (p < 0.01), histological grade (p < 0.05), histological type (p < 0.05), T stage (p < 0.001), and poor overall survival (OS) (P = 0.003) of GC. The constructed nomogram can well predict the 1-, 3-, and 5-years OS probability of patients with GC (C-index [95% confidence interval] = 0.725 [0.701–0.750]). Conclusion: THBS2 is closely related to the poor prognosis and immune infiltration of gastric cancer.


INTRODUCTION
Gastric cancer (GC) is the fifth most common malignancy and the third most common cause of death globally (Park et al., 2018). Currently, the treatment of GC includes surgery (Herrera-Almario and Strong, 2016), radiotherapy (Tey et al., 2017), neoadjuvant chemotherapy (Zhandossov et al., 2018), and immunotherapy (Niccolai et al., 2015). Due to the improvement in endoscopic screening technology, in countries with a high incidence of GC, such as Japan and South Korea, more than 50 percent of GC patients can be diagnosed at an early stage (So et al., 2021). However, in most countries with poor screening technology, most patients are diagnosed at an advanced GC stage, which leads to poor treatment outcomes (Choi et al., 2015;Inokuchi et al., 2018). Therefore, biomarkers related to the progression and GC prognosis should be explored to improve the therapeutic effect and reduce mortality owing to GC.
Thrombospondin-2 (THBS2) is a member of the matricellular calcium-binding glycoprotein family, which interacts with growth factors, cell receptors, and extracellular matrix (ECM), and THBS2 plays an important role in cell proliferation, adhesion and apoptosis (Bornstein et al., 2000;Zhuo et al., 2016). It has also been reported that THBS2 may be a serum biomarker in the diagnosis of colon cancer (Wang et al., 2016) and lung cancer (Weng et al., 2016). Bioinformatics analysis suggested that THBS2 might be a potential biomarker for GC (Cao et al., 2018). An in vitro study showed that THBS2 silencing inhibited the proliferation, migration and invasion of gastric cancer cells (Ao et al., 2018). The results of Chuanjun Zhuo showed that GC patients with low THBS2 expression had a better prognosis (Zhuo et al., 2016). However, there have also been contrasting reports. It has been reported that the expression of THBS2 is down-regulated in most GC patients, and the higher the expression of THBS2, the better the prognosis of GC patients, but the sample size was too small, only 14 cases (Sun et al., 2014). Therefore, the role of THBS2 in the prognosis of GC needs to be verified further.
Since immune checkpoint blocking therapy is used to treat several types of tumors, some patients have achieved significant clinical responses (Ribas and Wolchok, 2018). In September 2017, pembrolizumab (anti-programmed death 1 [PD1] antibody) was approved to treat GC or gastroesophageal junction cancer; however, the response rate was relatively low (Fashoyin-Aje et al., 2019). The effectiveness of tumor immunotherapy should be based on the premise that effector cells infiltrate into the tumor microenvironment (TME) (Kirkwood et al., 2012). Antictla-4 (T lymphocyte associated antigen 4) monoclonal antibody plays an antitumor role by prolonging T cell stimulation and restoring T cell proliferation (Engelhardt et al., 2006). T cells are the only type of tumor-infiltrating lymphocytes (TILs), and other common TILs cells include macrophages and NK cells. Therefore, understanding the tumor microenvironment and TILs of GC may lay a certain foundation for improving the effect of immunotherapy in GC. However, it has not been reported whether THBS2 can affect the TILs abundance of GC.
In this study, we analyzed the significance of THBS2 in GC in the TCGA database and Gene Expression Omnibus (GEO) database using bioinformatics analysis, including differential expressed genes (DEGs) analysis, functional enrichment analysis, gene set enrichment analysis (GSEA), immune cell infiltration analysis, clinical correlation analysis, and survival analysis. Moreover, a nomogram was also created to predict the overall survival (OS) of patients with GC.

Data Sources
We downloaded gene expression data and clinical information for GC, which included 32 cancer-adjacent tissues and 375 tumor tissues, from TCGA (https://portal.gdc.cancer.gov/, accessed time: July 31, 2021). Samples with incomplete clinical data were excluded. RNAseq data was converted by log2 in R software for subsequent analysis. The clinical characteristic of GC is presented in Table 1. Additionally, two microarray datasets (GSE54129 and GSE13911) containing GC and cancer-adjacent tissues were downloaded from GEO online database (https:// www.ncbi.nlm.nih.gov/geo/, accessed time: December 17, 2021).

THBS2 Differential Expression in Pan-Cancer and GC Tissues
Log2-converted TPM RNAseq data for cancer-adjacent and tumor tissues (TCGA) were obtained from the UCSC XENA. The differential expression between tumor and non-tumor tissues was tested using Wilcoxon Rank Sum Test and visualized through the box and scatter plots. Additionally, the receiver operating characteristic (ROC) curve was constructed to assess the predictive diagnostic value of THBS2 expression.

Real-Time PCR of THBS2 Expressions in GC and Adjacent Tissues
This study has been approved by the Ethics Committee of the First Affiliated Hospital of Guangxi Medical University and informed consent of all patients. From May to July 2021, cancer and para-cancerous biopsy specimens were incessantly collected from 16 patients diagnosed with GC at the Endoscopy Center of the First Affiliated Hospital of Guangxi Medical University (Nanning, Guangxi). All tissue specimens were immersed in RNA protective solution and quickly transferred to a -80°C refrigerator for preservation. All patients had not received any prior treatment for the tumor, including radiation or chemotherapy. Furthermore, patients who were complicated with other known tumors were excluded. The clinical information of these 16 patients is presented in Table 2.

Immunohistochemistry
We continuously collected tumor and para-tumor tissues from 80 GC patients without other known malignant tumors after surgery in Suqian First People's Hospital from January 2017 to April 2019. None had a history of chemotherapy or radiation therapy before surgery. All the tissues were modified into tissue chips for further IHC staining. Dewaxing, and hydration were performed before the THBS2 primary antibody (ab112543, Abcam, 1:1000) was added and incubated at 4°C for 12 h. After incubating with HRP labeled linked polymer (KIT-5009, MXB biotechnologies) at 26°C for 40 min, signal detection was performed using DAB (P0202, Beyotime). The expression of THBS2 was measured using ImageJ software and visualized with GraphPad Prism 8.0 after statistical analysis.

Identification of DEGs Between High and low Expression Groups of THBS2
According to the mean value of THBS2 expression, all the samples from TCGA datasets were divided into high and low expression groups. RNA-Seq expression data were processed using the DESeq2 package, and DEGs were defined with a p. adj <0.05 and an absolute logFC >1.5. Detailed gene expression is shown in the volcano map.

Functional Enrichment Analysis of DEGs
In this study, or. Hs.eg.db package was used for ID conversion, while the clusterProfiler package was used for enrichment analysis. The threshold conditions included: p. adj <0.05, and q value < 0.2. Statistical analysis and GSEA visualization are performed by the clusterProfiler package, using C2 collection from MSigDB. Permutations with 10,000 times were performed by gene set, and significance was set as an adjusted P < 0.05 and false discovery rate (FDR) < 0.25.

Immune Cell Infiltration
The relative tumor infiltration levels of immune cell types were quantified using ssGSEA of GSVA package (Hänzelmann et al., 2013) to interrogate gene expression levels in published signature gene lists (Bindea et al., 2013). To explore the correlation between THBS2 and the immune infiltration levels and the association of immune infiltration with the different expression groups of THBS2, Spearman's correlation and Wilcoxon signed-rank sum test was adopted. Furthermore, the average optical density (AOD) of THBS2 expression was calculated using ImageJ software based on immunohistochemistry results of 80 GC patients. According to the mean value of AOD of THBS2, all 80 GC patients were divided into high and low expression groups. Finally, 10 cases were randomly selected from the high expression and low expression groups for immunohistochemical staining of macrophages and NK cells. The immunohistochemical procedure was the same as part 2.5, and we used CD56 (MAB-0743, MXB biotechnologies) and CD68 (MAB-0041, MXB biotechnologies) to label NK cells and macrophages, respectively. The proportion of CD56/CD68 positive area was measured using ImageJ software and visualized with GraphPad Prism 8.0 after statistical analysis.

Clinical Correlation Analysis of THBS2 in GC
Wilcoxon signed-rank sum test and logistic regression were used to evaluate the correlation between THBS2 expression and clinicopathological variables. Further, univariate and multivariate Cox regressions were used to compare the effect of THBS2 expression and other clinicopathological variables on survival in patients with GC. Independent factors for GC prognosis were identified using multivariate Cox regression analysis. Furthermore, we collected clinicopathological data from 80 patients who underwent THBS2 immunohistochemistry to evaluate the relationship between THBS2 expression and clinicopathological variables. Chi-square tests were used to evaluate the relationship between gender, pathological type, residual tumor status, and THBS2 expression. Fisher's exact tests were used to evaluate the relationship between pathologic stage, T stage, N stage, primary treatment outcome, and THBS2 expression. Wilcoxon signed Rank-Sum test was used to evaluate the relationship between age and THBS2 expression.

Construction and Verification of Nomogram
A nomogram used to predict the probability of survival at 1, 3, and 5 years for patients with GC included all independent prognostic factors. The prognostic data were derived from published studies . Nomogram was formulated using R with the survival and rms package. We used the bootstrap method to calculate the parameters, and the number of repeated calculations for each group of samples was 200. The C-index is used to evaluate the prediction ability of the nomogram; the closer C-index is to 1, the stronger the prediction ability is.

THBS2 Differential Expression in Pan-Cancer and GC Tissues
Significant differential expression was observed in most of the 33 cancers, as illustrated in Figure 1A. Secondly, the THBS2 expression in 375 STAD samples and 32 para-cancer samples in the TCGA was compared. There was a significant difference in THBS2 expression between GC tissues and cancer-adjacent tissues (P < 0.001) ( Figure 1B). The THBS2 expression in GC tissues was also significantly higher than that in the corresponding adjacent tissues (P < 0.001) ( Figure 1C). The area under the curve (AUC) of the ROC curve was used to evaluate the predictive value of THBS2 for the diagnosis of GC. The result of the ROC curve ( Figure 1D) demonstrated that THBS2 had a high predictive value in distinguishing GC from cancer-adjacent tissues (AUC 0.864, 95%CI: 0.812-0.915).
The real-time PCR results further verified the reliability of our differential analysis at the transcriptional level ( Figure 2E Canceradjacent vs. Cancer, P 0.021). Moreover, the same results were also observed in IHC ( Figure 2F, Cancer-adjacent vs. Cancer, P 0.010).
We further analyzed THBS2 mRNA expression in GC and cancer-adjacent tissues RNA sequencing data from GSE13911 and GSE54129. The results showed that the THBS2 mRNA expression level was higher in GC than in cancer-adjacent tissues ( Figure 3A, B, p < 0.001). The result of the ROC curve demonstrated that THBS2 had a high predictive value in distinguishing GC from cancer-adjacent tissues ( Figure 3C, GSE13911, AUC 0.915, 95%CI: 0.843-0.987; Figure 3D, GSE54129, AUC 0.979, 95%CI: 0.957-1.000). Enriched biological processes (BP), cellular components (CC), and molecular function (MF) were used to comprehend the biological functions of DEGs better. The top three items enriched in BP, CC, MF, and KEGG of DEGs were visualized. As shown in Figure 4B, THBS2related genes were involved in KEGG (Protein digestion and absorption, Staphylococcus aureus infection, and ECMreceptor interaction); MF (extracellular matrix structural constituent, extracellular matrix structural constituent conferring tensile strength, and glycosaminoglycan binding); CC (collagen-containing extracellular matrix, cornified envelope, and collagen trimer); and BP (extracellular matrix organization, extracellular structure organization, and skin development). The details of GO and KEGG enrichment analysis results are illustrated in Table 3.

THBS2 Expression is Closely Related to the Immune cell Infiltration and Clinicopathological Variables of Gastric Cancer
An analysis based on online data showed that the higher the expression of THBS2, the more macrophages, NK cells, iDC, and eosinophils infiltrate, and the fewer Th17 cells, T helper cells, and NK CD56 bright cells infiltrate in GC tissues. (Figures 5A-G). We further verified the differential enrichment of macrophages and NK cells in the tumor and cancer-adjacent tissues by immunohistochemistry. The proportion of positive area of CD56 (NK cells) and CD68 (macrophages) in GC tissues was significantly higher than that in cancer-adjacent tissues (Figure 6; CD56, Cancer vs. Cancer-adjacent, P < 0.0001; CD68, Cancer vs. Canceradjacent, P < 0.0001). The proportion of the positive area of CD56 and CD68 in the THBS2 high expression group was significantly higher than that in the low expression group (Figure 7; CD56, High vs. Low, P < 0.01; CD68, High vs. Low, P < 0.0001). Finally, we verified that the proportion of CD56 and CD68 positive area was significantly correlated with THBS2 expression (Figure 7; CD56, r 0.636, P < 0.01; CD68, r 0.466, P < 0.05).
There were 375 patients with complete clinicopathologic data in the TCGA database who were included in the cohort. As illustrated in Figures 8A-D, the high expression of THBS2 was significantly correlated with the pathological grade (stage I vs. stage II and III & IV, P < 0.01), histological grade (G1 & G2 vs. G3, P < 0.05), histological type (Diffuse Type vs. Tubular Type, P < 0.05), and T stage (T1 vs. T2, T1 vs. T3, T1 vs. T4, P < 0.001) of patients. Moreover, there was no significant statistical correlation between the THBS2 expression level and N and M stages ( Figures 8E,F).
We further analyzed the relationship between THBS2 expression and clinicopathological data of 80 patients who underwent THBS2 immunohistochemistry. The results showed that THBS2 expression was significantly correlated with pathologic stage (P 0.044), T stage (P 0.003), histologic type (p < 0.001), and histological grade (P 0.027) in 80 GC patients who underwent THBS2 immunohistochemistry ( Table 4).

High Expression of THBS2 Was Associated With Poor Prognosis of GC
As can be seen from Figure 9A, the high THBS2 expression was significantly correlated with poor OS (P 0.003). Moreover, high THBS2 expression was also associated with worse PFI and DSS; however, the difference was not statistically significant ( Figures 9B,C)

Construction and Validation of Nomogram
All parameters persisting as independent predictors in the multivariate Cox models of the subgroups, including age, primary therapy outcome, THBS2 expression, pathological stage, and histologic grade, were integrated into the nomogram (Figure 10). The nomogram and calibration blots demonstrated that the nomogram-predicted 1-, 3-, and 5-years survival probabilities of GC were similar to the actual probabilities  0.725 [0.701-0.750]), indicating that the prediction was in good agreement with the actual survival probability of GC ( Figure 10B). These results indicated that the prediction ability of the nomogram for the survival probabilities might be clinically applicable also suggests that the nomogram was well-calibrated, with the mean predicted probabilities close to observed probabilities.

DISCUSSION
By analyzing the expression profile of GC patients in the open online database, we found that the expression level of THBS2 in GC tissues was significantly higher than that in canceradjacent tissues. There was also a significant difference in THBS2 expression in GC tissues and corresponding adjacent tissues, which was validated using RT-PCR and IHC. Additionally, the ROC curve demonstrated that THBS2 had a potential predictive value in distinguishing GC from canceradjacent tissues (AUC 0.864, 95% CI: 0.812-0.915). This observation suggests that THBS2 may serve as a potential indicator for GC diagnosis.
Totally, 599 DEGs based on THBS2 expression levels were screened, including 170 up-regulated and 429 down-regulated DEGs. The results of GO enrichment analysis showed that the  (Bonnans et al., 2014). Similarly, studies have also demonstrated that the ECM is an important component of all cancer markers that play a crucial role (Andreuzzi et al., 2020). ECM receptor interactions are critical in the tumor microenvironment (Chen et al., 2021), tumorigenesis, and tumor progression (Malik et al., 2015). The ECM protein regulates the metastasis of GC cells through the ITGB4/FAK/SOX2/HIF-1α signaling pathway induced by ECM receptor interaction (Gan et al., 2018). For KEGG pathway enrichment analysis, the DEGs are enriched in ECM-receptor interaction during S. aureus infection. The structural function of ECM is critical to maintaining normal cell activity (Giussani et al., 2019). ECM promotes tumorigenesis and progression by avoiding apoptosis, regulating cell growth, promoting tumor angiogenesis, and acquiring the ability of invasion and metastasis (Pickup et al., 2014;Malik et al., 2015;Poltavets et al., 2018;Eble and Niland, 2019). Signaling pathways associated with protein digestion and absorption are also critical for tumor progression (Mo et al., 2020). Therefore, treatment targeting ECM and tumor angiogenesis may be effective in preventing metastasis and recurrence of GC (Gan et al., 2018).
Additionally, GSEA analysis demonstrated that THBS2related enrichment pathways were as follows: focaladhesion, VEGF signaling, Wnt signaling, cytokinecytokine-receptor interaction, immunoregulatory lymphoid and a non-lymphoid cell, senescence, and autophagy in cancer. These pathways are closely related to cell adhesion, immune regulation, tumor autophagy, tumor progression, and tumor angiogenesis, respectively. Adhesion to ECM via specific focal adhesion points is an important step in cancer cell migration and invasion (Paluch et al., 2016;Ringer et al., 2017). VGFR related signaling pathway is important in the molecular pathogenesis of tumor growth and metastasis (Adamis and Shima, 2005). THBS2 may play a crucial role in the occurrence and progression of GC, so it is reasonable to speculate that THBS2 may be a potential therapeutic target for GC. Previous studies have confirmed that the interaction between tumor-infiltrating immune and malignant cells leads to the immune invasion of tumors, as the immune system plays a dual role by supporting both tumor progression and host defense (Bremnes et al., 2011). This study showed that the abundance of macrophages and NK cells in the tumor microenvironment of GC with high THBS2 expression was higher than that of GC with low THBS2 expression, while the abundance of Th17 cells was lower. Single tumorinitiating cells were detected to recruit polarized M2-like macrophages and assist evasion from immune clearance (Guo  (Doak et al., 2018). In some tumor types such as ovarian, prostate, and colorectal cancer, Th17 cells induce an antitumor immune response by recruiting cytotoxic effector T cells and producing effector cytokines, including interferon-gamma (Kryczek et al., 2009;Tosolini et al., 2011). Therefore, we speculated that the high expression of THBS2 might affect the disease progression and prognosis of patients with GC by inducing changes in the tumor microenvironment. THBS2, a member of the thrombospondin family of multidomain and secreted multicellular calcium-binding glycoproteins, mediates cell-to-cell and cell-to-matrix interactions (Ng et al., 2021). The expression level of THBS2 in colon cancer was significantly increased, and the higher the expression level of THBS2, the worse the OS of patients (Wang et al., 2016b). The high expression of THBS2 promotes the metastasis of colon cancer and is associated with an advanced clinical stage (Wang et al., 2016b). Our results suggest that the high expression of THBS2 was significantly correlated with the pathological grade (stage I vs. stage II and III & IV, P < 0.01), histological grade (G1 & G2 vs. G3, p < 0.05), histological type (Diffuse Type vs. Tubular Type, P < 0.05), T stage (T1 vs. T2, T1 vs. T3, T1 vs. T4, P < 0.001) of patients. The high THBS2 expression was significantly correlated with poor OS. Additionally, the subgroup survival analysis demonstrated that the prognosis of GC patients with THBS2 high expression was poor in T3 & T4, N1 & N2 & N3, M0, and stage I and II subgroups in OS, and stage III and IV subgroups in PFI. This study also indicates that the high expression of THBS2 may accelerate GC progression, leading to a poor prognosis. The higher the expression level of THBS2, the worse the degree of differentiation of GC.
Multivariate Cox regression analysis demonstrated that age, primary therapy outcome, THBS2 expression, pathological stage, and histologic grade were independent risk factors associated with poor OS in GC. A nomogram based on the above independent risk factors affecting OS of GC was also established to predict the survival probability of patients with GC. The survival probability of GC patients predicted by the constructed nomogram was similar to the actual survival probability, indicating that the created nomogram could be successfully used to predict the survival of GC patients. The accurate prediction of survival of cancer patients can provide essential information for individualized treatment of cancer patients (Erickson et al., 2014;Gratian et al., 2014). Our nomogram to predict the survival probability of GC has been verified using a calibration blot and should be clinically promoted.

CONCLUSION
Overall, the findings of the current research are summarized below: Firstly, the THBS2 expression level in GC was significantly higher than that in para-cancer tissues, and THBS2 may be a potential biomarker for GC diagnosis. Secondly, THBS2 may affect the progression and prognosis of GC by changing the tumor microenvironment and may be a potential therapeutic target for GC. Thirdly, histologic grade, primary therapy outcome, age, and THBS2 might be independent risk factors associated with poor OS in GC. Finally, the nomogram may provide more individualized prognostic information for patients with GC. However, neither the TCGA data nor the  hospital data collected were sufficiently large. Therefore, more information should be collected to verify the accuracy of the results.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of the First Affiliated Hospital of Guangxi Medical University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
GT contributed to the conception of the study; SZ performed the experiment and manuscript; HY contributed significantly to analysis and manuscript preparation; XX performed the data analyses and wrote the manuscript; LL and HH helped perform the analysis with constructive discussions.