Tumor-associated macrophage polarization promotes the progression of esophageal carcinoma

Distribution of TIICs in the ESCA microenvironment

To explore the distribution of TIICs in the ESCA microenvironment, we analyzed ESCA samples in The Cancer Genome Atlas (TCGA). We found that 70-80% of the TIICs in the ESCA microenvironment were macrophages and T cells with different phenotypes (Figure 1).

Identification of DEGs between macrophages of different phenotypes

Our previous research indicated that a high density of M2 macrophages in the TME was associated with the occurrence and progression of ESCC, and with a poor prognosis in ESCC patients [6]. Thus, we used various macrophage datasets from the Gene Expression Omnibus (GEO) to identify genes that were differentially expressed between macrophages of different phenotypes. As shown in the Venn diagram in Figure 2A, we identified 25 consistently DEGs between M0 macrophages and M2 macrophages in three datasets (GSE57614, GSE36537 and GSE5099). Coincidentally, all 25 DEGs were highly expressed in samples of M2 macrophages. Next, we used cBioPortal to assess the expression of these 25 DEGs in ESCA tissues. The results revealed that some of the DEGs were highly expressed in ESCA tissues (Figure 2B). In addition, we identified 91 consistently DEGs between M1 and M2 macrophages in four datasets (GSE57614, GSE36537, GSE5099 and GSE95405) (Figure 2C). A volcano plot of the 91 DEGs demonstrated that tumor necrosis factor receptor superfamily member 11a (TNFRSF11A) was the most upregulated gene and C-X-C motif chemokine ligand 10 (CXCL10) was the most downregulated gene in M2 macrophages compared to M1 macrophages (Figure 2D).

FGL2 expression correlates with the level of immune cell infiltration in ESCA

To investigate whether the DEGs between M0 macrophages and M2 macrophages were associated with the prognosis of esophageal cancer patients, we used the Kaplan-Meier Plotter online platform to assess the survival of ESCC patients based on their expression of the 25 DEGs. We found that FGL2, exoribonuclease 1 (ERI1) and WNT5B levels were associated with the overall survival (OS) and relapse-free survival (RFS) of ESCC patients. Specifically, higher FGL2 expression was marginally associated with a poorer prognosis in ESCC patients (OS hazard ratio [HR] = 2.57, 95% confidence interval [CI] = 1.05 to 6.28, P = 0.033; RFS HR = 3.78, 95% CI = 1.43 to 9.98, P = 0.0039) (Figure 3A, 3B). However, lower ERI1 expression was associated with a poorer prognosis (OS HR = 0.22, 95% CI = 0.08 to 0.59, P = 0.0011; RFS HR = 0.32, 95% CI = 0.12 to 0.84, P = 0.015) (Figure 3C, 3D). Likewise, lower WNT5B expression was associated with a poorer prognosis (OS HR = 0.41, 95% CI = 0.19 to 0.93, P = 0.027; RFS HR = 4.32, 95% CI = 0.99 to 18.93, P = 0.034) (Figure 3E, 3F). Information about the functions of FGL2, ERI1 and WNT5B can be found in Table 1.

Next, we used the Tumor IMmune Estimation Resource (TIMER) to assess whether FGL2, ERI1 or WNT5B expression correlated with infiltrating immune cell levels in ESCA. Tumor purity significantly influences the genomic analysis of immune cell infiltration in clinical tumor samples [20]. As shown in Figure 4, FGL2 expression correlated negatively with ESCA tumor purity and positively with the infiltrating levels of B cells (r = 0.331, P = 5.85e-06), CD8+ T cells (r = 0.147, P = 4.93e-02), CD4+ T cells (r = 0.323, P = 1.07e-05), macrophages (r = 0.559, P = 3.37e-16), neutrophils (r = 0.337, P = 3.66e-06) and DCs (r = 0.268, P = 2.73e-04). However, ERI1 expression did not correlate significantly with ESCA tumor purity, and only correlated weakly with the infiltrating levels of B cells (r = 0.243, P = 10.4e-03), CD8+ T cells (r = 0.179, P = 1.59e-02) and DCs (r = -0.185, P = 1.28e-02). Similarly, WNT5B expression did not correlate significantly with ESCA tumor purity, and only correlated weakly with the infiltrating levels of macrophages (r = 0.239, P = 1.21e-03). These findings strongly suggested that FGL2 promotes the infiltration of immune cells, especially macrophages, in ESCA.

