High expression of suppressor of cytokine signaling 1 (SOCS1) associates with favourable responsiveness to neoadjuvant chemoradiotherapy in esophageal squamous cell carcinoma

Background: The responsiveness to preoperative neoadjuvant therapy in esophageal squamous cell carcinoma (ESCC) is signicantly related to the surgical effect and long-term prognosis of patients. Biomarkers for predicting the effect of preoperative neoadjuvant therapy of ESCC are urgently needed in clinical practice. M2 macrophage in the tumor microenvironment (TME) has been conrmed to have a denite correlation with neoadjuvant therapy. However, the research on the potential functional genes of M2 macrophage has not been carried out. The purpose of this study is to systematically screen the potential functional genes of M2 macrophages and verify the correlation between the expression of screened genes and responsiveness to neoadjuvant therapy in ESCC. Methods: The Cancer Genome Atlas (TCGA) database was used to screen the potential functional genes of M2 macrophages systematically. The correlation between the screened genes and responsiveness to neoadjuvant therapy in ESCC was analyzed in the training set GSE45670 and the validation set GSE104958. The correlation of the screened genes was conrmed in a cohort of 27 patients using immunohistochemistry (IHC). Results: A total of 11 M2 macrophage potential functional genes were screened out. The suppressor of cytokine signaling 1 (SOCS1) and ORAI calcium release-activated calcium modulator 3 (ORAI3) were selected for subsequent verication. The expression of SOCS1 and ORAI3 in 27 cases of ESCC was evaluated by the H-score system. The result showed that the high expression of SOCS1 was signicantly correlated with favourable responsiveness to preoperative neoadjuvant therapy (AUC = 0.830, 95% CI: 0.647-1, P = 0.006). No signicant correlation was found between the expression of ORAI3 and neoadjuvant therapy in ESCC (AUC=0.688, 95% CI: 0.478-0.899, P=0.117). Conclusion: The high expression of SOCS1 in the local tumor microenvironment is signicantly correlated with favourable responsiveness to neoadjuvant therapy in patients with ESCC, which can be used as a biomarker to predict the responsiveness to neoadjuvant therapy, and has certain clinical value.


Introduction
Esophageal cancer is the sixth most common tumor in the world [1]. The symptoms of patients with early esophageal cancer are not typical, and the disease is often in the middle and late stages when swallowing obstruction occurs. For patients with locally advanced esophageal cancer, preoperative neoadjuvant therapy is often used to reduce the tumor stage, and the following radical resection can signi cantly improve the R0 resection rate and improve the prognosis [2]. However, the responsiveness of neoadjuvant therapy varies from patient to patient, according to the disease stage, age, immune status and tumor heterogeneity [3]. Patients with poor responsiveness to neoadjuvant therapy suffer from not only di culties in receiving effective treatment, but also unnecessary economic pressures. Therefore, it is crucial to predict the responsiveness to neoadjuvant therapy. Many studies have shown that tumor microenvironment (TME) is closely related to tumor treatment e cacy and prognosis [4,5]. TME, composed of immune cells and their subtypes, stromal cells, various small molecules and other components, is closely related to the occurrence, development and metastasis of tumor, and is an important part of tumor heterogeneity [6]. M2 macrophage is an important immune cell in TME, which can promote tumor chemoradiotherapy resistance. Previous studies have found that the higher the in ltration abundance of M2 macrophages in tumor tissue, the worse the responsiveness to preoperative neoadjuvant therapy, and the more di cult it is to achieve pathological remission [7]. Moreover, this study also found that the genes related to M2 macrophage function are also related to the responsiveness to preoperative neoadjuvant therapy for ESCC [7]. Therefore, M2 macrophage function-related genes were systematically screened using TCGA and Gene Expression Omnibus (GEO) database, and their correlation with the responsiveness to the preoperative neoadjuvant therapy of ESCC were analyzed in the training set and validation set. Veri cation was carried out using the preoperative specimens of ESCC in our hospital.

