Macrophage Infiltration and ITGB2 Expression in ESCC : A Novel Correlation

The gene chip data were obtained from four series in the GEO (Gene Expression Omnibus) database: GSE161533↗, GSE23400↗, GSE66274↗, and GSE67268↗. To enhance data reliability, all chip data underwent cleaning, standardization, inspection, and batch effect removal. Subsequently, probe annotation and gene ID conversion were performed.

The high‐throughput sequencing data were sourced from The Cancer Genome Atlas (TCGA↗) and the Genotype‐Tissue Expression (GTEx↗) Database. To ensure data reliability, we utilized the TCGA original count data from UCSC Xena↗ and systematically checked each TCGA sample using the “Merged Sample Quality Annotations” file from the TCGA website, removing samples with poor quality. Only high‐quality samples with stage “A” were retained. Moreover, for ESCC samples, we selected samples with “Sample Type Codes” of “01”, “02”, “05”, or “08”. For normal tissue samples, we only chose samples with “Sample Type Codes” of “11”. Considering the different sources of data from TCGA and GTEx, direct merging analysis could introduce significant bias and affect the accuracy of results. Therefore, we utilized the “UCSC Toil RNAseq Recompute Compendium” data [11] from UCSC Xena for merged analysis to minimize bias arising from different data sources. For the original count data, batch effect correction was performed using the “batch_number” provided by UCSC Xena. Since the data from “UCSC Toil RNA‐seq Recompute Compendium” did not provide batch effect reference, we explored and removed batch effects using the R package “RUVSeq(v1.32.0)”.

The single‐cell sequencing data were obtained from the GEO database under accession number GSE199619↗. We utilized the “GSE199619_ELN.integrated.rds.gz” file. To minimize bias, we selected ESCC samples that underwent surgical resection or endoscopic resection without neoadjuvant chemotherapy (NACT) (Table S2).

Summary information for all data is provided in Table. S1

We conducted immunoinfiltration analysis using ImmuCellAI [12, 13]. The preprocessed gene expression matrices from GSE161533↗ and GSE23400↗ were imported into ImmuCellAI, and the “Analysis” function was utilized to quantify the infiltration of various immune cells in tumor samples.

Subsequently, based on the gene expression matrices and corresponding immunoinfiltration analysis results, we constructed a weighted gene co‐expression network using the R package “WGCNA (1.71)” and screened for gene modules associated with various immune cell infiltrations. Specifically, we first selected 10,000 genes with large standard deviations and identified outlier samples through sample clustering, followed by their removal. Then, we chose an appropriate soft threshold to construct the scale‐free network and TOM matrix. Next, we calculated the correlation between gene modules and immune cell infiltration using the “moduleTrait” function. If multiple associations were found between immune cells and gene modules, we applied two criteria for selection: (1) selecting immune cells with higher average infiltration abundance and (2) conducting correlation tests with p < 0.01. Subsequently, we obtained immune infiltration‐related genes by intersecting the gene modules obtained from GSE161533↗ and those obtained from GSE23400↗.

We annotated the biological significance of intersecting genes using the R package “clusterProfiler (4.6.0)” to indirectly demonstrate the correlation between these genes and infiltrating immune cells. This annotation involved two enrichment analysis methods: Gene Ontology (GO) and Kyoto Encyclopedia of Genes (KEGG).

Subsequently, we imported the intersecting genes into the STRING database to construct a protein–protein interaction (PPI) network. Associations with the “combined score” less than 0.7 were excluded to obtain a reliable PPI network. The network was then imported into Cytoscape (3.9.2), and key genes in the network were calculated using the “MCC” method in the cytoHubba tool. Finally, the top ten ranked genes were selected as the core genes of the PPI network.

