Development and Validation of a Prognostic Gene Signature Correlated With M2 Macrophage Infiltration in Esophageal Squamous Cell Carcinoma

The level 3 gene expression profile and corresponding clinical information of ESCC patients were downloaded from UCSC Xena (dataset ID: TCGA-ESCA-sampleMap/HiSeqV2, https://xenabrowser.net/datapages/↗). The gene expression profile was measured experimentally using the Illumina HiSeq 2000 RNA Sequencing platform by the University of North Carolina TCGA genome characterization center. Gene expression was provided as gene-level transcription estimates with units as log2(x+1) transformed RSEM normalized count. Low-expression genes with mean expression values below 1 RSEM in all samples were filtered out using the “limma” package in R version 4.0.2 software. The original data included 185 tumor tissues and 11 adjacent tissues. Among these, 96 samples were histologically diagnosed as squamous cell carcinoma, including 1 paired metastasis tissue. The clinical information of the patients is shown in Supplemental Table 1. The risk score model was further validated using the GSE53624↗ dataset from the Gene Expression Omnibus (GEO, http://www.ncbi↗. nlm.nih.gov/geo/). The screening process of the validation dataset is provided in Supplemental Figure 1. The GSE53624↗ dataset included 119 paired tumor and normal ESCC tissues based on GPL18109↗ platform (Agilent-038314 CBC Homo sapiens lncRNA + mRNA microarray V2.0).

Immune, stromal, and ESTIMATE scores of the samples were calculated using the estimate R package. We determined the optimal cutpoint based on the function “surv_cutpoint” from the survminer R package. Kaplan–Meier analyses were performed using the survival and survminer packages in R to illustrate the correlation of immune/stromal scores and patient overall survival (OS). The log-rank test was applied to verify the results.

According the optimal cutpoint, immune and stromal scores were divided into high/low groups, respectively. DEGs were identified using the R package, limma. The threshold set for up- and down-regulated genes was a |log 2 foldchange (FC)| > 1 and false discovery rate (FDR) < 0.05. Heatmaps were plotted with the package pheatmap.

Functional enrichment analyses of the DEGs with the Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene ontology (GO) were performed with the R package “clusterProfiler”, “org.Hs.eg.db”, “enrichplot” and “ggplot2”. GO enrichment includes biological processes (BP), cellular component (CC) and molecular function (MF). Categories with a p- and q-value of <0.05 were considered significantly enriched. All of the DEGs were uploaded into the STRING (https://string-db.org/↗) database (v 11.0) to obtain PPI networks, with a combined score > 0.4 considered statistically significant. Cytoscape (version 3.7.1) was used to reconstruct the network. Network nodes represent proteins and edges represent protein-protein associations.

Univariate cox regression analysis was performed to examine the prognostic value among ESCC patients. 36 genes with p < 0.05 were identified as prognostic DEGs and were visualized using the forest diagram. Multivariate analyses were performed to develop the 10-gene risk score model. The model was based on expression data multiplied by Cox regression coefficients. The final risk score model formula was as followed: Risk score = [Expression level of C1QA * (0.32596)] +[Expression level of C1QB * (1.40234)]+[Expression level of C1QC* (-1.36687)+ [Expression level of CD86 * (0.35249)]+[Expression level of C3AR1* (-0.07155)+ [Expression level of CSF1R * (0.28413)]+[Expression level of ITGB2 * (-0.52918)]+[Expression level of LCP2* (-0.22934)+ [Expression level of SPI1 * (1.05328)]+[Expression level of TYROBP* (-0.87910)]. Patients were divided into low-risk and high-risk groups according to the optimal cutpoint. The K-M survival curves for the groups with low or high risk were performed. The predictive ability of the model was assessed by the survival receiver operating characteristic (ROC) package in R software and was used to compare the area under the curve (AUC) of our gene signature and those derived in two other published studies. To confirm the risk model’s independent prognostic value, univariate and multivariate Cox survival analyses were performed with select clinical factors. Finally, external data from GSE53624↗ was applied to verify the reliability of the gene signature’s impact on the prognosis of the patients. The differential expression analysis was performed based on 119 paired tumor and normal samples.

CIBERSORT in combination with the LM22 method was carried out to quantify the abundances of immune cell types in the TME. The 22 types of infiltrating immune cells inferred by CIBERSORT include B cells, T cells, natural killer cells, macrophages, dendritic cells, eosinophils and neutrophils. The CIBERSORT p-value reflects the statistical significance of the results, only tumor samples with p<0.05 were used for further analysis.

A nomogram was established based on the risk score and select clinicopathologic characteristics including age, gender and stage to predict the survival probability of 1-, 2-, and 3-year OS of ESCC patients. The nomogram and calibration plots were generated based on the rms R package. The calibration curve of the nomogram was plotted to evaluate the prediction possibilities against the observed rates.

GSEA was carried out to evaluate associations between immune pathways and 10-gene signature using the software GSEA-4.1.0. Hallmark (h.all.v7.4.symbols.gmt) and C7 gene sets (c7.all.v7.4.symbols.gmt, Immunologic Signatures) were downloaded from Molecular Signatures Database (MSigDB) as the target sets. Only gene sets with NOM p < 0.05 and FDR q < 0.25 were considered significant.

Tissue cDNA chips including cDNA from 67 cases of esophageal squamous cell carcinoma tissue and 28 peri-carcinoma tissues with complete clinical and survival information were purchased. A cDNA microarray chip (cDNA-HEsoS095Su01,Outdo Biotech Company, Beijing, China) was used for the tumor or peritumor tissue samples in this study. The mRNA expression levels of hub genes and immune cell surface biomarkers were detected by Hieff qPCR SYBR Green Master Mix (Low Rox Plus) (YEASEN Biotech Co., Ltd). The qPCR protocol was 95°C for 5min, 40 cycles at 95°C for 10 s, and 60°C for 60 s. Primers used in this study are presented in Supplemental Table 2. The relative expression levels of hub genes were determined by the 2^−ΔΔCT method.

The correlation analysis was performed using the Spearman method. Survival curves were compared using the Kaplan-Meier method and the log-rank test. R version 4.0.2 and GraphPad 5.0 were used to perform statistical analysis. The results are presented as mean ± standard deviation (mean ± SD). Differences between groups were evaluated by the Wilcoxon rank-sum test. All tests were two-sided and p < 0.05 indicates statistical significance.

The analysis process for this study is presented in Figure 1. We estimated the immune/stromal/ESTIMATE scores of 95 ESCC tumors using the ESTIMATE algorithm. The immune scores varied from -1389.05 to 3395.19, the stromal scores ranged from -1859.96 to 1301.57, the ESTIMATE scores ranged from -2672.1 to 3970.4. All samples were categorized to high/low groups with the optimal cutpoint. Using Kaplan-Meier analysis, we found that patients with higher immune scores experienced poorer overall survival (OS) compared to those with low scores (p=0.015, Figure 2A). In agreement, stromal scores were inversely correlated with OS (p=0.012, Figure 2B). The high ESTIMATE score group also showed poorer OS in comparison to the low score group (p=0.057, Figure 2C). These findings demonstrate that the immune/stromal components in TME are significant in predicting the prognosis of ESCC patients. We further explored the potential relationship between the clinicopathological characteristics and the immune/stromal scores, but found no significant associations.

To investigate the potential relationship between gene expression profiles and immune/stromal scores, we performed differential expression analysis of high- and low-score groups. According to the analysis (high immune score vs low immune score), a total of 747 DEGs were selected, which contained 674 up-regulated and 73 down-regulated DEGs. Similarly, 1027 DEGs were obtained based on differential analysis of stromal scores, consisting of 855 up-regulated and 172 down-regulated DEGs. Hierarchical clustering (Figures 3A, B) showed that DEGs were significantly dysregulated between the two groups. Furthermore, after intersection of the two lists of genes (Figures 3C, D), we obtained 133 up-regulated DEGs shared by immune and stromal groups. These DEGs can be regarded as candidate TME-related genes.

Furthermore, functional enrichment analysis was performed to find the potential mechanism of 133 DEGs in the TME. The top 10 significant results of the enrichment analysis for BP, CC and MF are displayed in Figure 3E. For BP, DEGs were mainly enriched in immune effector processes, negative regulation of immune system processes, and positive regulation of cell activation. In the CC group, DEGs were mainly enriched in the external side of plasma membrane, protein complex involved in cell adhesion, and secretory granule membrane. In the MF classification, the top terms included immunoglobulin binding, immune receptor activity, chemokine activity, cytokine receptor activity, and chemokine receptor binding. For the KEGG analysis (Figure 3F), DEGs were mainly enriched in the Natural killer cell mediated cytotoxicity, Toll−like receptor signaling pathway, viral protein interaction with cytokine and cytokine receptor, cell adhesion molecules, Th17 cell differentiation, and chemokine signaling pathways. Collectively, these findings suggest that the functions of these genes are immune-related.

To better illustrate the interrelationship among these DEGs, we used the STRING database and the Cytoscape software to construct a PPI network. As depicted in Figure 4A, the color from light to dark represents the ascending logFC value. The barplot (Figure 4B) shows the top 20 hub genes ranked by the number of nodes. Next, we carried out univariate cox regression analysis and identified 36 genes significantly associated with the poor OS of ESCC patients (HR>1, p<0.05). Finally, the 10 independent prognostic genes in the PPI and univariate cox regression analysis were overlapped to identify the common hub genes, including C1QA, C1QB, C1QC, CD86, C3AR1, CSF1R, ITGB2, LCP2, SPI1 and TYROBP (Figure 4C). The forest diagram (Figure 4D) illustrates the relationships between these 10 genes and prognosis. The Kaplan-Meier curve confirmed that high expression of 9 of the genes (C1QA, C1QB, CSF1R, C3AR1, ITGB2, LCP2, SPI1, TYROBP) were significantly related to poor OS (Supplemental Figures 2A–I, p<0.05). Moreover, for those variables in which survival curves intersected, we carried out landmark analysis discriminating between events occurring before and after 3 years of follow-up to eliminate immortal time bias, which confirmed that all genes are significantly related to overall survival before 3 years (Supplemental Figure 3B, p-value<0.05).

The CIBERSORT algorithm was then used to quantify the proportion of tumor-infiltrating immune subsets and further understand the correlation of hub genes with the immune TME. Figure 5A shows the relative proportions of immune cells in each ESCC sample. Positive and negative correlations between immune cells were obtained, as displayed in Figure 5B. There was a moderate correlation between macrophage M0 and macrophage M2 (r=-0.49) in ESCC. CD8 T cells were moderately and positively correlated with activated CD4 memory T cells (r=0.27), while negatively correlated with resting CD4 memory T cells (r=-0.62). M1 macrophages were positively associated with Tregs (r=0.23) and M2 macrophages (r=0.24), and negatively associated with activated dendritic cells (r=-0.44). These results reveal that different kinds of immune cells interfere with each other in TME. Pearson correlation analysis was further performed to compare the expression of ten genes and the results of ESTIMATE. All ten genes were strongly positively correlated with ImmuneScore but negatively correlated with TumorPurity (Figure 5C, p <0.05), demonstrating lower the purity of the tumor and increased immune cells in the TME. Moreover, we searched the CCLE (Cancer Cell Line Encyclopedia) database to help clarify whether these genes are overexpressed in immune cells or tumor cells. The expression values for the ten genes and classic epithelial markers in esophageal cancer cell lines was downloaded. All ten genes showed extremely low expression in esophageal cancer cell lines (Supplemental Figure 4), with average expression values (RPKM) of: C1QA 0.0005, C1QB 0.0002, C1QC 0, C3AR1 0.0385, CD86 0.2466, CSF1R 0.1609, ITGB2 1.0324, LCP2 0.0123, SPI1 0.0255, TYROBP 0.0086. Comparably, the epithelial markers selected as controlled genes showed relatively high expression CDH1 57.414, CLDN4 60.703, CLDN7 50.385, MUC1 13.366, TJP3 4.5215. Therefore, these combined evidences lead us to believe that these genes are likely to be specifically expressed in immune cells. To investigate this further, we performed a single-gene immune infiltration analysis for each gene to illustrate its relationships with various immune cells. Most of the genes were related to certain types of immune cells. As depicted in Figure 5D, three kinds of tumor infiltrating immune cells (TICs) were positively correlated with all genes, including M1 macrophages, M2 macrophages, and regulatory T cells (Tregs). M2 macrophages had a higher correlation with C1QA (r=0.54), C1QB (r=0.55), C1QC (r=0.53), C3AR1 (r=0.48), CSF1R (r=0.46) and TYROBP (r=0.44). M1 macrophages were moderately linked with LCP2 (r=0.46), ITGB2 (r=0.43), CSF1R (r=0.44), C1QA (r=0.43), C1QB (r=0.44), C1QC (r=0.42) and C3AR1 (r=0.41). All genes were positively associated with M2 macrophages (Supplemental Figures 5 A–J, p<0.05). These findings further confirm that these prognostic genes are related to the immune activity of the TME, especially M2 macrophages.

Utilizing the multivariate Cox regression analysis, we established a 10-gene risk score model. The risk model met the proportional hazard assumption based on the Schoenfeld Individual Test results which showed that each covariate is not statistically significant (Supplemental Figure 2A, p>0.05). We then calculated the risk score for each sample and divided patients into high- or low-risk group according to the optimal cutpoint by maximally selected rank statistics. Kaplan-Meier analysis suggested that patients in the high-risk group have significantly shorter OS than those in the low-risk group (Figure 6A, p<0.001). The distribution of risk scores, survival status, and ten-gene expression levels among patients in the high- and low-risk group are given in Figure 6C. To evaluate the independent predictive value of our risk model, univariate and multivariate analyses were performed. The results showed gender, stage, N-stage and risk score were significantly associated with OS in univariate analysis (p<0.05). In the multivariate analysis, only risk score was associated with OS (HR=2.104, 95% CI=1.343-3.295; p=0.001) (Table 1). To verify the prognostic value and reliability of our results, the risk score model was further validated using the GEO dataset, which includes 119 ESCC patients. All patients were divided into high- and low-risk group according to the previous formula. In agreement with the training cohort, patients in the high-risk group had significantly worse OS than the low-risk group (Figure 6B, p-value=0.008). The univariate and multivariate analyses of risk score and other clinical characteristics confirmed that the risk score model was an independent prognostic indicator (HR=1.6915, 95% CI:1.053-2.717; p=0.0297, Table 2). We also carried out differential expression analysis of the complete ten-gene signature in the validation dataset, The Wilcoxon rank sum test revealed that the expression of 9 genes (C1QA, C1QB, C1QC, C3AR1, ITGB2, LCP2, SPI1, TYROBP) in the tumor samples were significantly higher than that in matched normal adjacent tissue (Supplemental Figures 6A–J, p<0.05). Time-dependent ROC curve analysis demonstrated that during 1-, 2-, and 3-year follow-up, the area under the curve (AUC) values were 0.672, 0.854, and 0.81 respectively (Figures 6E–G) Recently, Sun et al. identified a prognostic gene signature among patients with ESCC (23). Zhang et al. constructed a prognostic model based on immune-related genes to predict prognosis of esophageal cancer (24). We calculated C-indexes to compare the prognostic values of our model and theirs. As shown in Figures 6E–G, the concordance index of the risk score model for the 1-, 2-, and 3-year OS was higher than the other two studies, indicating that our risk score model may have better performance in predicting prognosis. Taken together, the results confirm that the risk score model is an reliable and independent prognostic factor for ESCC patients.

We estimated the difference of immune infiltration between high- and low-risk ESCC patients in 22 subpopulations of immune cells using the CIBERSORT algorithm. Levels of M1 macrophages (p <0.001) and M2 macrophages (p < 0.001) were markedly higher in the high-risk compared to the low-risk group. Additionally, the proportion of M2 macrophages were significantly higher than M1 macrophages (Figure 7A). In contrast, a high fraction of resting CD4 memory T cells, activated dendritic cells, and M0 macrophages mainly infiltrated low-risk ESCC patients. The correlation analysis further confirmed the risk score was moderately correlated with M1 macrophages (r = 0.43, p<0.0001) and M2 macrophages (r = 0.47, p <0.0001). Moreover, there were weak associations between risk score and CD8 T cells (r=0.21, p=0.046), and a negative correlation between risk score and M0 macrophages (r=-0.36, p<0.001, Figure 7B). These results indicate that differences in immune infiltration in high- and low-risk patients with ESCC might be used as a prognostic indicator and target for immunotherapy.

To establish a more convenient and applicable clinical prognostic approach, we developed a nomogram based on our risk score and other clinical characteristics including age, gender and pathologic stage (Figure 8A). The concordance index of the nomogram was 0.734. The calibration plot for the possibility of 1-, 2- and 3-year survival showed good agreement between the prediction and actual observations (Figure 8B). These findings illustrate that the nomogram may be a more effective method to predict prognosis of ESCC patients for clinicians.

To help illustrate the underlying mechanism of how these genes impact patient outcomes, GSEA was utilized to evaluate different expression profiles among the two groups. Results showed that genes in the high-risk group were mainly enriched in several immune-related pathways, such as allograft rejection, IL2_STAT5_SIGNALING, IL6_JAK_STAT6_ SIGNALING, KRAS_ SIGNALING_UP and PI3K_AKT_MTOR_SIGNALING (Figure 9A). However, no significant gene sets were enriched in the low-risk group. For the C7 immunologic gene sets defined by MSigDB, multiple immune functional gene sets were enriched in high-risk group (Figure 9B) while only three gene sets were enriched in the low-risk group (Figure 9C). These results suggest that these genes might affect tumor immune status through these pathways.

Finally, we validated five rarely reported genes among these ten genes using an ESCC cDNA microarray by qRT-PCR. As demonstrated in Figures 10A–E, the relative expression level of C1QA, C3AR1, LCP2 and TYROBP in ESCC samples were significantly higher relative to normal samples (p < 0.05). The expression of SPI1 was also up-regulated although not statistically significant. The survival curve showed that they were all significantly associated with the overall survival of ESCC (log-rank p < 0.05; Figures 10F–J). In addition, based on the CIBERSORT results, we further examined the correlation between 5 genes and macrophage M1 and M2 surface markers. As depicted in Figure 11, compared with CD8 and CD86, they had stronger correlation with CD206 and CD4, which further verified they may be involved in the activity of M2 macrophages and play an immunosuppressive role in TME.