Materials And Methods
Datasets and preprocessing TCGA Esophageal Cancer (ESCA) RNA expression matrix in the format of log2(x + 1) transformed RSEM normalized count and phenotype le, as well as survival data, were downloaded from https://xenabrowser.net/datapages/. The series matrix les of two GEO datasets, GSE45670 and GSE104958, which contain RNA expression pro les and pathological responsiveness information, were obtained from the NCBI website (https://www.ncbi.nlm.nih.gov/geo/). The RNA expression pro les of GSE45670 and GSE104958 were log2 scale. The only probe set was selected to represent the expression level of corresponding genes based on the maximum interquartile range (IQR). All probe sets were annotated into ENTREZID. There were 11854 genes retained common to all three cohorts.
The GSE45670 dataset included the biopsies from 28 patients with ESCCs who received neoadjuvant chemoradiotherapy (nCRT) and surgery, in which 11 cases reached pathological complete remission (pCR) and the other 17 cases did not. All 41 tissue samples in the GSE104958 dataset were biopsies during the endoscopic examination before administering the rst course of neoadjuvant chemotherapy with docetaxel, cisplatin and 5-uorouracil, in which 10 cases reached pCR.

Development of pathological response predicting model
To identify genes that in uence the function of M2 on responsiveness to nCRT, an approach that comprised the two main steps was adopted. Firstly, candidate genes that potentially in uence the function of M2 were identi ed with the method similar to Peng's approach [9]. A total of 79 patients had squamous cell carcinoma with R0 resection and had a complete record of disease-speci c survival (DSS) intervals and censored events in the TCGA ESCA data set. The baseline information for this cohort of TCGA is summarized in Table 1. Cox regression analysis had shown that the abundance of M2 macrophage was strongly associated with DSS in these 79 patients (HR = 4.562, 95% CI: 1.539-13.52, P = 0.00618, Supplementary Table 1). A linear model hazard = a×M2 + b×V + d×M2×V + c was solved using Cox regression for DSS. In the above formula, a, b and d represent the regression coe cients of the abundance of M2 macrophage, gene expression level, and the interaction of M2 macrophage and gene expression level respectively. The regression analysis of 11854 genes showed that the three coe cients (a, b and c) of 279 genes were signi cant (crude P < 0.05, Supplementary Table 2). Hence, 279 genes were selected as potential functional genes of M2 macrophages. In the GSE45670 dataset, the ratio of 279 potential functional genes to the abundance of M2 macrophage genes was used as the candidate variable, and the Least Absolute Shrinkage and Selection Operator (LASSO) regression was used to construct the classi er predictive of pCR. Among the 279 potential functional genes, 11 genes were nally selected with the minimum λ value. A weighted linear model was used as the predictive score for pathological response. However, to further verify the genes in the clinical setting, only the rst two genes with the largest absolute regression coe cient were selected for subsequent veri cation, including suppressor of cytokine signaling 1 (SOCS1) and ORAI calcium release-activated calcium modulator 3 (ORAI 3). All the above analysis was done with R language (Version 3.6.2). The ow chart of the analysis was shown in Fig Evaluation of responsiveness to nCRT Iodine water angiography and chest enhanced CT were used to record the changes of lesion size before and after neoadjuvant therapy. RECIST solid tumor e cacy evaluation standard 1.1 was used to evaluate the responsiveness to neoadjuvant chemoradiotherapy [10]: complete remission (CR): the target lesion disappeared; partial remission (PR): the sum of the longest single diameter of the lesion decreased by more than 30%; progression (PD): the sum of the longest single diameter of the lesion increased by more than 20% or new lesions appeared; stability (SD): the sum of the longest single diameter of the lesion decreased but did not reach PR or increased but did not reach PD.

Immunohistochemistry
Para n-embedded specimens of ESCC were sectioned with a thickness of 3 µm. The slices were baked at 60 ℃ for 2h and dewaxed in xylene for 15 minutes. 3% hydrogen peroxide was used to remove endogenous peroxidase activity. Statistical analysis SPSS 19.0 statistical software was used for statistical analysis. Receiver operating characteristic (ROC) curve analysis was used to evaluate the predictive e cacy of Score1 and Score2 in response to nCRT. In univariate analysis, the chi-square test was used to analyze the correlation between remission and gender, stage, chemotherapy regimen and other classi cation variables, and the Kruskal Wallis test was used to analyze the correlation between remission and age, Score1, Score2. P < 0.05 was considered statistically signi cant.

Results
Screening and identi cation of potential function genes of M2 macrophages The linear risk model was constructed in TCGA data, and Cox regression was used for DSS. After regression analysis of 11854 genes, 279 genes had signi cant effects on M2 macrophage in ltration abundance, gene expression level, and the three regression coe cients of their interaction (P < 0.05). These 279 genes can be used as potential functional genes of M2 macrophages (Supplementary  Table 2).