Correlation analysis between FGL2 expression and immune cell marker sets

To explore the relationship between FGL2 and various infiltrating immune cells, we used the TIMER and Gene Expression Profiling Interactive Analysis (GEPIA) databases to examine the correlations between the levels of FGL2 and markers of various immune cells (monocytes, TAMs, M1 and M2 macrophages, B cells, T cells [general], CD8+ T cells, neutrophils, natural killer cells and DCs) and different functional T cells (Tregs, exhausted T cells, Th1, Th2, Tfh and Th17 cells) in ESCA samples (Table 2). After adjusting for purity, we found that FGL2 levels correlated significantly with the marker levels for most immune cells and different T cells in ESCA samples (Table 2). Higher FGL2 expression in M2 macrophages was associated with greater DC and CD8+ T cell infiltration in ESCA tissues, which was notable because DCs can promote tumor metastasis by reducing CD8+ T cell cytotoxicity and increasing Treg numbers [21]. We also found significant correlations between the levels of FGL2 and markers of exhausted T cells and Tregs, such as programmed cell death 1 (PD-1), cytotoxic T-lymphocyte associated protein 4 (CTLA4), lymphocyte activation gene 3 (LAG3), T-cell immunoglobulin mucin 3 (TIM-3), forkhead box P3 (FOXP3), C-C motif chemokine receptor 8 (CCR8) and signal transducer and activator of transcription 5B (STAT5B) (Table 2). These results confirmed that FGL2 expression correlates with infiltrating immune cell levels in ESCA.

Next, we constructed scatterplots to depict the correlations between the levels of FGL2 and markers of monocytes (e.g., CD86, CSF1R and CD14) and different macrophage phenotypes (e.g., CD68, CCL2 and IL-10 of TAMs; IRF5, NOS2 and PTGS2 of M1 macrophages; CD163, VSIG4 and MS4A4A of M2 macrophages) in ESCA samples, as shown in Figure 5. Interestingly, the results from the TIMER database demonstrated that the levels of all the monocyte, TAM and M2 macrophage marker sets correlated strongly with FGL2 levels (Figure 5A, 5B, 5D), while M1 macrophage marker levels did not correlate with FGL2 levels in ESCA (Figure 5C). Similar correlations between FGL2 expression and monocyte, TAM and M2 macrophage marker levels in ESCA samples were found using the GEPIA database (Table 3). These data verified the results from the GEO database, indicating that FGL2 may promote the polarization of TAMs into M2 macrophages in ESCA.

Together, the above findings suggested that M2 macrophages may regulate immune cell infiltration by secreting the immunosuppressive factor FGL2, thereby producing a microenvironment that promotes the occurrence and progression of ESCA.

FGL2 is crucial in promoting the initiation and progression of ESCA

The upregulation of FGL2 has been reported in different cancers [22]. Thus, we evaluated FGL2 expression in esophageal cancer tissues from TCGA and esophageal cancer cell lines from the Cancer Cell Line Encyclopedia database. As shown in Supplementary Figure 1, FGL2 expression level in esophageal cancer tissues was higher than normal tissues. However, FGL2 expression was low in esophageal cancer cell lines. We speculate that FGL2 is highly expressed in the tumor stroma, and that M2 macrophages in the tumor stroma secrete FGL2 to alter the TME.

To assess the association between the levels of FGL2 and other genes in ESCA tissues, we analyzed mRNA sequencing data from 184 ESCA patients using the function module of LinkedOmics. As shown in the volcano plot in Figure 6A, 4,078 genes exhibited significant positive correlations with FGL2, whereas 1,177 genes exhibited significant negative correlations with FGL2 (false discovery rate [FDR] < 0.01). The 50 most significant gene sets that correlated positively and negatively with FGL2 are shown in the heatmaps in Figures 6B, 6C. These results revealed the widespread impact of FGL2 on the transcriptome.

To better understand whether FGL2 is a crucial determinant of the TME and tumor metastasis in esophageal cancer, we compared the cytokine levels among EC109/9706 (ESCC) cells cultured alone, co-cultured with M0 macrophages or co-cultured with M2 macrophages (Figure 7A). We considered a 20% change in cytokine expression to be statistically significant. The levels of interleukin (IL)-10, matrix metalloproteinase 9 (MMP9), C-C motif chemokine ligand 5 (CCL5), TIM-3, IL-13, vascular cell adhesion molecule 1 (VCAM1), macrophage colony-stimulating factor (M-CSF) and fibroblast growth factor 7 (FGF-7) were significantly greater in EC109/9706 cells co-cultured with M2 macrophages than in EC109/9706 cells cultured alone or co-cultured with M0 macrophages (Figure 7B–7E). These results indicated that M2 macrophages promote the secretion of these cytokines in the TME, thereby inducing the occurrence and progression of ESCA.