For chip data, we used the R package “limma (3.54.2)” to perform paired comparisons between ESCC samples and normal esophageal tissue samples from the GSE161533↗ and GSE23400↗ datasets to analyze gene expression differences. For high‐throughput sequencing data, we extracted RNA counts data for ESCC and normal esophageal tissue from the “UCSC Toil RNAseq Recompute Compendium” dataset and conducted differential analysis using the R package “DESeq2 (1.38.1)”.

As for survival analysis, we obtained data from TCGA. Utilizing the R package “survminer (0.4.9)”, we determined the optimal cutoff point and then assessed whether there were differences in overall survival (OS), disease‐specific survival (DSS), and progression‐free interval (PFI) between patients with high expression of core genes and those with low expression using the R package “survival (3.4.0)”.

Finally, we selected the core genes from the PPI network that were significant in both all expression differential analyses and all survival analyses as target genes. The immune cells corresponding to the gene module where the target genes were located were identified as target cells.

First, we compared the expression differences of the target genes in different cells using the R package “Seurat(4.3.0)”. Subsequently, we explored whether the expression of the target genes was related to the differentiation and development of target cells using the R packages “IOBR(0.99.9)” and “Monocle(2.26.0)”. We also analyzed whether different tumor stages affect the expression of target genes in the target cells.

We conducted western blot (WB) and immunohistochemistry (IHC) experiments to examine the protein expression of the target genes. The samples used in the experiments were obtained from ESCC patients, with detailed information provided in Table. All patients were unrelated Asians who were hospitalized and treated at the First Affiliated Hospital of Guangxi Medical University (Guangxi Province, China). ESCC was confirmed by histopathological examination of tissue obtained from surgical resection of the tumor or biopsy. None of the patients received chemotherapy or radiotherapy before tumor resection. Biological samples were collected immediately after tumor resection from patients and analyzed according to the specified procedures. This study was approved by the Ethics Committee of the First Affiliated Hospital of Guangxi Medical University, and informed consent was obtained from each patient. The experimental procedures are briefly described as follows. S3

For the WB, β‐actin (Affinity, AF7018, 1:9000) was used as the loading control, and ITGB2 (Abclonal, RA2173, 1:1000) as the target protein. SDS‐PAGE gels were first prepared, followed by protein extraction from tissue samples, sample denaturation, and electrophoresis. After electrophoresis, membrane transfer was performed using the wet transfer method. Finally, the membranes were scanned using a scanner, and the relative protein expression levels were assessed with ImageJ software.

For the IHC, ITGB2 (Abclonal, RA2173, 1:100) was used as the primary antibody. Formalin‐fixed, paraffin‐embedded tissue sections were first deparaffinized, followed by antigen retrieval. After removing endogenous peroxidase activity, the sections were blocked. Then, primary antibody incubation, secondary antibody incubation, and DAB staining were performed. Subsequently, counterstaining (nuclear staining), mounting, and imaging were carried out. Finally, ImageJ software was used to quantify the relative expression levels of ITGB2.

For chip data, we extracted the infiltration abundance of target cells from the previous ImmuCellAI results, and then performed correlation analysis between the expression values of target genes extracted from the gene expression matrix and the infiltration abundance of target cells. For high‐throughput sequencing data, we first downloaded the immune infiltration data for esophageal cancer from the ImmuCellAI website, extracted the infiltration abundance of target cells in ESCC, and then converted the count data of ESCC to TPM data. Subsequently, we conducted correlation analysis between the TPM data of target genes and the infiltration abundance of target cells.

To enhance the reliability of the results, we additionally utilized six different algorithms from the TIMER2↗ website (TIMER, CIBERSORT, CIBERSORT.ABS, QUANTISEQ, XCELL, EPIC) to calculate the immune infiltration abundance of target cells in ESCC samples derived from TCGA. Subsequently, we performed correlation analysis between the results of these additional algorithms and the TPM data of target genes.

Subsequently, we utilized gene expression data from the GSE16533↗, GSE23400↗, and TCGA datasets to perform correlation analysis and assess the relationship between the target genes and the markers of target cells.