Establishment and validation of the prediction model for nCRT
The result of LASSO regression constructing the classi er predictive of pCR was presented in Supplementary Table 3. In the GSE45670 dataset, 11 genes were screened out from 279 potential functional genes to predict the response of neoadjuvant therapy of ESCC (Table 2). To meet the needs of actual detection, only SOCS1 and ORAI3 genes with larger regression coe cients were selected as followup validation genes. The sum of their regression coe cients multiplying by the ratio of the two genes expression to the abundance of M2 macrophage was calculated as the nal index, namely Score for predicting pCR. Using the GSE45670 dataset as the training set, the area under the ROC curve (AUC) predicted by Score was 0.904 (95% CI: 0.792-1, P < 0.001, Fig. 2a). The GSE104958 dataset was used as the validation set, and the Score was calculated by the same method. The AUC was 0.706 (95% CI: 0.5033-0.9096) (Fig. 2b). The Score of pCR group was higher than that of the non-pCR group (H = 3.775, P = 0.052, Fig. 3a). There was no signi cant difference in the in ltration abundance of M2 macrophages (H = 2.304, P = 0.129) (Fig. 3b). In this small sample size cohort, using the multiple gene expression ratios is more effective than using M2 abundance alone. To explore the correlation between the in ltration abundance of M2 macrophage and the expression of SOCS1 and ORAI3, MCPcounter score that represented the abundance of 10 immune cell types got from TCGA ESCA data set was used for regression analysis (Fig. 4a). The negative correlation between the expression ratio of SOCS1 to M2 macrophage, namely Score1 and the in ltration abundance of M2 macrophage (Fig. 4b). A similar result of ORAI3 was represented in Fig. 4c. SOCS1 overexpression is associated with favourable responsiveness to preoperative neoadjuvant therapy To con rm the correlation between the expression of SOCS1 and ORAI3 in tumor tissues and preoperative neoadjuvant responsiveness, 27 cases of ESCC para n sections were stained by immunohistochemistry (Fig. 5). ROC curve was used to analyze the predictive effect of CD163 H-score, Score1 and Score2 on the e cacy of neoadjuvant therapy. The results showed that the AUC of CD163 H-score for predicting curative effect was 0.556(95% CI: 0.307-0.804, P = 0.643). Score1 that re ected the expression of SOCS1, was signi cantly correlated with favorable responsiveness to preoperative neoadjuvant therapy (AUC = 0.830, 95% CI: 0.647-1, P = 0.006, Fig. 6a). There was no signi cant correlation between Score2 that re ected the expression of ORAI3 and neoadjuvant therapy in ESCC (AUC = 0.688, 95% CI: 0.478-0.899, P = 0.117). Univariate analysis showed that Score1 in patients with PR was signi cantly higher than that in patients with PD + SD (P = 0.006) (Fig. 6b). There was no signi cant difference in CD163 Hscore between patients with PR and PD + SD (H = 0.216, P = 0.642). Remission was signi cantly correlated with gender, stage (P < 0.05) ( Table 3).