Next, we used cBioPortal to evaluate the correlations between the levels of FGL2 and the cancer-promoting cytokines measured above (IL-10, MMP9, CCL5, TIM-3, IL-13, VCAM1, M-CSF and FGF-7) in ESCA tissues. As shown in the scatterplots in Figure 8, FGL2 expression exhibited a strong positive correlation with the expression of each of these cytokines. The positive correlation between FGL2 and TIM-3 (HAVCR2), a crucial promoter of T cell exhaustion [23], suggested that high FGL2 expression accompanies TIM-3-induced T cell exhaustion. We also assessed the correlation of FGL2 expression with PD-1, CTLA4, FOXP3 and CCR8 expression in ESCA tissues in cBioPortal. FGL2 levels correlated strongly and positively with the levels of PD-1 and CTLA4, which are induced during T cell exhaustion (Supplementary Figure 2), consistent with the results in Table 2.

These results confirmed that increased FGL2 expression in M2 macrophages alters the TME in ESCA. The correlation between the expression of FGL2 and the presence of specific infiltrating immune cells in ESCA suggested that FGL2 is a vital contributor to immune escape and immunosuppression in the esophageal cancer microenvironment.

Reversing TAM polarization may provide new targets for tumor immunotherapy

We previously demonstrated that TAMs directly promote the survival, invasion and metastasis of esophageal cancer cells [24]. To identify an effective immunotherapeutic target for ESCA, we used the WEB- based GEne SeT AnaLysis Toolkit (WebGestalt) to conduct an enrichment analysis of the GO biological functions (Biological Process, Molecular Function and Cellular Component) and KEGG pathways of the 91 DEGs identified between M1 and M2 macrophages above. The GO analysis indicated that the DEGs were primarily enriched in the Immune response Type I interferon signaling pathway, Immune effector process, Tumor necrosis factor receptor binding, Major histocompatibility complex (MHC) protein binding, STAT family protein binding, MHC class I peptide loading complex, Cytosol, Cytoplasm, etc. (Figure 9A). Furthermore, the KEGG pathway enrichment analysis demonstrated that these DEGs were enriched in the Toll-like receptor (TLR) signaling pathway, Tumor necrosis factor signaling pathway, RIG-l-like receptor signaling pathway, Pertussis, NOD-like receptor signaling pathway, Nuclear factor κB signaling pathway, Necroptosis, Measles, Influenza A, Human papillomavirus infection, etc. (Figure 9B). These results suggested potential treatment targets involved in TAM polarization that could be used to improve anti-tumor immunity.

We also performed a STRING database screening to analyze the functional interactions between the proteins encoded by the 91 DEGs. We used Cytoscape software to visualize closely connected regions in the protein-protein interaction (PPI) network, and used the MCODE plugin to identify the top hub genes constituting the most closely connected module. The PPI network included 54 nodes and 422 edges (Figure 9C). The top 30 genes among the 54 hub genes are shown in Figure 9D. These key genes could be considered during the development of novel molecular drugs for TAM-related immunotherapies to treat ESCA.

FGL2 is already considered a potential molecular target for glioblastoma treatment [25]. Therefore, we performed a gene set enrichment analysis (GSEA) to identify functions or pathways that varied between ESCA tissues with low and high FGL2 levels. We found that the mitogen-activated protein kinase (MAPK), Janus kinase (JAK)/STAT, TLR and other immune-related pathways were enriched in ESCA tissues with high FGL2 levels (Supplementary Figure 3). Our results indicated that inhibiting FGL2 could potentially weaken the immunosuppressive activity of ESCA tumors to prevent their progression. A schematic diagram depicting the potential function of FGL2 in ESCA is shown in Figure 10.

Evaluation of TIICs in ESCA

