MEOX2 Serves as a Novel Biomarker Associated with Macrophage Inltration in Esophageal Squamous Cell Carcinoma

Background Esophageal squamous cell carcinoma (ESCC) is a kind of digestive system malignant tumor with high morbidity and mortality worldwide. With the rise of immunotherapy applied in cancer treatment, immune-related mechanisms in ESCC and other digestive system carcinomas are in urgent need of being detected. In our study, single-sample gene set enrichment analysis (ssGSEA) was performed at rst to analyze the expression prole downloaded from NCBI Gene Expression Omnibus (GEO) database. Then via a series of bioinformatic analyses, including Mann-Whitney test, weighted gene co-expression network analysis (WGCNA), functional enrichment analysis and differentially expressed genes (DEGs) analysis, we identied targeting immunocyte and related genes. Finally, we validated the results in TIMER database. activation; GO:0050900: leukocyte migration; ko04640: Hematopoietic cell lineage; GO:0009617: response to bacterium; and GO:0070661: proliferation. TIMER database to verify the correlation between MEOX2 expression level and degree of immunocytes inltration in esophageal carcinoma As shown in gure, MEOX2 expression levels have a highly positive correlation with inltrating levels of macrophage (r = 0.449, P = 2.52e-10). Then we further assessed the correlations between MEOX2 expression level and degree of immunocytes inltration in other digestive system carcinomas, including stomach adenocarcinoma (STAD), pancreatic adenocarcinoma (PAAD), colon adenocarcinoma (COAD) and rectum adenocarcinoma (READ) (Figure 6B-6E). The results indicated that MEOX2 expression levels had a signicantly positive correlation with inltrating levels of macrophages in all kinds of carcinomas mentioned above, including STAD (r = 0.634, P = 4.53e-43), PAAD (r = 0.629, P = 3.24e-20), COAD (r = 0.646, P = 5.40e-49) and READ (r = 0.495, P = 6.07e-10).


Abstract Background
Esophageal squamous cell carcinoma (ESCC) is a kind of digestive system malignant tumor with high morbidity and mortality worldwide. With the rise of immunotherapy applied in cancer treatment, immunerelated mechanisms in ESCC and other digestive system carcinomas are in urgent need of being detected.

Methods
In our study, single-sample gene set enrichment analysis (ssGSEA) was performed at rst to analyze the expression pro le downloaded from NCBI Gene Expression Omnibus (GEO) database. Then via a series of bioinformatic analyses, including Mann-Whitney test, weighted gene co-expression network analysis (WGCNA), functional enrichment analysis and differentially expressed genes (DEGs) analysis, we identi ed targeting immunocyte and related genes. Finally, we validated the results in TIMER database.

Results
Our analyses showed that macrophage in ltrating level is obviously higher in advanced stages in ESCC compared with other types of immunocytes. MEOX2 was detected as a biomarker correlated with macrophage in ltration in ESCC and other types of digestive system carcinomas. And MEOX2 expression was strongly associated with mRNA expression of colony-stimulating factor 1 (CSF-1) and CSF-1 receptor (CSF-1R) in these kinds of carcinomas.

Conclusion
We speculated that MEOX2 could facilitate macrophage in ltration via CSF-1/CSF-1R signaling in ESCC and other kinds of digestive system carcinomas, which might serve as a novel target in prospective tumor immunotherapy.

Background
As a common type of upper digestive tract malignant tumor, esophageal carcinoma (EC) is the sixth most common cause of cancer-related death around the world, and account for 509,000 cancer-related deaths each year.(1) Esophageal carcinoma is composed of two main histological types, esophageal adenocarcinoma (EAC) and esophageal squamous cell carcinoma (ESCC), of which squamous-cell carcinoma makes up about 90% of cases worldwide.(2) Although large amounts of breakthroughs have been achieved in diagnosis, screening, surgical resection and chemoradiotherapy of ESCC, the prognosis remain unsatisfactory.(3) Now the immunotherapy has gradually become a hot topic in oncology treatment, and its application in melanoma, pulmonary carcinoma and many other malignancies has been proven a great success.(4) However, immunotherapy in ESCC treatment still remain limited due to absence of reliable biomarkers related to immunocytes in ltrated in tumor microenvironment (TME) of ESCC. Therefore, it is urgent to detect novel biomarkers of immune interaction with ESCC and identify new targets for immune-related therapy.
Previous studies have reported parts of immune microenvironment characteristics of ESCC. Esophageal carcinoma was considered to have an adequate mutational load for checkpoint inhibitor treatment, which opened up possibilities for immunotherapy in this disease.(5) Schumacher et al. found that the CD8(+) T cell in ltration in esophageal carcinomas was a favorable prognostic factor and had potential therapeutic value.(6) Nishimura et al. detected that mature LAMP-3 (lysosome-associated membrane glycoprotein 3) DCs (Dendritic cells) was correlated with increasing in ltrating CD8(+) T cells in ESCC, and that were identi ed to be associated with a better prognosis. (7) But in consideration of complicated factors implicated in the regulation of tumor immunity, more studies are required to complete the immune atlas in ESCC. Currently, with the development of bioinformatics, a great deal of public data containing gene expression and clinicopathologic information can be utilized to explore the tumor microenvironment of ESCC. Therefore, on the basis of public data and bioinformatic analysis, our study aimed to explore the association between in ltrating immunocytes and clinicopathologic characteristics, and to detect novel immune-relevant genes. In this study, via singlesample gene set enrichment analysis (ssGSEA), weighted gene co-expression network analysis (WGCNA), differentially expressed genes analysis (DEGs) and online database analysis, we identi ed a novel biomarker associated with immunocytes in ltrating in ESCC. The detailed work ow is shown in Fig. 1.