Discussion
Previous studies have con rmed a signi cant correlation between the in ltration degree of M2 macrophages in TME and the responsiveness to neoadjuvant therapy [7,11]. However, the mechanism of M2 macrophages participating in neoadjuvant therapy is still unclear. In this study, all the functional genes of M2 macrophages were systematically screened by the TCGA database for the rst time. The functional genes, namely SOCS1 and ORAI3, which may participate in the preoperative neoadjuvant therapy of ESCC, were screened and veri ed by the GEO database.
In this study, there are two main reasons that Score was used to be the nal evaluation criterion for predicting pCR in the GEO database. Firstly, Score was the sum of the two genes' regression coe cients multiplying by ratios of the expression of the two genes to the abundance of M2 macrophage genes, re ecting the relationship between the two genes and M2 macrophage better. Secondly, the expression of genes as the numerator in Score may be the leading factor of predicting the response of nCRT, not denominator. So it further indicated that the expression of potential genes had better predictive value than the abundance of M2 macrophage.
27 para n-embedded specimens of ESCC were used to verify the results of GEO database. We found no signi cant correlation between CD163 H-score and remission, which was consistent with the results of database analysis. Therefore, the predictive value of Score1 for neoadjuvant therapy can re ect the correlation between SOCS1 expression in the local tumor microenvironment and the curative effect.
Finally, we found that the high expression of SOCS1 is signi cantly correlated with the favorable responsiveness to preoperative neoadjuvant therapy for ESCC.
SOCS1, as a cell signal suppressor, is mainly involved in the regulation of two signaling pathways, namely Janus tyrosine kinase / signal transducer and activator of transcription (JAK/STAT) signaling pathway and toll-like receptors (TLR) signaling pathway [12,13,14,15,16]. In JAK / STAT signaling pathway, SOCS1 can interact with JAK kinase domains, such as Jak1, JAK2, Tyk2, to inhibit the activity of the kinase and affect the downstream STAT protein phosphorylation, leading to the inactivation of the signaling pathway and playing a negative feedback regulation role [17,18,19]. In the TLR signaling pathway, SOCS1 can interact with TLR adaptor proteins, such as MAL and IRAK, leading to their degradation, thereby inhibiting this signaling pathway [16,17,20,21]. As a tumor suppressor, SOCS1 gene promoter methylation leads to the inactivation of SOCS1 gene expression, which occurs in 53-71% of patients with pancreatic cancer and occurs in patients with acute myeloid leukemia and liver cancer [22,23,24]. McCormick found that in the signal pathway of IL-4 activating tyrosine phosphorylation of insulin receptor substrate-2 (IRS-2) to induce M2 macrophage polarization, SOCS1 protein can inhibit M2 macrophage differentiation by mediating IRS-2 ubiquitination and promoting protease degradation of phosphorylated IRS-2 [25]. Whyte and Liu also found that SOCS1 was closely related to the differentiation of M2 macrophages [26,27]. This study shows that SOCS1 can be used as a biomarker for the favorable responsiveness to preoperative neoadjuvant therapy for ESCC. The reason may be that SOCS1 can negatively regulate the differentiation of M2 macrophages in the tumor microenvironment, reduce the in ltration of M2 macrophages in the microenvironment, and thus reduce the immunosuppressive effect of M2 macrophages [28, 29,30].
There are several limitations in this study: (1) The number of samples in training set GSE45670 and validation set GSE104958 is small, and the accuracy of the screened genes needs to be further studied and veri ed; (2) The number of para n-embedded specimens used for validation is small, and the results need to be further con rmed by expanding the sample size; (3) The ESCC specimens were from preoperative endoscopic biopsy, and the samples were taken from the surface of the tumor. So they may not re ect the whole tumor microenvironment. However, compared with the prediction e ciency of tumor microenvironment of preoperative endoscopic biopsy specimens and postoperative resection specimens in previous studies, it is found that the tumor microenvironment of preoperative endoscopic biopsy specimens can also better predict the prognosis and curative effect of patients [31,32]. (4) The function of SOCS1 in M2 macrophage had not been studied, and further research is needed to verify the function.
In conclusion, the high expression of SOCS1 in the local tumor microenvironment is signi cantly related to the favorable responsiveness to neoadjuvant therapy in patients with ESCC, which can be used as a biomarker to predict the responsiveness to neoadjuvant therapy and has certain clinical value.

Declarations
Ethics approval and consent to participate: The ethical approval for obtained from Ethics Committee of the Army Medical Center of the PLA. The committee's reference number was 2020#91. Informed consent was obtained from all participants. All of those studies from databases previously were approved by their respective institutional review boards. The need for informed consent of patients involved in three independent cohorts (TCGA ESCA, GSE45670, GSE104958) and our hospital was waived by the Ethics Committee of the Army Medical Center of the PLA. This study was conducted in accordance with the Declaration of Helsinki.

Consent for publication:
Not applicable Availability of data and materials: The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request. Datasets of TCGA ESCA were downloaded from https://xenabrowser.net/datapages/. GSE45670 and GSE104958 were obtained from the NCBI website (https://www.ncbi.nlm.nih.gov/geo/).

Competing interests:
The authors declare that they have no competing interests. Funding:

Not applicable
Author contributions XY, HX, KX and YY performed design of the study. JL and LL retrieved all the data for ESCC patients information. HX and KX analyzed and interpreted the patient data. KX and YY completed the experiment and were major contributors in writing the manuscript. All authors read and approved the nal manuscript.

Figure 1
The ow chart of analysis.

Figure 3
The correlation between score and M2 macrophage in ltration abundance in GSE104958 and neoadjuvant response. Figure 3a The Score of patients with PR in remission was higher than that of patients with PD+SD in remission, with the median values of -0.19 and 9.56 (P=0.052), respectively.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.