Finally, we performed dual immunofluorescence (IF) to simultaneously detect the expression of the target genes and the target cell markers in ESCC, providing further evidence for the association between the target gene and the target cells. Co‐localization analysis of the target genes expression and the target cell marker expression was conducted using the Coloc 2 plugin in ImageJ software. The experimental procedures are briefly described as follows.

The tissue section samples were obtained from ESCC patients, with detailed information provided in Table. The sources of the patients have been described previously. ITGB2 (Abcam, 10,554‐1‐AP, 1:50) and CD163 (Abcam, 68,218‐1‐lg, 1:50) were used as primary antibodies (target proteins). First, the tissue sections were deparaffinized, followed by antigen retrieval in either citrate buffer (pH 6.0) or Tris‐EDTA buffer (pH 9.0). Blocking was then performed using a blocking agent/antibody dilution solution. The sections were subsequently incubated with the primary antibodies. Immunostaining was visualized using tyramide signal amplification (TSA) with fluorophore‐conjugated secondary antibodies specific to each primary antibody. After each round of antibody staining, a microwave‐based stripping procedure was employed. Following staining, the sections were scanned and imaged using a microscope. Finally, ImageJ software was used to quantify the relative expression levels of ITGB2 and CD163. Co‐localization analysis of ITGB2 and CD163 expression was performed using the Coloc 2 plugin in ImageJ software. S3

We identified the primary known targets of the target cells through a review of relevant literature. Using gene expression data from GSE161533↗, GSE23400↗, and TPM data from TCGA, we conducted correlation analyses to determine the associations between these targets and the target genes.

We selected TME signature gene sets related to immunotherapy response based on relevant literature and gene sets collected by R package “IOBR (0.99.9)”. Utilizing the ESCC data from the GSE23400↗ dataset and the TPM data from TCGA, we quantified these TME signatures using the “ssGSEA” algorithm in the R package “IOBR(0.99.9)”. Subsequently, correlation analysis was conducted to understand the association between the target genes and these signatures.

Initially, we utilized miRWork↗ to predict miRNAs that may potentially bind to the target genes. Subsequently, we filtered out miRNAs with a binding probability greater than 0.9. Differential expression analysis was conducted using data from GSE66274↗ and GSE67268↗ to identify miRNAs exhibiting differential expression between ESCC and esophageal normal tissues. Following this, we further refined the selection by filtering for miRNAs with an adjusted p < 0.01 and a fold change (log₂FC) less than −2. Finally, we intersected the miRNAs obtained from GSE66274↗, those obtained from GSE67268↗, and those selected from miRWork to obtain the final miRNAs.

The immune cell profiling using ImmuneCellAI revealed a significant presence of dendritic cells, B cells, macrophages, NK cells, and other immune cell types within ESCC (Figure 1A). However, the abundance of other immune cell types is comparatively lower. Consequently, these ten immune cell types were included for subsequent analysis.

For GSE161533↗, we set the “cutHeight” parameter of WGCNA to 95, removing 3 outlier samples (GSM4909616↗, GSM4909627↗, GSM4909622↗), and then selected a soft threshold of 6 to fit the optimal scale. At this point, R² = 0.901 (Figure S1A), and the average network connectivity was 33.5 (Figure S1C). After constructing the weighted gene co‐expression network and merging modules with correlation coefficients greater than 0.75, all genes were assigned to 16 different colored gene modules (Figure S1E). Through correlation analysis, we found multiple associations between gene modules and immune cell infiltration (Figure S2A), such as a significant correlation between the magenta module and macrophage infiltration (p = 0.001), and between the yellow module and NK cell infiltration (p = 0.001).