Data Collection and Preprocessing
The gene expression pro le of GSE69925 was searched in NCBI Gene Expression Omnibus (GEO) with "esophageal squamous cell cancer" and "Homo sapiens" as keywords. Probes in this series were annotated by gene symbol information from soft formatted family les, the probe annotated by multiple genes was deleted, and the average value was recorded as the nal expression value for the gene matched by multiple probes. Implementation of Single-Sample Gene Set Enrichment Analysis (ssGSEA) Immune-related gene sets were obtained from the cancer immunome database (TCIA).(9) The literature from TCIA contains 28 immunocyte-related gene sets, including innate immunity and adaptive immunity. The innate immunity consists of activated dendritic cell, immature dendritic cell, plasmacytoid dendritic Page 4/17 cell, CD56dim natural killer cell, CD56bright natural killer cell, eosinophil, mast cell, macrophage, MDSC, monocyte, natural killer cell, natural killer T cell and neutrophil. While the adaptive immunity comprises of activated CD4 T cell, activated CD8 T cell, central memory CD4 T cell, central memory CD8 T cell, effector memory CD4 T cell, effector memory CD8 T cell, gamma delta T cell, activated B cell, immature B cell, memory B cell, T follicular helper cell, regulatory T cell, type 1 T helper cell, type 2 T helper cell and type 17 T helper cell.
Related in ltration levels for 28 immunocytes were quanti ed by the ssGSEA algorithm in R package GSVA, which could apply gene signatures expressed by immunocytes to every cancer samples.(10) Then the R package ESTIMATE was utilized to calculate the immune score, stromal score, tumor purity and estimate score of every tumor sample. (11) Identi cation of Differentially In ltrated Immunocytes between Different Stages Based on T stages of samples containing tumor stage information, 40 samples were classi ed into two groups, T2+T3 and T4. Then we could identify immunocytes in ltrating differentially between two groups, which was analyzed by the Mann-Whitney test, and differences were considered statistically signi cant when P Value<0.05. The immunocyte with the greatest difference and clinical signi cance was selected as the target immunocyte for next analysis.
Identi cation of In ltrating Immunocyte Related Gene List by Weighted Gene Co-expression Network Analysis (WGCNA) The expressing pro ling of all 274 samples was applied to this section. In this section, 5,000 genes with the highest average expression were selected, and the gene expression was calculated using the WGCNA algorithm.(12) First, the soft threshold for network construction was selected by calculating the scale-free topology t index for several powers. Then, the scale-free network was constructed and the module eigengene (ME) was calculated. Each module contained a gene list with the most appropriate accounting of gene expression pro le. And all 274 samples were divided into highly-in ltrating group and lowlyin ltrating group according to in ltrating scores of the targeting immunocyte, which was utilized as the clinicopathologic features for next step. Finally, the correlation between ME and clinicopathologic features in each module was calculated. According to module-trait relationships, the module displaying the most obvious correlation was considered a key module and selected for further calculation. (13) Functional Enrichment Analysis Metascape is a meta-analysis website integrating gene annotation and functional enrichment analysis. In this section, by using the "Express Analysis" module in Metascape, we explored the enriched biological processes and pathways of the gene list contained in the key module. (14) Differentially Expressed Genes (DEGs) and Correlation Analysis By DEGs analysis of gene list contained in the key module between highly-in ltrating and lowly-in ltrating group, we could further screen targeting genes associated with immunocyte in ltrating. DEGs were analyzed using package limma of R software with adjusted P Value<0.05 and |logFC| >1. (15) Spearman's correlation test was performed to evaluate potential correlations between mRNA expression of target genes and immunocyte in ltration levels. (16) TIMER Database Analysis Gene expression levels against tumor purity is a critical factor that have an effect on the analysis of immunocyte in ltrating level and is shown in the left-most panel. (11) The correlations between mRNA expressions of target genes and that of related colony-stimulating factors (CSFs) were also analyzed in TIMER.