In this study, we used expression data from TCGA ESCA cohorts. The CIBERSORT analytical tool contains 547 genes that distinguish 22 human hematopoietic cell phenotypes, including seven T cell types, naive and memory B cells, macrophages of different phenotypes, plasma cells, natural killer cells and myeloid subsets [48]. We used this tool to estimate the proportions of the 22 TIICs in both ESCA samples and normal samples. To assess the accuracy of the deconvolution results, CIBERSORT derives a P-value for each sample. A CIBERSORT P-value < 0.05 was considered statistically significant and was used as the sample inclusion criterion. In total, 78 samples (2 normal samples and 76 ESCA samples in TCGA) were included. The CIBERSORT software package was obtained from the developers and combined with R (R.3.6.1). We used the default signature matrix with ESCA permutations for analysis.

Microarray data

Genes were screened using the following datasets from the GEO database (http://www.ncbi.nlm.nih.gov/geo): the GSE57614 series [49] and the GSE36537 series [50] on the GPL6480 platform (Agilent-014850 Whole Human Genome Microarray 4x44K G4112F), the GSE5099 series [51, 52] on the GPL96 platform (Affymetrix Human Genome U133A Array) and the GPL97 platform (Affymetrix Human Genome U133B Array), and the GSE95405 series on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array). According to the annotation information in each platform, probes were converted into corresponding gene symbols. The GSE57614 dataset contained nine samples each of M0 macrophages, M1 macrophages and M2 macrophages. GSE36537 contained nine samples each of M0 macrophages, M1 macrophages and M2 macrophages. GSE5099 contained six samples each of M0 macrophages, M1 macrophages and M2 macrophages. GSE95405 contained three samples each of M1 macrophages and M2 macrophages.

Identification of DEGs

GEO2R is an interactive web tool that allows users to compare two or more groups of samples in a GEO series in order to identify genes that are differentially expressed across experimental conditions. We used GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) to identify DEGs between M0 macrophages and M2 macrophages samples, and between M1 and M2 macrophages samples. The parameters used for the comparison of M0 macrophages and M2 macrophages samples were ‘‘adjusted P-value < 0.05, |logFC| ≥ 1.0’’. The parameters used for the comparison of M1 and M2 macrophages samples were ‘‘adjusted P-value < 0.01, |logFC| ≥ 1.0’’. The Venn online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to identify significantly DEGs, and overlapping DEGs in the Venn diagram were retained for further analysis.

Selection and analysis of hub genes in macrophage and M2 macrophage samples

