The SOCS3 expression and corresponding clinical follow-up information from different tumors were collected from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/). Normal tissue data from the GTEx database (http://gtexportal.org) was the normal control. The SOCS3 expression difference between tumor tissues and normal tissues was displayed using a boxplot (R package, ggplot2, v 3.4.0). Clinical follow-up information consisted of overall survival (OS). The relationship between SOCS3 and the prognosis was examined using the Cox proportional hazards regression model (R package, survival, coxph, v 3.2-7) and Kaplan–Meier analysis (R package, survival, survfit, v 3.2-7). Median levels of SOCS3 expression divided patients into low and high expression groups in each tumor and performed differentially expressed genes (DEGs) analysis (R package, TCGAbiolinks, v 2.24.3). Subsequently, we performed gene set enrichment analysis (GSEA) of DEGs for KEGG and Hallmark pathway in the pan-cancer, of which the top ten were shown (R package, clusterProfiler, GSEA, v 4.4.4). Circular barplot displayed the result of enrichment analysis (R package, ggplot2, v 3.4.0). Cancer types and corresponding abbreviations were available in the official documentation (https://www.cancer.gov/types).

TIMER (Tumor Immune Estimation Resource) database (https://cistrome.shinyapps.io/timer/) was used to evaluate the abundance of tumor infiltrating immune cells (TIICs) in different cancer types, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells (17). The ratio of the immune-stromal component in the tumor microenvironment (TME), the tumor mutation burden (TMB), and Microsatellite instability (MSI) were calculated using the Sangerbox online platform (http://www.sangerbox.com/tool) (18). ImmuneScore and StromalScore are the scores of each sample’s immune and matrix components in the TME, respectively, while ESTIMATEScore is the sum of the two, representing the combined ratio of the two components in the TME, which means that the higher score matched the higher the corresponding component in the TME. These scores were also calculated using the Sangerbox online platform. The correlation analysis between EXO5 and immune-gene signatures was performed by the TIGER (Tumor Immunotherapy Gene Expression Resource) database (http://tiger.canceromics.org/#/) (19).

We collected colon cancer patients who had lung metastasis at the First Affiliated Hospital of Zhejiang University from 2010 to 2020. All patients were histopathologically diagnosed by two pathologists and the attending physician. Clinical information was gathered, including gender, age at diagnosis and stage. The Tumor-Node-Metastasis (TNM) -based staging of colon cancers was defined by the new 8th editions of the relevant Union for International Cancer Control (UICC) and American Joint Committee on Cancer (AJCC) publications. Patients without clinical information and pathological section were excluded. The gene expression data and corresponding clinical information of GSE68468, GSE41258, GSE22834, and GSE6988 data sets were collected from GEO data sets (https://www.ncbi.nlm.nih.gov/geo/).

A total of 32 colon cancer patients who had lung metastasis were included in this study. Normal intestinal mucosal samples from normal adjacent tissue were used as control samples. The CD68, CD163, and SOCS3 status of all primary tumors and metastases and normal intestinal mucosal samples were conducted by IHC staining through standard staining procedures (SOCS3 antibody: ab16030, Abcam, Cambridge, UK; CD68 antibody: ZM0464; ZSGB-BIO, Beijing, China; CD163 antibody: ZM-0428; ZSGB-BIO, Beijing, China). Brownish-yellow cytoplasmic staining was determined to be positive on target tissue cells. Each sample was assessed with a minimum of 500 cells from 5 representative high-power fields (200 ×) and the proportion of positively colored cells was computed. Positive cell density was divided using a four-point method: 1 point for the proportion of positively stained cells < 25%, 2 points for those between 25%-50%, 3 points for those between 50%-75%, and 4 points for those > 75%. The color intensity of positively stained cells was divided according to the three-point method: 1 point for weakly positive cells, 2 points for moderately positive cells, and 3 points for strongly positive cells. The IHC staining score was calculated by multiplying the density value and color intensity of positive cells, and was determined by the four-point method: 0-1, negative; 2-4, weak (+); 5-8, moderate (++); and 9-12, strong (+++). Two pathologists who were blind to the patient’s clinical characteristics underwent independent IHC staining assessment. All patients were divided into two groups based on IHC staining score: low expression (negative and +) and high expression (++ and +++).

DEGs between the low SOCS3 expression group and the high SOCS3 expression group in colon primary tumor and lung metastasis were compared with |log₂FC| > 1 and P < 0.05, and the shared DEGs were displayed through the heat map (R package, pheatmap, v 1.0.12). STRING database (https://string-db.org/) (20) was used to conduct the relationship between DEGs, and the PPI network was displayed using the Cytoscape 3.7.1 software. DAVID (https://david.ncifcrf.gov/tools.jsp) was used to perform the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis for DEGs (21). The Venn diagram and enrichment factor bubble chart were performed using the OmicStudio tools at https://www.omicstudio.cn/tool.

Pearson’s chi-squared test and student’s t-test were used to evaluate the association of CD68, CD163, and SOCS3 status with gender, age at diagnosis, and stage. The student’s t-test was used to evaluate the differences in SOCS3 expression of the GEO data set. A P-value <0.05 was considered statistically significant. Data management and analysis were performed using R version 4.2, GraphPad Prism 8.0 Software (GraphPadSoftware, Inc., SanDiego, CA) and SPSS version 19.0 (SPSS Inc., Chicago, IL, USA). All methods were carried out by relevant guidelines.

To figure out the difference in SOCS3 expression between normal tissue and tumor tissue in different cancer types, we explored the expression data in TCGA and GTEx databases. The results indicated that urologic neoplasms (including BLCA, KICH, KIRP) and respiratory tumors (including LUAD and LUSC), and ACC, LAML, and THCA were more inclined to low SOCS3 expression, whereas reproductive system neoplasms (including OV, TGCT, UCS) and gastrointestinal cancer (including ESCA, STAD, and PAAD) were more likely to high SOCS3 expression. In addition, the expression of SOCS3 was significantly downregulated in PRAD and BRCA, which belong to reproductive system neoplasms, and COAD and LIHC, members of gastrointestinal cancer (Figure 1A).

According to the Cox proportional hazards regression model, the SOCS3 expression was less obviously associated with prognosis and tumor progression. High SOCS3 expression patients harbored a poor prognosis in ESCA, GBM, HNSC, KIRC, LGG, STAD, and UVM, whereas harbored a better prognosis in BRCA (Figure 1B). Subsequently, we divided patients into low expression group and high expression group by the median of SOCS3 expression, and performed the Kaplan–Meier analysis. The prognosis of the low SOCS3 expression group in LGG and KIRC was significantly better than that of the high SOCS3 expression group, whereas the opposite was in BRCA (Figure 1C). However, over time, the prognosis of the two groups converged.

To shed further light on the potential function of SOCS3 in pan-cancer, we divided patients into low and high expression groups in each tumor according to median levels of SOCS3 expression, and performed DEGs analysis. Subsequently, we performed GSEA through the GO functional annotation, KEGG pathway and the Hallmark pathway and summarized the pathways that were common to different tumors in the enrichment analysis results (Figure 1D). For GO enrichment analysis, the high SOCS3 expression group was chiefly enriched in developmental regulation, regulation of cell differentiation, and immune response, whereas the low expression group was in organonitrogen compound biosynthetic process and oxidative phosphorylation. The results of the KEGG enrichment analysis revealed that patients with high SOCS3 expression mainly exhibited enrichment of the ECM receptor interaction, immune-related disease, and antigen processing and presentation. Those patients showed enrichment of epithelial mesenchymal transition, TNFα signaling via NFκb, and interferon response for the Hallmark pathway. Therefore, we speculated that there was a possible association between SOCS3 and the immune response and microenvironment in different tumors.

Based on the above discovery, we explored the relationship between SOCS3 expression and immune cell infiltration in different tumors through the TIMER database, including B cell, CD4+ T cell, CD8+ T cell, neutrophil, macrophage, and dendritic. It is observed that a large scale of tumors was positively related to immune cell infiltration, in which the four most significantly associated tumors were BRCA, LIHC, PRAD, and COAD (Figure 2A). Through the ESTIMATE algorithm, immune scores were calculated in different tumors (Figure 2B), for which the top three significantly associated tumors were BLCA (R= 0.54, P< 0.001), BRCA (R= 0.43, P< 0.001), and COAD (R= 0.69, P< 0.001).

Besides, TMB and MSI, two biomarkers associated with the immunotherapy response, were analyzed concerning the SOCS3 expression in different tumors. The expression level of SOCS3 positively and significantly correlated with TMB in several tumors, including COAD, LAML, LGG, SARC, and THYM, whereas negatively and significantly correlated with that in UCEC, ESCA, LIHC, PAAD, PRAD, and THCA (Figure 3A). COAD exhibited positive correlations between SOCS3 expression and MSI, and DLBC, HNSC, STAD, and UCEC exhibited negative correlations (Figure 3B). Subsequently, we analyzed the relationship between SOCS3 expression and 11 published immune-gene signatures, including T cell-inflamed GEP, CAF, TAM M2, IFNG, CD8, CD274, TLS, TLS-melanoma, T cell Dysfunction, T cell exclusion, and MDSC. We found that the SOCS3 expression in the majority of tumors positively correlated with immune-gene signatures, except for TAM M2 (Figure 3C).

The function of SOCS3 is mainly focused on immune response, which can regulate macrophage and TH cell functions. Also macrophages in tumors can promote cancer development and metastasis through various mechanisms. Therefore, based on the previous research and the distinct relationship between SOCS3 status and the immune response in COAD, we collected tumor tissue and adjacent normal tissue from 32 colon cancer patients who had lung metastasis, and performed IHC staining with several markers, including CD68, CD163, and SOCS3. CD68 and CD163 are currently the most widely used macrophage marker and are predominantly located in cell membranes and lysosomal membranes. IHC staining indicated that there was a clear aggregation distribution effect on the leading edge of the colon primary tumors and lung metastasis focus for CD68 and CD163 markers, and this effect was more remarkable in lung metastasis than in colon primary tumor (Figures 4A, B). Compared with adjacent normal intestinal mucosal tissue, SOCS3 status in colon primary tumor tissue was slightly weak, mainly with moderate positive staining in the cytoplasm (Figure 4C). Nevertheless, SOCS3 and CD163 staining in lung metastasis showed more intense membrane and cytoplasm staining than in colon primary tumors (Figure 4D).

To search out the relationship between SOCS3 expression and clinical characteristics, we gathered clinical information from these patients during the first visit and calculated the IHC staining score of SOCS3. Patients who had lymph node metastasis and distant metastasis at the initial diagnosis preferred high SOCS3 expression in lung metastasis focus, whereas had no significant difference in colon primary tumors focus (Table 1). Besides, the high IHC staining score of SOCS3 and CD163 mostly correlated with lung metastasis focus (Table 2). Compared with the colon primary tumor, lung metastasis harbored higher CD163 and SOCS3 expression (Table 3). Subsequently, we analyzed the correlation between SOCS3 expression and macrophage infiltration, and found that SOCS3 expression significantly correlated with macrophage markers and high SOCS3 expression was more inclined to high CD163 expression in lung metastasis (Table 4). The results suggested that SOCS3 might be associated with macrophage infiltration and tumor invasion in lung metastasis.

Based on IHC staining results, we investigated transcription data in the GEO database to explore whether the transcript level and protein level have the same characteristics. The most common site of blood metastasis of colon cancer is the liver, followed by the lung, brain and bone. Transcript data sets consisting of the colon primary tumor and metastasis in the GEO database were used, including GSE68468 and GSE41258 (colon primary tumor and lung metastasis), GSE22834 and GSE6988 (colon primary tumor and liver metastasis). The CD68 and CD163 expression levels in metastasis were higher than in the colon primary tumor. The SOCS3 expression of lung metastasis was significantly higher than in the colon primary tumor, whereas that in liver metastasis was the opposite (Figure 5). These results revealed that SOCS3 shows different manifestations in different metastasis.

To explore the molecular mechanisms of SOCS3 in lung metastasis, we divided the colon cancer and lung metastasis patients into two groups in the GSE68468 dataset respectively via the median of SOCS3. Subsequently, DEGs were compared between low SOCS3 expression and high SOCS3 expression in colon primary tumor and lung metastasis with |log₂FC| > 1 and p < 0.05. 367 DEGs were screened out in colon primary tumors and 459 DEGs in lung metastasis, and 51 genes were shared by two sets of DEGs (Figure 6A). There were 353 upregulated genes and 4 downregulated genes in the colon primary tumor, and 247 upregulated genes and 212 downregulated genes in lung metastasis (Figure 6B). Interestingly, four genes that were originally highly expressed in the high SOCS3 expression group in colon cancer became low expressed after lung metastasis, including TTC40, PCDHA2, SP3P and ZNF471. Through STRING and Cytoscape, we conducted the interplay network of 51 shared DEGs and screened out two genes as hub genes, including CCL2 and PTGS2 (Figure 6C). Furthermore, the enrichment analysis of GO and KEGG pathways was performed to explore the underlying interplay of the exceptional DEGs in lung metastasis. For GO analysis, 20 biological processes (BPs) were enriched, such as immune response, MHC class II protein related process, antigen processing, T cell activation, and inflammatory response (Figure 7A). Regarding the KEGG pathway analysis, Rheumatoid arthritis, Th17 cell differentiation, and Graft−versus−host disease were enriched (Figure 7B). These results showed that immune responses and regulations were significantly enriched in both GO and KEGG analysis.