Data Information
GSE69925 was submitted by Japan National Cancer Center Research Institute, and was constitutive of 274 pretreatment biopsy specimens of esophageal squamous cell carcinoma, which were performed by comprehensive gene expression pro ling. After preprocessing, 21,654 genes were annotated and retained for following analysis.

In ltrated Immunocytes Relevant to Different Stages in ESCC
The immune in ltration levels of 40 samples with stage information is visualized as a heat map ( Figure  2A). Between two groups (T2+T3 and T4), our analysis showed that 7 types of immunocytes had signi cantly different in ltrating levels, including central memory CD8 T cell, gamma delta T cell, immature B cell, macrophage, monocyte, natural killer T cell and type 2 T helper cell. Among these types of immunocytes, macrophage had a greater statistic difference (P<0.001, Figure 2B) and had been reported to participate in the tumor promoting activities in esophageal carcinoma. (18)  Based on 5,000 genes with the highest average expression in 274 samples, gene co-expression modules were analyzed by WGCNA algorithm. In our analysis, when power value is 3, the scale-free co-expression network can be achieved ( Figure 3B). Then 6 modules containing the 5,000 genes above-mentioned were constructed. All the modules are marked with different colors and shown in gure 3C. Using macrophage in ltrating levels as clinicopathologic features, we could plot module-feature relationships ( Figure 3D).
From the gure, we could detect that high macrophage in ltrating group (Macro-H) was strongly associated with the yellow module (r = 0.43, P = 1E-13). Gene signi cance of 125 genes contained in the yellow module versus module membership is depicted in a scatter plot ( Figure 3E). Therefore, the gene list contained in the yellow module was considered to be signi cant for macrophage in ltration in ESCC and selected as the hub gene list.

Functional Enrichment Analysis of Hub Gene List
Biological processes of the gene list were explored utilizing GO and KEGG analysis in Metascape ( Figure  4). As shown in this gure, the top 5 enrichment items are: GO:0046649: lymphocyte activation; GO:0050900: leukocyte migration; ko04640: Hematopoietic cell lineage; GO:0009617: response to bacterium; and GO:0070661: leukocyte proliferation. Results of enrichment analysis indicated that the hub gene list was highly connected with immune system functions.

Identi cation of Macrophage In ltration Related Genes
The hub gene list was further screened by differentially expressed genes (DEGs) analysis between high macrophage in ltration group and low macrophage in ltration group. In this step, 29 genes expressing differentially were extracted via limma package in R (Table 1, Figure 5A). Considering that the selected gene list was strongly correlated with high level of macrophage in ltration, every DEGs was upregulated gene (adj.P.Val<0.05, logFC>1).
The DEGs contained some genes that were members of TCIA immune gene sets used for ssGSEA, which were excluded rstly. Then, the genes that expressed speci cally in immune cells were eliminated. Finally, we focused on MEOX2, which was also known as GAX and MOX2 according to genebank in NCBI. (19) Expression level of MEOX2 in two groups was visualized as a violin plot in gure 5B. The correlation between MEOX2 expression and macrophage in ltrating level was calculated by Spearman's correlation test, which is shown in gure 5C (r=0.31, P=1.16e-7).  Colony-stimulating factor 1 (CSF-1) is known as essential for proliferation and differentiation of macrophages. Then we explored correlations between mRNA expressions of MEOX2 and mRNA expressions of colony-stimulating factor 1 (CSF-1) and colony-stimulating factor 1 receptor (CSF-1R) in TIMER database. The results showed that MEOX2 expression had positive correlations with mRNA expressions of CSF-1 (r = 0.336, P = 2.84e-6) and CSF-1R (r = 0.475, P = 8.07e-12) in esophageal carcinoma ( Figure 7). Further analysis indicated that in all the other digestive system carcinomas MEOX2 expression had positive correlations with mRNA expressions of CSF-1 and CSF-1R as well, correlations in these kinds of carcinomas were even more obvious than that in esophageal carcinoma (Figure 7).