The cBioPortal (https://www.cbioportal.org/) is a web resource for exploring, analyzing and visualizing multidimensional cancer genomics data [53]. We used cBioPortal to analyze the expression of hub genes in ESCA samples.

Kaplan-Meier Plotter (http://kmplot.com/analysis/) can be used to assess the effect of any gene or gene combination on survival in breast cancer, lung cancer, ovarian cancer, gastric cancer, colon cancer, prostate cancer, melanoma and 14 other tumor types. The database includes over 50,000 samples that were assessed using gene arrays and RNA sequencing [54]. We used Kaplan-Meier Plotter to determine the correlation of hub gene expression (FGL2, ERI1 and WNT5B) with the OS and RFS of ESCC patients, together with the HR and log-rank P-value.

The TIMER database is a comprehensive resource for the systematic analysis of immune infiltrates across diverse cancer types (https://cistrome.shinyapps.io/timer/) [55], and can be used to infer the abundance of TIICs from gene expression profiles according to a previously published deconvolution method [56]. Using the gene modules in TIMER, we assessed the correlation of hub gene expression with the abundance of immune infiltrates such as CD4+ T cells, CD8+ T cells, macrophages, B cells, DCs and neutrophils. We also analyzed the association of hub gene expression with ESCA purity [57].

Correlation analysis between the levels of FGL2 and TIICs

We used the correlation modules in TIMER to explore the correlations between the levels of FGL2 and genetic markers of TIICs (monocytes, TAMs, M1 macrophages, M2 macrophages, T cells [general], CD8+ T cells, B cells, Tregs, Th1 cells, Th2 cells, exhausted T cells, etc.). These genetic markers have been referenced in previous studies [58–60]. The correlation module generates an expression scatterplot between a pair of user-defined genes, as well as the Spearman’s correlation coefficient and statistical significance. Gene levels were calculated as log2 RSEM data.

The online database GEPIA (http://gepia.cancer-pku.cn/index.html) [61], a web server for cancer and normal gene expression profiling and interactive analyses, was used to confirm the correlations between the levels of FGL2 and markers of macrophages of different phenotypes. ESCA and normal tissue datasets were analyzed, and Spearman’s correlation coefficients were determined.

Correlation analysis between the levels of FGL2 and co-expressed genes and cytokines, and GSEA based on FGL2 expression

The publicly available LinkedOmics database (http://www.linkedomics.org/login.php) includes multi-omics data from all 32 TCGA cancer types, as well as mass spectrometry-based proteomics data generated by the Clinical Proteomics Tumor Analysis Consortium for TCGA breast, colorectal and ovarian tumors [62]. We used this database to identify DEGs in correlation with FGL2 in the ESCA cohort from TCGA (n=184). The Pearson method was used to determine the correlation coefficients. The LinkFinder also generates statistical expression plots for individual genes, including volcano plots and heatmaps. The correlations between the levels of FGL2 and cytokines in samples from TCGA were analyzed using cBioPortal. The Cancer Cell Line Encyclopedia (CCLE) have provided a rigorous framework with which to study genetic variants (https://portals.broadinstitute.org/ccle) [63]. FGL2 mRNA levels in esophageal cancer cell lines were determined using the CCLE.

GSEA can be used to determine whether a predefined gene set exhibits consistently significant differences in two biological states [64]. We used GSEA software (version: 4.0.0) to identify functions or pathways that differed between ESCA samples with low and high FGL2 mRNA levels. Differences were screened based on a normalized P-value < 0.05 and an FDR < 0.05, and 1000 gene set permutations were used in each analysis.

Macrophage and ESCC cell co-culture

Human ESCC (EC109/9706) cells and leukemia THP-1 monocytes were obtained from Fuxiang Biological Company (Shanghai, China). All cells were maintained in Roswell Park Memorial Institute 1640 medium (Gibco, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (Gibco), and were cultured in a humidified incubator containing 5% CO₂ at 37°C. Phorbol-12-myristate-13-acetate (Beijing 4A Biotech Co., Ltd) was added to the THP-1 monocytes to a final concentration of 40 ng/mL. Approximately 36 h later, the cell morphology was observed under an inverted microscope to confirm that the THP-1 monocytes had been induced to differentiate into M0 macrophages. Then, IL-4 and IL-13 (PeproTech, USA) were added to the M0 macrophages to final concentrations of 8 ng/mL and 4 ng/mL, respectively. After 48 h, the cell morphology was observed under an inverted microscope to confirm that the macrophages had been induced to differentiate into M2 macrophages. The EC109/9706 cells and M2 macrophages were respectively cultured in the upper and lower chambers of six-well plates (Corning Inc., Corning, NY, USA), each with 2 mL of culture medium. The cells were placed in an incubator for 48 h to establish a non-contact co-culture. The non-contact co-culture of EC109/9706 cells and M0 macrophages was established in the same way.

Human cytokine antibody array

We used the Human XL Cytokine Array Kit (R & D; Cat ARY022B) to detect the expression of various cytokines in cell co-culture samples. For each Array unit, 200-500 μg of cell culture supernatant was used. The protein chip was blocked with 2 mL of Array Buffer and incubated for 1 h on a shaker. The samples were diluted with Array Buffer to a volume of 1.5 mL per membrane and then incubated with the Array protein membrane at 4 °C overnight. The membrane was washed three times with wash buffer, treated with a diluted detection antibody and incubated at room temperature for 1 h. The membrane was then washed, treated with streptavidin-horseradish peroxidase, incubated for 30 min and shaken at room temperature. After the membrane was washed, 1 mL of the chemical color reagent was added to each well, and the color was developed on a chemiluminescence imager. Professional software was used to process and analyze the data.

Analysis of DEGs between M1 and M2 macrophages samples

WebGestalt (http://www.webgestalt.org) [65], one of the most widely used GSEA tools, was used to analyze the KEGG pathways and GO terms (Biological Process, Molecular Function and Cellular Component) of the DEGs between M1 and M2 macrophages. The rank criterion was an FDR < 0.05. The most important module in the PPI network was built using the STRING online database (http://string-db.org, Version: 11.0) [66]. A comprehensive Gt score > 0.7 was considered statistically significant. Cytoscape [67] was used to cluster the resulting network to reveal closely connected regions.

Data availability

All datasets generated for this study are included in the manuscript.

Ethics statement

As the data (TCGA and GEO datasets, etc.) are publicly available, no ethical approval was required.