For GSE23400↗, we set the cutHeight to 65, removing 4 outlier samples (GSM573938↗, GSM573951↗, GSM573901↗, GSM573906↗), and then selected a soft threshold of 5 to fit the optimal scale. At this point, R² = 0.851 (Figure S1B), and the average network connectivity was 26.9 (Figure S1D). After constructing the weighted gene co‐expression network and merging modules with correlation coefficients greater than 0.75, all genes were assigned to 12 different colored gene modules (Figure S1F). Through correlation analysis, we also found multiple associations between gene modules and immune cell infiltration (Figure S2B), such as a significant correlation between the yellow module and neutrophil infiltration (p = 0.001), and between the red module and macrophage infiltration (p = 0.002).

Therefore, we screened using the two predefined filtering criteria and determined one gene module from each of the two sets of WGCNA results for further analysis. The magenta module was identified for GSE161533↗, and the red module for GSE23400↗ (Figure S2). The corresponding immune cell was macrophages. Subsequent analysis revealed a significant positive correlation between module membership and gene significance for both groups (Figure 1B,C). Thus, taking the intersection of the two gene modules yielded 108 genes associated with macrophage infiltration (Figure 1D).

The outcomes of KEGG analysis suggested that the intersecting genes were linked to biological processes like “Phagosome,” “Osteoclast differentiation,” “Lysosome,” “NF‐kappa B signaling pathway,” and “Antigen processing and presentation” (Figure). GO analysis results indicated that the intersecting genes were associated with biological characteristics including “macrophage activation,” “phagocytosis,” “cell activation involved in immune response,” “myeloid leukocyte activation,” “myeloid leukocyte migration,” “ficolin‐1‐rich granule,” “ficolin‐1‐rich granule lumen,” “phagocytic cup,” “endocytic vesicle membrane,” “integrin complex,” “mannose binding,” “immunoglobulin binding,” “immune receptor activity,” “pattern recognition receptor activity,” and “complement binding” (Figure). These findings further emphasize the probable correlation between the expression of intersecting genes and macrophage infiltration in ESCC. S3A S3B

As shown in Figure, using the STRING database, we obtained the protein–protein interaction (PPI) network of the relevant genes. Associations with a “combined score” less than 0.7 were excluded (Figure), and using the “MCC” method of the cytoHubba tool for calculation, we ultimately identified ten core genes in the network (Figure). Ranked from high to low based on their scores, they are: CSF1R, ITGB2, ITGAM, TYROBP, FCGR2A, CD14, FCER1G, ITGAX, TLR4, and C1QB. S4A S4B S4C

Gene expression differential analysis and survival analysis revealed significance in ITGB2, ranked second among these ten core genes, in all analyses (Figure 2, Tables S4,S5). Results from three different datasets consistently showed higher expression of ITGB2 in ESCC compared to esophageal normal tissues (all p < 0.05, Figure 2A–C). Patients with high expression of ITGB2 in ESCC had shorter OS, DSS, PFI times compared to those with low expression of ITGB2 (all p < 0.05, Figure 2D–F).

Furthermore, both IHC and WB results demonstrated significant upregulation of ITGB2 protein in ESCC compared to esophageal normal tissues (all p < 0.05, Figure 2G–J).

These findings suggest that ITGB2 is a promising and worthy gene for further investigation. The overexpression of ITGB2 is likely associated with the occurrence and progression of ESCC, and this correlation may be closely related to macrophage infiltration. Therefore, we identified ITGB2 as the target gene for this study. Correspondingly, macrophages were identified as our target cells.

The single‐cell sequencing samples comprised 20,793 cells. We classified cells using DCN, PDPN, CD2, CD3D, ITGAX, CD68, CD79A, KRT5 EPCAM, and TPSB2 as markers (Figure S5C). Subsequently, we isolated macrophages from myeloid cells using CD14 as a marker (Figure S5D), obtaining a total of 3491 macrophages (Figure S5A,B). We observed that ITGB2 expression in macrophages was higher than in other cell types (Figure 3A, Table S6).