Discussion
In this study, we found that macrophage in ltrating level is obviously higher in advanced stages in ESCC. Then by a series of bioinformatic analysis, MEOX2 was detected as a biomarker with a positive correlation with macrophage in ltration in ESCC and other types of digestive system carcinomas. As far as we know, in our study MEOX2 is rst reported as a biomarker associating with immunocytes in ltration in carcinoma.
Tumor associated macrophages (TAMs) are de ned as macrophages settled in the tumor tissues, which are essential members of immunocytes that constitute tumor microenvironment (TME) and play complicated roles in tumorigenesis, tumor development and metastasis. (20) In initial stage of tumor, attracted by tumor cell derived colony-stimulating factors (CSFs), macrophages from surrounding tissues would in ltrate in the tumor tissue and construct supporting microenvironment for tumor cells. (21) Then macrophages secreted various kinds of pro-tumorigenic cytokines, which could promote tumor cell proliferation, inhibit anti-tumor immunity, and upregulate angiogenesis in tumor tissues. (22)(23)(24)(25) With the progress of tumor, some collagenases like MMP2 and MMP9 derived from macrophages could degrade basement membrane and promote angiogenesis, which played indispensable roles in tumor invasion and metastasis. (26) As an indispensable factor for macrophage differentiation and survival, colony-stimulating factor 1 (CSF-1) has been detected to be largely expressed in several kinds of carcinoma cells. (22,27,28) The CSF-1 receptor (CSF-1R) belongs to transmembrane tyrosine kinase receptor and is specially expressed in monocytic lineage. (27) And CSF-1/ CSF-1R mediated signaling was reported to be associated with worse prognosis in a variety of carcinomas, including breast and hepatic cancer. (29,30) Therefore, blockades targeting CSF-1/ CSF-1R has drawn more and more attention recently, such as humanized monoclonal antibody RG7155 and small molecular compound PLX3397, which has displayed e cacy in several preclinical models. (31,32) Role of TAMs in esophageal carcinoma has also been reported in many literatures. Manabu Shigeoka et al. found that high TAM density was correlated with poor prognosis in ESCC. (33) In our study, macrophage in ltrating level in stage T4 is apparently higher than that in stage T2-3, which is consistent with the previous study. Taisuke Yagi et al. detected that TAMs could upregulate programmed death ligand 1 (PD-L1) expression in esophageal carcinoma cells, which would be likely to have clinical implications on immune checkpoint therapy. (34) In addition, TAMs were identi ed as a biomarker for predicting chemotherapy response in esophageal carcinoma as well. (35) Mesenchyme Homeobox 2 (MEOX2) is a protein coding gene, which is also known as GAX (growth arrest homeobox) or MOX2. (36) Encoding production of MEOX2 is a homeodomain transcription factor, which has traditionally been regarded as essential for bone and muscle growth in embryos. (37) In recent years, MEOX2 has been detected that it could suppress cell proliferation and EMT (epithelial-to-mesenchymal transition) in cell models. (38,39) In addition, MEOX2 has also been reported to paly vital roles in platinum-based drug resistance and EGFR-TKI-based (epidermal growth factor receptor-tyrosine kinase inhibitor) treatment response in non-small cell lung cancer (NSCLC) patients, which means poor prognoses in these patients. (40) Whereas the role of MEOX2 in immune in ltration of TME has never been described.
In our study, mRNA expression of MEOX2 was positively correlated with macrophage in ltrating level in esophageal squamous cell carcinoma, which was con rmed by the following analysis in TIMER database. Such correlations could be detected in other kinds of digestive system carcinomas as well, which were even more obvious. Further analysis told us that MEOX2 expression was strongly associated with mRNA expression of CSF-1 and CSF-1R in esophageal carcinoma and all the other digestive system carcinomas. Collectively, we speculated that MEOX2 could promote macrophage in ltrating into tumor microenvironment via CSF-1/CSF-1R signaling in esophageal carcinoma and other digestive system carcinomas.
We acknowledge that there are several limitations in this study. Histological type of samples contained in GSE69925 is squamous cell carcinoma, whereas the majority of esophageal samples utilized in TIMER database are adenocarcinoma. In addition, our study was based on bioinformatic analysis of expression pro le from public databases, further biological experiments were required to con rm the results.

Conclusions
In summary, our study provides a novel perspective into the function of MEOX2 in macrophage in ltration in ESCC and other kinds of digestive system carcinomas. In carcinomas mentioned above, MEOX2 expression was not only associated with macrophage in ltrating level but also positively correlated with mRNA expression of CSF-1 and CSF-1R. Therefore, we supposed that MEOX2 could facilitate macrophage in ltration via CSF-1/ CSF-1R signaling in these kinds of carcinomas.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The datasets analyzed during the current study are available in the NCBI Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/). cancer drug resistance, overall survival and therapy prognosis in lung cancer patients. Oncotarget