We then performed pseudotime analysis on the macrophages. Following the widely accepted cancer immune editing theory [14], we assessed whether pseudotime conformed to an objective pattern; namely, macrophages in the TME would progressively acquire pro‐tumor characteristics. Over time, macrophages would exhibit increasingly significant M2 features. We quantified the strength of M2 features in each macrophage using the “ssGSEA” algorithm in the R package “IOBR.” Details of the M2 feature gene set are provided in Table S7.

The results confirm the reliability of our pseudotime analysis. As time progresses, the M2 characteristics of macrophages become more prominent (Figure 3B,C). We found a positive correlation between the expression of ITGB2 within macrophages and pseudotime, as well as M2 characteristics (all p < 0.05, Figure 3E,F). Further analysis revealed that in the microenvironment of advanced ESCC, ITGB2 expression within macrophages is higher compared to that in the microenvironment of early ESCC (adjusted p = 0.54e‐43, Figure 3D).

In all three different datasets, the expression of ITGB2 was significantly positively correlated with macrophage infiltration (p < 0.05, Figure 4A–C). In the GSE161533↗ dataset, the Spearman correlation coefficient (rho) was 0.7 with a p‐value of 4e‐05 (Figure 4D); in the GSE23400↗ dataset, the Spearman correlation coefficient was 0.39 with a p‐value of 3.7e‐03 (Figure 4E); in the TCGA dataset, the Spearman correlation coefficient was 0.64 with a p‐value of 1.9e−10 (Figure 4F). Results from multiple additional methods also indicated a significant positive correlation between ITGB2 expression and ESCC macrophage infiltration (Figure S6).

Subsequently, we analyzed the correlation between ITGB2 and M2 macrophage markers MRC1 (CD206) and CD163. We found that ITGB2 expression was significantly positively correlated with MRC1 and CD163 across different datasets (Table). S8

Finally, we performed dual IF to simultaneously detect the expression of ITGB2 and CD163 proteins in ESCC tissues (Figure 4G). Correlation analysis revealed that the results were consistent with those from the three datasets mentioned above, showing a significant positive correlation between ITGB2 and CD163 expression (Figure 4H). Using the Coloc 2 plugin in ImageJ software, we further identified a co‐localization relationship between ITGB2 and CD163 expression (Figure 4I,J). These findings further support the association between ITGB2 and macrophage infiltration in ESCC.

Through a review of relevant literature, we identified ten major known targets of macrophages, namely CCL2, CCR5, CSF1R, CD47, CD40, TLR3, TLR7, and TREM2. Clinical trials related to these targets are currently underway, with some having achieved certain results (Table S11). In the GSE161533↗ dataset, we found that ITGB2 expression was positively correlated with these targets (all correlation coefficients greater than 0, all p < 0.05), except for TLR3 (p > 0.05) (Figure 5A). In the GSE23400↗ and TCGA datasets, ITGB2 expression was positively correlated with all targets (all correlation coefficients greater than 0, and all p < 0.05) (Figure 5B,C).

By reviewing relevant literature and referring to various gene sets collected by “IOBR”, we screened 16 gene sets associated with tumor immunotherapy response (Table S9). Through correlation analysis, we found that the expression of ITGB2 is positively correlated with these 16 TME features related to immunotherapy response, in both the GSE23400↗ and TCGA datasets (all p < 0.05, correlation coefficients > 0, Figure 5D, Table S10).

Utilizing miRWork, we identified 2796 miRNAs with a binding probability greater than 0.9 to ITGB2 (Figure 5E). Through the aforementioned screening methods, we ultimately identified three key miRNAs (hsa‐miR‐18a, hsa‐miR‐196a, hsa‐miR‐21, Figure 5F). In the GSE66274↗ dataset, compared to esophageal normal tissues, these three key miRNAs were downregulated in ESCC (all p < 0.001, Figure 5G). Similar results were observed in the GSE67268↗ dataset (all p < 0.0001, Figure 5H).