Spatial transcriptome profiling identifies DTX3L and BST2 as key biomarkers in esophageal squamous cell carcinoma tumorigenesis

Abstract

Background

Understanding the stepwise progression of esophageal squamous cell carcinoma (ESCC) is crucial for developing customized strategies for early detection and optimal clinical management. Herein, we aimed to unravel the transcriptional and immunologic alterations occurring during malignant transformation and identify clinically significant biomarkers of ESCC.

Methods

Digital spatial profiling (DSP) was performed on 11 patients with early-stage ESCC (pT1) to explore the transcriptional alterations in epithelial, immune cell, and non-immune cell stromal compartments across regions of distinct histology, including normal tissues, low- and high-grade dysplasia, and cancerous tissues. Furthermore, single-cell spatial transcriptomics was performed using the CosMx Spatial Molecular Imaging (SMI) system on 4 additional patients with pT1 ESCC. Immunohistochemical (IHC) analysis was performed on consecutive histological sections of 20 pT1 ESCCs. Additionally, public bulk and single-cell RNA-sequencing (scRNA-seq) datasets were analyzed, and in vitro and in vivo functional studies were conducted.

Results

Spatial transcriptional reprogramming and dynamic cell signaling pathways that determined ESCC progression were delineated. Increased infiltration of macrophages from normal tissues through dysplasia to cancerous tissues occurred. Macrophage subtypes were characterized using the scRNA-seq dataset. Cell–cell communication analysis of scRNA-seq and SMI data indicated that the migration inhibitory factor (MIF)-CD74 axis may exhibit pro-tumor interactions between macrophages and epithelial cells. DSP, SMI, and IHC data demonstrated that DTX3L expression in epithelial cells and BST2 expression in stromal cells increased gradually with ESCC progression. Functional studies demonstrated that DTX3L or BST2 knockdown inhibited ESCC proliferation and migration and decreased M2 polarization of tumor-associated macrophages.

Conclusions

Spatial profiling comprehensively characterized the molecular and immunological hallmarks from normal tissue to ESCC, guiding the way to a deeper understanding of the tumorigenesis and progression of this disease and contributing to the prevention of ESCC. Within this exploration, we uncovered biomarkers that exhibit a robust correlation with ESCC progression, offering potential new avenues for insightful therapeutic approaches.

Background

Esophageal cancer (EC) remains a major cause of cancer-related morbidity and mortality worldwide [1, 2]. Esophageal squamous cell carcinoma (ESCC) is the most common histological subtype of EC in Asian populations [3]. ESCC has an extremely poor prognosis and high mortality rate due to its rapid progression, treatment resistance, and high metastasis rate [4, 5]. Therefore, strategies to treat ESCC in its early or premalignant stages are warranted. However, the cellular and molecular mechanisms underlying early tumorigenesis remain largely unexplored, and the evolution of normal tissues to low-grade dysplasia, high-grade dysplasia, and eventually invasive cancers, is unclear.

Research targeting ESCC tumorigenesis has predominantly focused on ESCC tumor cell behavior [6,7,8,9]. However, the importance of the tumor microenvironment (TME) and its interaction with tumor cells during tumor initiation and progression has recently been recognized [10, 11]. The evolutionary advantage of cancer is majorly influenced by the overall potency of the cellular ecosystem rather than solely by tumor cells [12]. Furthermore, the constantly evolving TME should be comprehensively understood to develop effective therapeutic strategies as it may unveil the mechanisms of immune suppression and potential targets for investigation. However, the impact of the TME on ESCC progression remains unclear. Existing related studies have relied on bulk expression analysis [13,14,15] or single-cell RNA sequencing (scRNA-seq) [16,17,18,19]. These methods have inherent limitations, such as an inability to capture the comprehensive cellular composition and spatial context across distinct regions within tumor tissues. Consequently, a deeper understanding of the in situ intercellular communication and spatial niches that orchestrate tissue development and homeostasis has been hindered [20]. Advancements in spatially resolved transcriptomic techniques, such as in situ hybridization and spatial barcoding provide promising avenues for overcoming the limitations of bulk RNA-seq and scRNA-seq [21]. These techniques enable unbiased mapping of transcriptomic profiles while simultaneously capturing morphological changes within tissues [22]. NanoString digital spatial profiling (DSP) facilitates extensive characterization of tumor-immune landscapes and high-plex quantification of mRNA and proteins within a single formalin-fixed paraffin-embedded (FFPE) tissue section, specifically targeting regions of interest (ROIs) or areas of illumination (AOIs) [23]. A recent study utilizing DSP characterized ESCC progression [24]; however, a deeper investigation of molecular changes within each tissue compartment and a comprehensive exploration of immune-related alterations during ESCC progression were absent.

Analyses using DSP are particularly advantageous when studying early-stage ESCC, such as pT1 ESCC. Various histological entities, including normal tissues, distinct levels of dysplasia, and cancer, can coexist within a single lesion. This remarkable scenario presents a distinct opportunity to investigate the underlying biological processes that drive the initiation of ESCC tumorigenesis. Notably, histological entities originate from a shared source and possess an identical genetic background (germline) when exposed to comparable environments. Consequently, this model enables comprehensive exploration of the intricate mechanisms underlying ESCC development.

In this study, we investigated the transcriptional changes during ESCC onset using DSP to the matched normal tissues, low- and high-grade dysplasia, and cancers present within the same pT1 ESCC lesion. Genes and pathways that were consistently deregulated during ESCC progression were identified in the epithelial, immune cell, and non-immune cell stromal compartments. The TME composition exhibited marked changes during ESCC onset. A correlation between DTX3L, BST2, macrophages, and ESCC progression was revealed, respectively. In vitro and in vivo functional studies elucidated the roles of candidate genes in ESCC progression and M2 macrophage polarization. Overall, these findings provide valuable insights into the molecular characteristics underlying ESCC progression, identify potential targets for further investigation, and offer novel perspectives for developing effective therapeutic strategies.

Methods

Patient eligibility

We retrospectively analyzed patient FFPE tissue samples collected from 24 patients with pT1 ESCC who underwent complete surgical resection at the Department of Thoracic Surgery, First Affiliated Hospital of Soochow University. The eligibility criteria for inclusion were as follows: (a) absence of any prior preoperative treatment; (b) complete excision of the tumor via Ivor Lewis or McKeown esophagectomy, with routine lymph node dissection; (c) postoperative pathology confirmed as pT1 ESCC; and (d) no history of other intricate or malignant neoplasms. All patients underwent pathological staging according to the 8th edition of the tumor-node-metastasis (TNM) staging system established by the American Joint Committee on Cancer. Additional file 1: Table S1 summarizes the clinical characteristics of the patients.

Patient cohorts

Out of the eligible 24 patients, 11 were randomly selected for Digital spatial profiling (DSP). Immunohistochemical (IHC) analysis was performed on the consecutive sections of the same 11 pT1 lesions and an additional 9 pT1 ESCCs. Four additional patients were chosen to perform single-cell spatial transcriptomics using the CosMx Spatial Molecular Imaging (SMI) system. Please see Additional file 1: Table S1 for detailed information.

Public datasets including bulk RNA-seq and single-cell RNA-seq (scRNA-seq) data were downloaded. A total of 95 transcriptome profiles from patients with ESCC were obtained from The Cancer Genome Atlas (TCGA) cancer sample cohorts via the Xena data portal (https://xenabrowser.net/datapages). Furthermore, an ESCC dataset comprising ten tumor samples and ten adjacent normal tissues (GSE213565) [25] was obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). These two datasets were used to characterize the candidate genes. Additionally, three publicly available scRNA-seq datasets for ESCC were downloaded (Additional file 1: Table S2), namely GSE145370 [19], GSE160269 [17], and GSE221561 [26]. The GSE145370 dataset encompassing cells obtained from seven adjacent normal samples and seven ESCC tumor samples was used to characterize macrophage subpopulations. GSE160269 and GSE221561 datasets comprising six adjacent normal samples and eight ESCC tumor samples were used to investigate the crosstalk between macrophages and tumor epithelial cells.

Digital spatial profiling

DSP was performed using well-established methodologies (NanoString) [23]. FFPE Sects. (5-µm-thick) from eleven different patients (n = 11) were subjected to incubation with immunofluorescent antibodies and UV-photocleavable barcode-conjugated RNA in situ hybridization probes sourced from GeoMx Cancer Transcriptome Atlas (CTA; v.2.0; Additional file 1: Table S3). The prepared slides were stained with pan-cytokeratin (PanCK) (AE1 + AE3, Novus Biologicals, Cat# NBP2-33200 AF488) to identify epithelial cells, CD45 (D9M8I from CST, Cat# 13917, internally conjugated) to identify immune cells, and SYTO13 (Cat# 121303303) to label nuclei. Subsequently, the slides were loaded onto the GeoMx instrument, and regions of interest (ROIs) were selected based on the immunofluorescence images. Pathologists from the Department of Pathology of the First Affiliated Hospital of Soochow University confirmed the accuracy of the selected ROIs, which were meticulously chosen according to their histological features, including normal tissue, low- and high-grade dysplasia, and cancer. Auto-segmentation was achieved using tailored UV illumination masks, strategically creating areas of illumination (AOIs) that specifically released photocleavable tags. Subsequently, the segmented ROIs were categorized into distinct molecularly defined tissue compartments through fluorescent co-localization as follows: tumor (PanCK + /CD45 −), immune cell (PanCK − /CD45 +), and non-immune cell stromal (SYTO13 + /PanCK − /CD45 −) compartments. The cleaved barcodes from each AOI were collected in 96-well plates. The library preparation was performed per the manufacturer’s guidelines (NanoString). Sequencing was conducted on a NextSeq 500 platform using a PE50 kit (Illumina, San Diego, CA, USA). FASTQ files for all AOIs were processed as previously described [23]. In brief, the files underwent demultiplexing based on unique molecular identifiers (UMIs) and were subsequently aligned with spatial barcode sequences.

Spatial transcriptomic data analysis

The GeoMx NGS Pipeline [23] (v2.2.0.2) was used to transform the sequenced FASTQ files into DCC files. An extensive quality control (QC) process was applied to the data, including technical signal, technical background, probe evaluation, and normalization. To ensure data integrity, the technical signal QC was implemented, excluding AOIs with a < 80% alignment rate of the reads against the template sequence. The technical background involved three indicators: count of the no template control (NTC), count of negative probes, and AOI parameters. The NTC count was employed to identify and eliminate potential template contamination. AOIs with an NTC count > 1000 were excluded from the analysis. The overall technical signal level was determined by assessing the count of negative probes, with the threshold set at four. AOI parameters were evaluated based on the number of nuclei and surface area. To meet the QC standards, an AOI was required to possess a nuclei count > 20 or a surface area > 1600 µm². To ensure consistency across the AOIs, their sizes were standardized by cell number and area normalization. Thereafter, high-quality data were normalized using Quantile 3 (Q3) [27] and utilized for downstream analysis. Dimension reduction analysis was conducted using t-Distributed Stochastic Neighbor Embedding (tSNE). The edgeR [28] (v3.34.0) package was utilized for differential gene expression analysis with adherence to specific threshold criteria: false discovery rate (FDR) < 0.05 and an absolute log2 value (fold change [FC]) > 0.5. To gain insights into the biological processes and pathways associated with the differentially expressed genes (DEGs), enrichments with an adjusted p-value < 0.05 were considered significant. To determine the abundance of immune and stromal cells within the tumor (immune) microenvironment, deconvolution analysis was conducted using the SpatialDecon (v1.2.0) [29] and CIBERSORT algorithms. Spearman’s correlation was used to identify genes/pathways with increasing or decreasing trends during the stepwise transition from a normal esophagus to carcinoma. The regional parameter was treated as an ordinal variable (1: normal, 2: low-grade dysplasia, 3: high-grade dysplasia, and 4: carcinoma). Differential immune cells between different regions were identified using t-test. Additionally, the Pearson correlation was used to assess the correlation between stepwise transition-related genes and immune cell infiltration levels.

Gene set enrichment analysis

The log2 transformed normalized gene expression matrix underwent single-sample gene set enrichment analysis (ssGSEA) using Wiki pathways [30] (date 2023.10.10) downloaded from the WikiPathways database (https://data.wikipathways.org/current/gmt/). To ensure accuracy, pathways with < 2 genes available in the GeoMx CTA probe set were excluded. The WikiPathways database [30] was used to extract non-redundant biological pathways suitable for ssGSEA [31]. To investigate the transformation from normal tissue to carcinomas in PanCK and stromal segments, the Spearman correlation test was employed to correlate the pathways and distinct tissue regions. Here, the region parameter was treated as an ordinal variable (1 = normal; 2 = low-grade dysplasia; 3 = high-grade dysplasia; 4 = carcinoma).

Public RNA-seq data analysis

CIBERSORT and xCell algorithms were used to determine the potential infiltrating scores of immune and stromal cell types. The ESTIMATE algorithm [32] was utilized to calculate the immune and stroma scores. Pearson correlation was used to identify genes associated with BST2/DTX3L, with an adjusted p-value < 0.05 and correlation coefficient > 0.4. The Gene Ontology (GO) database was used to annotate the possible biological functions of BST2/DTX3L-related genes using the ClusterProfiler R package.

Single-cell RNA-seq data processing

All scRNA-seq analyses were performed using the Seurat package [33] (v5.1.0) in R (v4.2.2). Seurat default parameters were used unless otherwise specified. Genes expressed in < 3 cells, cells with > 20% mitochondrial genes, < 200 genes, relevant low quality, and potential doublets were filtered out. Next, the data matrix was normalized to 10,000 reads per cell using the NormalizeData function. The FindVariableGenes function was used to select variable genes, and the datasets collected from different samples were integrated using FindIntegrationAnchors and IntegrateData functions with the parameter “dims = 1:20” to remove batch effects. Subsequently, principal component analysis (PCA) was performed, and uniform manifold approximation and projection (UMAP) dimensionality reduction and visualization were performed based on the PCA results. The top 50 principal components were used for UMAP projection and clustering analysis. Resolution parameter set to 2 was applied during clustering with the aim to make sure that markers of different cell types were clearly expressed in different clusters. Cell types of different clusters were annotated according to the specific genes of different subgroups. Macrophages were rerun the PCA, UMAP projection, and clustering analyses. The top 20 principal components were used, and the resolution parameter was set as 0.4. The epithelial cells were extracted from the raw expression matrix, following which the same preprocessing steps were applied, and subsequent clustering was performed to obtain the subpopulation structures. Utilizing the top 10 principal components for UMAP projection, a resolution parameter of 0.4 was employed for clustering analysis. To avoid overclassification, raw clusters exhibiting minimal DEGs were merged. DEGs characterizing the clusters were identified through the default Wilcoxon rank-sum test in the FindAllMarkers function with default parameters.

Single-cell trajectory analysis

To infer the hierarchical organization of macrophages and establish their pseudo-time trajectory, the Monocle2 R package [34] (v.2.20.0) was used to calculate their differentiation trajectory. Macrophage subpopulations were extracted from the Seurat dataset, and shared nearest neighbor clustering and differential expression analyses were performed. Genes with adjusted p values below 0.05 and a logFC of the average expression above 0.5 were selected for Monocle2 to order the cells using DDRTree and reverse graph embedding. After determining the pseudo-time value arrangement and differentiation trajectory, the plot_cell_trajectory was used to illustrate the pseudo-time values along the differentiation trajectory. Branch-dependent gene regulation was identified using the BEAM function. The plot_genes_branched_heatmap function generated clustered heatmaps displaying the expression patterns of the top 50 most frequently occurring genes (hub genes) in branches 1, 2, 3, and 4.

Single-cell spatial transcriptomics and data processing

Single cell spatial transcriptomic profiling of FFPE tissue Sects. (4-μm-thick) from 4 pT1 ESCC patients was performed using the CosMx Spatial Molecular Imager (SMI), as previously described [35]. In brief, the CosMx Human 6 K Discovery Panel was used (Additional file 1: Table S4). Each reporter set contains 16 readout rounds featuring four distinct fluorophores, creating a 64-bit barcode design with a Hamming distance of 4 (HD4) and a Hamming weight of 4 (HW4) to ensure minimal error rates. Probe fluorescence was detected at subcellular resolution by the CosMx SMI instrument, and the signal was aggregated to identify the specific RNA molecule in each location.

SMI data analyses were performed using the Seurat [33] (v5.1.0), harmony [36] (v 1.2.1), and scDblFinder [37] (v 1.18.0) packages in R (v4.4.1). Cells with ≤ 10 genes or ≤ 20 UMI counts or ≥ 98% UMI counts were filtered out. The scDblFinder was used to remove potential doublets. The NormalizeData function normalized the data with the normalization.method set to RC. The RunBanksy function was employed to incorporate neighborhood information that distinguishes subtly different cell types stratified by microenvironment for SMI data, with the lambda set as 0.2, use_agf set as T, and the group set as fov. The RunHarmony was utilized to remove the batch effect. Dimension reductions were performed using PCA and UMAP. The top 20 principal components were used for UMAP projection and clustering analysis. The resolution parameter was set to 1 during clustering analysis. The Vlnplot function was used to show the expression of key genes.

Cell–cell communication analysis

The R package CellChat [38] v.1.6.1 was used to infer the interplay between epithelial cells or subpopulations of epithelial cells and macrophages across tumor and normal samples. CellChat objects for tumor and normal samples were separately created, subsequently merging them for comparative analysis to uncover key ligand-receptor pairs and signaling pathways between two groups and to detect and visualize cell-state-specific cell–cell interactions. The signaling information was obtained from the “Secreted Signaling” module of CellChatDB.human. Following the official procedure, the standardized counts were inputted into CellChat, and standard preprocessing steps were performed, including functions with standard parameter settings, such as identifyOverExpressedGenes, identifyOverExpressedInteractions, and projectData. Subsequently, the computeCommunProb, computeCommunProbPathway, and aggregateNet functions were employed to calculate the strength of information flow and communication probability between different cell groups for each ligand-receptor pair. The visualization methods utilized include rankNet (with the mode set to comparison and stacked set to T), netVisual_bubble, netVisual_aggregate (with the layout set to circle), and plotGeneExpression.

Immunohistochemistry

FFPE tissue Sects. (5-μm-thick) from 20 pT1 ESCC patients were deparaffinized and rehydrated and heat-mediated antigen retrieval was performed using tris-ethylenediaminetetraacetic acid. To mitigate potential non-specific binding, the slides were blocked with phosphate-buffered saline blocking buffer for 30 min at room temperature. The primary antibodies used in this study were anti-DTX3L (PA552708, 1:100, Thermo Fisher, USA), anti-BST2 (ab243229,1:200, Abcam, USA), anti-CD68 (M0876, 1:200, Dako, USA), anti-CD163 (ZM-0428, 1:200, ZSGB-BIO, China), anti-SLC1A5 (20350–1-AP,1:100, Proteintech, China), anti-HAVCR2 (ab241332, 1:1000, Abcam, USA), anti-SIRPA (ab191419, 1:200, Abcam, USA), and Anti-CD276 (ab227670, 1:100, Abcam, USA). Tissue sections were incubated with the primary antibody overnight at 4°C. Next, the slides were incubated with poly-HRP (Cat#21140, Thermo Scientific, USA) for 1 h, followed by development with DAB chromogen (Cat#K3468, DAKO, 1:50 dilution, Denmark). Hematoxylin was used as a counterstain. The slides were then dehydrated, mounted using Micromount (Cat#3801731; Leica, Germany), and evaluated by two pathologists. Consecutive slides were used to detect different antibodies. The simplified semi-quantitative IHC scoring method [31, 39] was implemented for DTX3L and BST2. Briefly, ten high-power fields (HPFs) were randomly selected and assessed. The percentage of positive cells was graded as follows: 1, 0–24% positive cells; 2, 25–49% positive cells; and 3, 50–74% positive cells, and 4, 75–100% positive cells. The staining intensity was estimated and scored as follows: 0, absence of staining; 1, weak staining; 2, moderate staining; and 3, strong staining. Subsequently, the IHC score was derived by combining of these two scales (staining intensity score × the proportion score) and categorized as follows: 0–1, negative expression; 2–3, weak expression; 4–6, intermediate expression; and 8–12, strong expression. For BST2, within the stromal compartment of the tumor regions, scores exceeding the mean value (8.85) were classified into the high-expression group, those below the mean value were assigned to the low-expression group. The CD68, CD163, HAVCR2, SIRPA, and CD276 positive stromal cells were quantified in ten random HPFs from different regions and presented as the average number per HPF.

Cell lines and culture

ESCC cell lines (Eca-109, TE-1 and AKR) and the human monocytic leukemia cell line (THP-1) were obtained from the Shanghai Institute of Cell Biology Cell Bank (Chinese Academy of Medical Science, Shanghai, China) and cultured in RPMI-1640 medium (Gibco, USA) supplemented with 10% fetal bovine serum (Corning) and 1% penicillin/streptomycin (Gibco). All cell lines were cultivated in a humidified incubator at 37 °C with 5% CO₂ according to standard protocols. All cell lines were routinely tested for mycoplasma contamination, and the results were consistently negative.

Cell transfection and treatment

Cells were transfected with siRNAs (RiboBio, Guangzhou, China) using Lipofectamine RNA iMAX (ThermoFisher Scientific, USA) transfection reagent with OptiMEM (ThermoFisher Scientific, USA) according to the manufacturer’s instructions. Cells were collected 24 h after transfection. The lentiviral shRNA plasmid was co-transfected with packaging vectors (pMD2G and psPAX2) into HEK293T cells (CoruesBio, Nanjing, China). After transfection, virus particles were harvested at 24 and 48 h and filtered using a 0.45-µM filter. To establish knockdown stable cell lines, cells were infected with the designated virus particles along with 8 µg/mL polybrene. Virally infected cells were selected with puromycin to obtain stable DTX3L knockdown ESCC cell lines. THP-1-derived macrophages were treated with 20 ng/mL anti-BST2 (1221901–33-2, MCE, USA). Additional file 1: Table S5 lists the siRNA and shRNA sequences.

Transwell assay

To assess cell migration, a Transwell chamber (Millipore, Billerica, MA, USA) was positioned within a 24-well plate. The lower chamber received 600 μL of medium containing 10% fetal bovine serum, whereas the upper chamber was seeded with 20,000 cells in 200 μL of serum-free medium. After 24 h, cells adhering to the membrane were fixed with methanol and stained with crystal violet. Images were captured using an inverted fluorescence microscope for cell counting. Three random fields were selected and observed under the microscope to quantify the migratory cell population.

Wound healing assay

Transfected Eca-109 and TE-1 cells were planted in 6-well plates and incubated at 37 °C until full confluence. Upon reaching 80–90% confluence, a sterile pipette tip was utilized to create a linear scratch in one direction. Any detached cells were washed away from the monolayer. The cells were exposed to serum-free 1640/DMEM for 24–48 h, and images of the healing scratch wounds were captured using an inverted microscope.

Colony formation, cell counting Kit-8, and 5-ethyl-2'-deoxyuridine assays

For the colony formation assay, transfected Eca-109 and TE-1 cells were seeded in 6-well plates (500 cells/well) and cultured for 2 weeks. Cells were then fixed with methanol for 30 min and stained with crystal violet for 30 min. Cell colonies were then quantified.

For the Cell Counting Kit-8 (CCK-8) assay, transfected Eca-109 and TE-1 cells were seeded into 96-well plates (4000 cells/well), and cell viability was examined at different time points using CCK-8 assay (Beyotime Biotechnology, Shanghai, China). Cell absorbance at 450 nm was measured using a microplate reader (Bio-Rad, Hercules, USA).

The EdU DNA Proliferation Kit (KeyGene, China) was utilized to measure cell proliferation. Eca-109 and TE-1 cells were grown in 96-well plates in full medium until 80% confluence and then treated with 50 μM EdU for 6 h before analysis.

Human macrophage activation

THP-1 cell differentiation was induced via 6-h exposure to 185 ng/mL phorbol 12-myristate 13-acetate (PMA, S1819 Beyotime) in dimethylsulfoxide (DMSO). Cells were then polarized toward the M2 phenotype via incubation with 20 ng/mL interleukin (IL)-4 (P5129, Beyotime) and 20 ng/mL IL-13 (P5178, Beyotime) for 48 h in the presence of PMA.

RNA extraction and quantitative real-time PCR (qRT-PCR) analysis

Total RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The RNA concentration and purity were detected using a NanoDrop spectrophotometer (ND-2000, Thermo). Reverse transcription was performed using the ABScript II RT Mix for qPCR (ABclonal, China). Quantitative real-time fluorescence PCR was performed using the SYBR Green reagent (Vazyme Biotech, Nanjing, China). Expression results were detected on an ABI StepOne Plus qRT-PCR instrument. The 2ΔΔCt method was used to calculate the relative expression of target genes. The mRNA level was normalized to the housekeeping gene GAPDH. The sequences of qRT-PCR primers are listed in Additional file 1: Table S5.

Flow cytometry analysis

For cell surface flow cytometry staining, cells were stained with fluorescently labeled antibodies specific for surface proteins (1:100) in Cell Staining Buffer (BioLegend, USA) for 30 min at 4°C. Flow cytometry analysis (BD FACSCalibur, USA) was performed for data acquisition, and the data were analyzed using FlowJo software (version 10.5.3). CD163-FITC was purchased from BioLegend (333617, San Diego, California, USA).

Animal experiment

To assess the in vivo impact of DTX3L, female C57BL/6 mice (4–6 weeks old) were housed according to the Institutional Animal Guidelines of First Affiliated Hospital of Soochow University and randomly divided into three groups. Animal experiments were conducted with the approval of the Institutional Animal Care and Use Committee of the First Affiliated Hospital of Soochow University. Mice were subcutaneously injected with AKR cells that had been stably transfected with either sh-NC or sh-dtx3l-1/sh-dtx3l-2 (1 × 10⁶ cells). Tumor growth was monitored over a period of 28 days. Tumor volume was calculated using the standard formula a² × b × 0.5 (where a and b represent the short and long diameters of the tumor, respectively). Tumor measurements were taken every 7 days. The mice were kept in a positive pressure barrier facility at 20–26 ℃ temperature and 40–70% humidity, with 12-h light and dark lighting. After 4 weeks of subcutaneous tumor loading, or the weight loss of the mice reached 20–25%, or the appearance of cachexia and wasting symptoms, or the size of the solid tumor exceeded 10% of the animal’s body weight, the mice were euthanized with carbon dioxide, and the tumors were removed for IHC analysis.

Results

Spatially resolved transcriptomic profiling of ESCC tumorigenesis

To examine the molecular alterations occurring during ESCC onset and progression, GeoMx DSP technology was employed to profile 44 multistage esophageal resection specimens derived from 11 patients with pT1 ESCC. The specimens included normal tissues (N) (n = 11), low-grade dysplasia (L) (n = 11), high-grade dysplasia (H) (n = 11), and tumors (T) (n = 11). Representative DSP-analyzed tissue samples with hematoxylin and eosin staining or fluorescent labeling are shown in Fig. 1A and B. Three ROIs were selected for each sample (Fig. 1A, B), and each was segmented into molecularly defined tissue compartments via fluorescent co-localization using a masking and segmentation strategy to specifically examine the transcriptional alterations within each compartment, including the epithelial (PanCK + /CD45 −), immune cell (PanCK − /CD45 +), and non-immune cell stromal (PanCK − /CD45 − /SYTO13 +) compartments (Fig. 1C). In total, 373 AOIs, including 129 epithelial (33 N, 33 L, 33 H, and 30 T), 118 immune cell (29 N, 30 L, 30 H, and 29 T), and 126 non-immune cell stromal (31 N, 32 L, 32 H, and 31 T) AOIs, met the QC criterion (Additional file 1: Table S6 and S7). Upon combining the immune cell and non-immune cell stromal compartments as “stroma”, a total of 244 (60 N, 62 L, 62 H, and 60 T) AOIs were used for analysis (Additional file 1: Table S7).

Fig. 1

Digital spatial profiling analysis of early-stage ESCC. A Representative example of a pT1 ESCC specimen stained with H&E. Selected ROIs with different histological features are annotated. B Immunofluorescent detection of PanCK and CD45 on the same sample shown in A. C Representative ROI and AOIs annotated based on histology by pathologists and immunofluorescence staining with morphological markers and the compartmentalized image created by fluorescence co-localization. One 300 μM ROI per core was selected for DSP; panCK (green), CD45 (yellow), and SYTO13 (blue). D Dimensionality reduction visualization of all AOIs according to overall gene expression profiles by tSNE. tSNE plots are annotated by segment. E Heatmaps of gene expression (n = 1825) using unsupervised clustering for PanCK AOIs (n = 129, left), CD45 AOIs (n = 118, middle), and non-immune cell stromal AOIs (n = 126, right). Heatmaps are annotated by histological region and sample ID. ESCC, esophageal squamous cell carcinoma; H&E, hematoxylin–eosin staining; ROI, region of interest; AOI, areas of illumination; PanCK, pan-cytokeratin; DEGs, differentially expressed genes; tSNE, t-Distributed Stochastic Neighbor Embedding; FDR, false discovery rate

The tSNE for data exploration via dimension reduction showed that the AOIs were distinctly segregated based on segments (Fig. 1D; Additional file 2: Fig. S1A). The genes in the CTA collection (Additional file 1: Table S3) exhibited differential expression between segments within one sample type or all samples (Additional file 1: Table S8 and S9; Additional file 2: Fig. S1B, C, and S2). For a comprehensive overview, the expression of all 1825 CTA genes in unsupervised clustered heatmaps was subsequently visualized per segment (Fig. 1E). During tumorigenesis, 240, 54, and 238 genes were upregulated in epithelial, immune cell, and non-immune cell stromal compartments, respectively (Spearman’s correlation; FDR < 0.05; Additional file 2: Fig. S3; Additional file 1: Table S10). Notably, a larger number of genes demonstrated decreased expression during ESCC progression, particularly in the epithelial segment, and were generally consistent across samples (Fig. 1E; Additional file 2: Fig. S3), reflecting the shared loss of physiological characteristics in tissues during malignant transformation.

Distinct biological processes associated with ESCC onset between the epithelium and stroma

To comprehensively evaluate the changes that occur from normal tissue to carcinoma during progression, pathway enrichment analysis was conducted using a curated dataset obtained from the WikiPathways database [30]. Based on the pathway enrichment scores (Additional file 1: Table S11), the majority of cancer and high-grade dysplasia AOIs of the epithelial segment clustered together (Additional file 2: Fig. S4), highlighting clear differentiation from the histology of normal tissues and low-grade dysplasia. Seven pathway clusters (C1–C7) were identified with different degrees of association with malignant transformation (Fig. 2). Two clusters (C1 and C2) grouped pathways that were mainly active in the epithelial segment. C1 pathways included biological processes associated with cell proliferation and DNA damage repair, which consistently increased during ESCC onset. Most pathways clustered in C2 were associated with mitochondrial function, which decreased during malignant transformation in the epithelial segment. Pathways in C5 and C6 were associated with the stromal compartment. Pathways in C5 reflected immune-related pathways that were deregulated during tumorigenesis. Pathways clustered in C3, C4, and C7 were active in the epithelium and stroma, demonstrating variable associations with tumor progression.

Fig. 2

Key biological pathways associated with advancing histology in epithelial and stromal compartments. AOIs are ordered by segment and region. Unsupervised clustering of ssGSEA was calculated using gene sets from the WikiPathways database. All pathways with a significant association with histology, either in PanCK or stromal segments, are included (FDR < 0.05). Main clusters of identified pathways are indicated with white horizontal lines. The corresponding Spearman correlation coefficient, Rho, between histology as ordinal variable and the enrichment scores are indicated for each pathway in PanCK and stromal segments separately (dotted heatmap). FDR is calculated using Benjamini–Hochberg method. AOIs, areas of illumination; ssGSEA, single sample gene set enrichment; PanCK, pan-cytokeratin; FDR, false discovery rate

Immune-related alterations during ESCC progression

To delineate immune-related alterations, we examined the expression of 33 functionally characterized immunomodulatory genes associated with cancer immunity [40, 41]. The majority of immunomodulatory gene expression was derived from stromal segments (Additional file 2: Fig. S5A). Notably, inhibitory immune checkpoint genes (e.g., HAVCR2, CD276, and SIRPA) expressed by stromal cells, exhibited increased expression in the stroma of carcinoma regions (Additional file 2: Fig. S5A, B). These findings were further corroborated by IHC data (Additional file 2: Fig. S5C, D). These upregulated inhibitory immune checkpoint genes may serve as potential therapeutic targets.

Next, we re-assessed the altered enrichment of WikiPathways, specifically targeting the gene sets associated with immunity (Additional file 1: Table S12). Remarkably, immune-related signatures underwent alterations during malignant transformation, reflecting changes in the immune microenvironment (Fig. 3A). Cytokine-related (e.g., overview of interferons, type II interferon, interleukin (IL)-6, and IL-10) and inflammatory response pathways were upregulated in the stromal compartment during the transition from normal mucosa to carcinoma. Two pathways linked to innate immunity, including the “toll-like receptor signaling pathway” and “fibrin complement receptor 3 signaling pathway,” were consistently upregulated in the non-immune cell stromal segment.

Fig. 3

Immune-related alterations related to ESCC histology. A Forest plot of Spearman’s Rho and corresponding 95% CI for the correlation between the enrichment score of immune-related pathways from WikiPathways and the histology as the ordinal variable. Correlations were assessed in PanCK, CD45, and non-immune cell stromal segments separately. B Boxplots showing the proportions of distinct immune and stromal cell populations calculated by SpatialDecon across histologies in the CD45 and non-immune cell stromal segments. The centerline represents the median, and box limits represent the upper and lower quartiles. Each dot represents an AOI. C Boxplots showing the proportions of M2 macrophages calculated by CIBERSORT across histologies in the stromal segments. D Representative images of IHC-stained slides of CD68 and CD163 in distinct histologies. Scale bar = 50 μm. E Staining statistics of CD68 and CD163 by IHC staining in the four pathological stages in a validation set of 20 pT1 ESCC samples. CI, confidence interval; AOI, area of illumination; IHC, immunohistochemistry

Intrigued by the modifications observed in immune-related genes and pathways during ESCC onset, alterations in the tumor (immune) microenvironment were further evaluated between regions of distinct histology. First, immune cell deconvolution was performed using SpatialDecon to estimate the relative abundance and proportions of specific cell subsets. The relative proportions of CD4 + T cells, neutrophils, and mast cells in the immune cell segment decreased during the transformation from normal tissues to carcinomas (Fig. 3B). Conversely, the relative proportions of macrophages and regulatory T cells in the immune cell compartment and fibroblasts in the non-immune cell compartment increased during tumorigenesis (Fig. 3B). Macrophages potentiate the seeding and establishment of metastatic cells and function in tumor initiation, progression, and metastasis [42]. An increase in the proportion of pro-tumor M2 macrophages was observed during ESCC progression (Fig. 3C). The IHC results further demonstrated increased expression of CD68 (macrophage marker) and CD163 (M2 macrophage marker) across various disease stages, from normal tissues to ESCC (Fig. 3D, E; Additional file 1: Table S13). Overall, these findings strongly suggested the emergence of an immuno-suppressed microenvironment coinciding with tumor progression.

Characterization of macrophage subpopulations and crosstalk between macrophages and tumor epithelial cells

Tumor onset and progression may be influenced by macrophage phenotypes, with distinct macrophage subtypes playing pivotal roles in facilitating tumor-promoting activities [43]. Thus, a publicly available single-cell transcriptomic dataset (GSE145370) [19] was analyzed (see Additional file 1: Table S2 for detailed clinical information). The tumor group had a significantly higher proportion of total macrophages than the normal group (p < 0.05; Fig. 4A), consistent with our DSP deconvolution results (Fig. 3B). Five macrophage subclusters were identified (Fig. 4B). According to the review by Ma et al. [44] and based on their highly expressed signature genes and predicted functions (Additional file 2: Fig. S6A and S7), these subsets were classified as inflammatory cytokine-enriched (Macro.0), uncharacterized (Macro.1), interferon-primed (Macro.2), lipid-associated (Macro.3), and proliferating (Macro.4) macrophages. The composition ratio of the macrophage subclusters varied greatly between the normal and tumor groups (Fig. 4C). Notably, the tumor group had a higher ratio of Macro.0 and Macro.4, whereas the normal group had a higher ratio of Macro.1. The proportions of the remaining macrophage subclusters did not differ significantly between the groups.

Fig. 4

Characterization of macrophage subtypes and crosstalk between macrophages and tumor epithelial cells. A Total macrophages in tumor and normal tissues. B UMAP visualization of macrophage subclusters. C Percentages of identified macrophage subclusters in normal/tumor samples. D Trajectory order of the macrophage subpopulations by pseudotime value. E Distribution of macrophage subpopulations on the developmental tree by clusters. F Top representative significantly differentially expressed ligand-receptor signaling pathways between normal and tumor samples. G Bubble plot showing the selected ligand-receptor interactions between epithelial cells and macrophages with p-values < 0.01. H Bubble plot showing the MIF-(CD74 + CXCR4) and MIF-(CD74 + CD44) interactions between epithelial cells and macrophages using SMI data. I UMAP visualization of epithelial cell subclusters. J Bubble plot showing the MIF-(CD74 + CXCR4) and MIF-(CD74 + CD44) interactions between epithelial cell subclusters and macrophages. UMAP, uniform manifold approximation and projection; MIF, macrophage migration inhibitory factor; SMI, Spatial Molecular Imaging

Next, the transitional states of macrophages were explored through trajectory analysis. Based on macrophage gene expression dynamics, a pseudo-time developmental tree was constructed, revealing four distinct branch points (Fig. 4D). Five macrophage subclusters were dispersed across various branches within the developmental tree (Fig. 4E). Macro.4 exhibited the lowest pseudo-time value and occupied the initial position in the developmental tree, indicating its role as a developmental origin for the other subclusters. This was consistent with the GO enrichment analysis of Macro.4, labeling it as a population with proliferative characteristics (Additional file 2: Fig. S7). Macro.1 and Macro.2 were at similar branches, whereas Macro.0 and Macro.3 exhibited more widespread distributions. The genes with the most significant changes at the four branches were clustered (Additional file 2: Fig. S6B).

A comprehensive understanding of the interactions between tumor epithelial cells and macrophages within the TME will provide insight into the mechanisms underlying ESCC progression. Recent studies have highlighted the significance of communication between tumor and immune cells in ESCC development [18]. As the aforementioned scRNA-seq dataset focused on CD45-sorting cells, data pertaining to epithelial cells was lacking. Accordingly, two additional public ESCC scRNA-seq datasets, namely GSE160269 [17] and GSE221561 [26] (see Additional file 1: Table S2 for detailed clinical information) were used to investigate the signaling network between epithelial cells and macrophages in normal and tumor tissues. However, subclustering the macrophage populations was not feasible due to the limited number of macrophages. Thus, the signaling network between epithelial cells and total macrophages was elucidated. Complex cell–cell interaction networks were observed between epithelial cells and macrophages (Fig. 4F). Ten signaling pathways associated with crosstalk between epithelial cells and macrophages were exclusively identified in tumor samples. In tumor tissues, seven signaling pathways exhibited significant engagement in the communication between epithelial cells and macrophages compared with normal tissues. Importantly, all these pathways have been reported to be related to tumorigenesis. Additionally, single-cell spatial transcriptomics was performed on FFPE specimens from an additional 4 pT1 ESCCs using the CosMx SMI technique (Additional file 1: Table S1 and S4). Notably, scRNA-seq and SMI data revealed that the ligand-receptor pairs macrophage migration inhibitory factor (MIF)-(CD74 + CXCR4) and MIF-(CD74 + CD44) exhibited particularly high interaction scores (Fig. 4G, H).

To identify the epithelial cell subset(s) responsible for the crosstalk between the MIF signaling pathway and macrophages, the epithelial cells were further subdivided into 6 subsets (Epi.0-Epi.5; Fig. 4I; Additional file 2: Fig. S8A). The distribution of these subclusters showed a marked distinction between the tumor and normal groups (Additional file 2: Fig. S8B). Specifically, the proportions of Epi.0, Epi.2, Epi.3, and Epi.5 were higher in the tumor group, whereas those of Epi.1 and Epi.4 were higher in the normal group. The MIF signaling pathway showed robust interplay between macrophages and each epithelial subcluster (Fig. 4J; Additional file 2: Fig. S8C). Subpopulations of epithelial cells highly expressed MIF (Additional file 2: Fig. S8D). MIF, an inflammatory cytokine, is involved in the tumorigenesis of various cancer types [45,46,47]. CD74 promotes tumor cell proliferation by interacting with MIF [48, 49]. Consequently, inhibiting the MIF-CD74 axis holds promise for impeding tumorigenesis and may serve as a viable therapeutic strategy for treating ESCC.

Identification of potential biomarkers associated with stepwise ESCC progression

To elucidate the gene expression changes associated with the stepwise progression from normal tissue to carcinoma, the differential gene expression was analyzed between regions with distinct histology types. Within the epithelial region, comparisons between normal tissue and low dysplasia areas and between low- and high-grade dysplasia areas yielded 364 and 89 DEGs, respectively (Fig. 5A; Additional file 2: Fig. S9A; Additional file 1: Table S14). Only minor discrepancies were observed between normal tissue and low-grade dysplasia and between low- and high-grade dysplasia in the stromal compartment (Fig. 5B, C; Additional file 2: Fig. S9B; Additional file 1: Table S14). Notably, C-X-C Motif Chemokine Ligand 1 (CXCL1) and CCL18 expression were increased in the stroma of low-grade dysplasia (Fig. 5C; Additional file 2: Fig. S9B), indicating early inflammatory processes associated with enhanced tumorigenesis [50, 51]. In the tumor areas, substantial differences in gene expression were observed in the epithelial and stromal compartments (Fig. 5A–C, right panels; Additional file 2: Fig. S9A, B, right panels; Additional file 1: Table S14). In general, these findings endorse the notion that ESCC initiation is linked to transcriptional changes in (pre-) malignant epithelial cells and adjacent stromal cells.

Fig. 5

Identification of DEGs related to the stepwise progression of ESCC in spatially distinct compartments. A–C Volcano plots showing DEGs between distinct histologies in stepwise comparisons. Reference groups for comparisons are indicated. The analysis was performed for the epithelial segment (PanCK) (A), immune cell segment (CD45) (B), and non-immune cell stromal segment (C). − log10 (FDR-adjusted p-value) increased with statistical significance and is indicated on the y-axis. The dotted horizontal line represents the adjusted p-value cutoff (FDR < 0.05). The dotted vertical lines represent the log2 (fold change) cutoffs (|FC|> 1.5). Blue and red points are high-expression genes in the corresponding group that achieved statistical significance. DEGs with the highest and lowest log2FC are labeled. D Boxplots of log2-transformed normalized gene expression of four DEGs detected in epithelium and corresponding gene expression in immune cell and non-immune cell stromal segments. E Boxplots of four DEGs identified in immune cell or non-immune cell segments and corresponding gene expression in the epithelial segment. F Representative images of IHC detection for DTX3L and BST2 across regions with distinct histology (normal, low-grade dysplasia, high-grade dysplasia, and tumor). Scale bar = 50 μm. G Protein abundance of DTX3L and BST2 in a validation set of 20 pT1 ESCC samples. The stacked bar chart reflects the proportion of samples in each scored category. Significance is indicated by asterisks: *p < 0.05 and **p < 0.01. ESCC, esophageal squamous cell carcinoma; DSP, digital spatial profiling; IHC, immunohistochemistry; DEGs, differentially expressed genes; PanCK, pan-cytokeratin; FDR, false discovery rate; Log2FC, log2 fold change

Specific genes exhibited clear alterations in expression between neighboring regions of distinct histology, thereby emerging as potential indicators for early ESCC detection (Additional file 1: Table S15). To illustrate this, eight genes with markedly differential expression between regions were selected (Fig. 5D, E; Additional file 2: Fig. S9C, D), each with significant implications for tumor progression. ALCAM [52], PPL [53], SGK1 [54], SLC1A5 [55], TFRC [56], and TPM4 [57] are associated with ESCC tumorigenesis. ESCC specimens exhibited higher SLC1A5 expression in the epithelial compartment than in the corresponding adjacent normal tissue, whereas no such difference was observed in the stromal compartment (Fig. 5D; Additional file 2: Fig. S9C, E), consistent with the findings of Zheng et al. [58]. The expression of the potentially novel biomarkers, DTX3L in the epithelial segment and BST2 in the stromal segment, exhibited significant upregulation (Fig. 5D, E; Additional file 2: Fig. S9C, D). SMI data further confirmed that DTX3L expression in epithelial cells and BST2 expression in stromal cells increased with ESCC progression (Additional file 2: Fig. S9F). To ascertain whether the observed modifications in gene expression of these novel biomarkers correlated with alterations in protein levels, IHC staining was performed (Additional file 1: Table S1 and S13). Upregulation of DTX3L in the epithelial segment and of BST2 in the stromal segment during the stepwise progression of ESCC was validated across the samples (Fig. 5F, G; Additional file 1: Table S13).

To summarize, our approach successfully identified spatial markers that can be utilized as clinical biomarkers.

DTX3L and BST2 characterization

DTX3L, alternatively referred to as B-lymphoma and BAL-associated protein, functions in tumor progression in glioma, melanoma, and pancreatic cancer [39, 59, 60]. Bone marrow stromal cell antigen 2 (BST2, also known as HM1.24/CD317) is associated with tumor progression and metastasis in breast cancer, colorectal cancer, and pancreatic endocrine tumors [61,62,63]. Nevertheless, the functions of DTX3L and BST2 in ESCC have not previously been reported. We first investigated their expression across multiple cancer types using data from TCGA (Additional file 2: Fig. S10A, B). DTX3L and BST2 were broadly expressed, and significantly increased in tumor tissues of various cancer types. The GEO dataset (GSE213565) [25] demonstrated that DTX3L and BST2 expressions were significantly higher in ESCC tissues than in normal tissues (Additional file 2: Fig. S10C). TCGA dataset demonstrated DTX3L and BST2 expression did not significantly correlate with sex or the TNM stage of ESCC (Additional file 2: Fig. S10D, E). GO biological process analysis revealed that the top GO terms enriched in DTX3L-high and BST2-high ESCC were regulation of receptor catabolic process and interferon response, respectively (Additional file 2: Fig. S10F, G).

BST2 plays a crucial role in tumor progression and polarization of M2 tumor-associated macrophages (TAMs) in colorectal cancer [63]. Thus, whether BST2 had a similar biological function in ESCC was explored. ESCC patients with high BST2 expression had higher immune scores (Fig. 6A). xCell results showed that the BST2 high-expression group had more significant immune cell infiltration than the BST2 low-expression group (Fig. 6B). Total and M2 polarized macrophages were significantly enriched in the BST2 high-expression group. Additionally, analysis of the DSP data revealed a positive correlation between BST2 expression in the immune cell and stromal compartments and the total and M2 macrophages (Fig. 6C, D). The representative IHC images showed that high BST2 expression tended to be associated with increased M2 macrophages (CD163 positive) infiltration (Fig. 6E, F). Notably, scRNA-seq data revealed that total macrophages, Macro.1, and Macro.2 had significantly higher BST2 expression in tumor tissues than in normal tissues (p < 0.05) (Fig. 6G, H). Interestingly, interferon-primed macrophages (Macro.2) resembled a recently identified immunosuppressive macrophage subtype, which exerts immunosuppressive effects through tryptophan degradation and immunosuppressive Treg recruitment [44, 64]. Overall, BST2 may regulate the infiltration and state of TAMs in the ESCC microenvironment.

Fig. 6

Association between BST2 and TAM infiltration in ESCC. A Association between immune and stroma scores and BST2 expression in the TCGA cohort. B xCell analysis of immune cell infiltration between BST2 high- and low-expression groups in the TCGA cohort. C, D Correlation analysis between BST2 expression and proportions of total macrophages (C) and M2 macrophages (D) in the stromal segment using DSP data. E Representative IHC staining images of BST2, CD68, and CD163 (M2 macrophage marker) in the high- and low-BST2 expression groups. Scale bar = 50 μm. F Quantitative analysis of IHC staining. G, H Violin plots showing the distribution of BST2 expression in macrophages (G) and macrophage subclusters (H) of normal/tumor samples determined by analyses of the scRNA-seq dataset. Student’s t test is used to compare differences. ns p > 0.05, *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001. ESCC, esophageal squamous cell carcinoma; TAMs, tumor-associated macrophages; TCGA, The Cancer Genome Atlas; DSP, digital spatial profiling; IHC, immunohistochemistry

Identifying the roles of DTX3L and BST2 in accelerating ESCC progression

Functional investigations can offer insights into the molecular mechanisms underlying the influence of genes on promoting cancer progression. To estimate the impact of DTX3L in ESCC progression, DTX3L knockdown cell lines (Eca-109 and TE-1) were constructed to evaluate the effect of DTX3L on ESCC using the shRNA and siRNA techniques, respectively. Transwell and wound healing assays showed that DTX3L knockdown markedly inhibited ESCC cell migration (Fig. 7A, B; Additional file 2: Fig. S11A, B and S12A). Furthermore, DTX3L knockdown inhibited Eca-109 and TE-1 cell proliferation as determined by colony formation (Additional file 2: Fig. S11C) as well as CCK-8 and EdU assays (Fig. 7C, D; Additional file 2: Fig. S11D, E). To explore the biological role of DTX3L in vivo, Dtx3l knockdown and control AKR cells were inoculated in C57BL/6 mice, respectively. A signification reduction in tumor volume upon dtx3l knockdown was observed (Fig. 7E, F). Additionally, we utilized IHC to investigate the association between Dtx3l expression and proliferation (detected by Ki-67), as well as macrophage infiltration (detected by F4/80 and CD206). The results revealed a decrease in proliferation within the Dtx3l knockdown group, accompanied by a reduction in infiltrated macrophages, particularly M2 macrophages (Fig. 7G, H).

Fig. 7

Identifying the role of candidate genes in accelerating ESCC progression and M2 macrophage polarization in vitro and in vivo using the shRNA system. A Migration of DTX3L knockdown Eca-109 or TE-1 cells assessed by Transwell assay. Scale bar = 100 μm. B Migration of DTX3L knockdown Eca-109 or TE-1 cells assessed by wound healing assay. C CCK-8 assay of Eca-109 or TE-1 cells with DTX3L knockdown. D Proliferation of DTX3L knockdown Eca-109 or TE-1 cells, as determined by EdU incorporation assay. DAPI was used to stain the nuclei. Scale bar = 50 μm. E Pictures of tumors after mice euthanasia at 28 days. F Tumor volume of C57BL mice subcutaneously injected with AKR cells with Dtx3l knockdown or with sh-NC. G Representative IHC images displaying the detection of ki-67, F4/80, and CD206. Scale bar = 50 μm. H Statistics analysis of ki-67, F4/80, and CD206 by IHC staining. I Relative mRNA expression of macrophage markers (Arg1, Fizz1, and Mgl1) in THP-1-derived macrophages with BST2 knockdown. J Flow cytometry of CD163 in THP-1-derived macrophages with BST2 knockdown. K Migration ability of Eca-109 and TE-1 cells was assessed by Transwell assay (scale bar = 100 μm). Cells were cocultured with CM from THP-1-derived macrophages with BST2 knockdown. L Migration of DTX3L knockdown Eca-109 or TE-1 cells assessed by wound healing assay. Cells were cocultured with CM from THP-1-derived macrophages with BST2 knockdown. M Proliferation of Eca-109 and TE-1 cells was detected by CCK-8 cocultured with CM from THP-1-derived macrophages with BST2 knockdown. N Proliferation of Eca-109 and TE-1 cells detected by EdU assay. Cells were cocultured with CM from THP-1-derived macrophages with BST2 knockdown. Scale bar = 50 μm. O Relative mRNA expression of macrophage markers (Arg1, Fizz1, and Mgl1) in THP-1-derived macrophages cocultured with CM from Eca-109 and TE-1 cells with DTX3L knockdown. P Flow cytometry of CD163 in THP-1-derived macrophages cocultured with CM from Eca-109 and TE-1 cells with DTX3L knockdown. **p < 0.01, and ***p < 0.001. CCK-8, Cell Counting Kit-8; Edu, 5-Ethyl-2'-deoxyuridine; CM, conditioned medium

To assess the effect of BST2 on macrophage polarization, BST2 was knocked down in THP-1-derived macrophages, or THP-1-derived macrophages were treated with anti-BST2. As expected, M2 marker expression was significantly lower in the BST2-knockdown or anti-BST2-treated macrophages than in the controls (Fig. 7I; Additional file 2: Fig. S11F). These results were further supported by flow cytometry analysis, revealing a decrease in the number of CD163⁺ cells (Fig. 7J; Additional file 2: Fig. S11G). Considering that the phenotypic transition of macrophages regulates cancer cell proliferation, migration, and invasion capabilities, the effects of BST2 knockdown or anti-BST2 treatment of macrophages on ESCC progression were investigated. Conditioned medium (CM) collected from THP-1-derived macrophages with BST2 knockdown or treated with anti-BST2 was added to Eca-109 or TE-1 cell cultures, and the proliferation and migration abilities of Eca-109 or TE-1 cells were evaluated. ESCC cells cultured with the CM from BST2-knockdown or anti-BST2-treated macrophages exhibited inhibited migration (Fig. 7K, L; Additional file 2: Fig. S11H, I, and S12B) and proliferation (Fig. 7M, N; Additional file 2: Fig. S11J-L). Additionally, THP-1 cells cocultured with the CM from Eca-109 and TE-1 cells with DTX3L knockdown promoted M2 macrophage polarization (Fig. 7O, P; Additional file 2: Fig. S11M, N).

To ensure the reliability of our findings, we conducted a rescue experiment to address potential off-target effects resulting from the shRNA technique. This rescue experiment involved the overexpression of either DTX3L or BST2 in cells where DTX3L or BST2 has been knocked down, respectively. In the rescue group, we observed increased cell migration and proliferation, as evidenced by transwell, wound healing, CCK-8, and EdU assays, following the manipulation of DTX3L and BST2 (Additional file 2: Fig. S12C-F, H–K). Moreover, the expression levels of M2 macrophage markers in the rescue group were similar with those of the control group (Additional file 2: Fig. S12G, L).

Overall, the in vitro and in vivo data support the notion that DTX3L and BST2 participated in accelerating malignancy of ESCC cells, and tumor-associated M2 macrophages served as critical characteristics in ESCC progression.

Discussion

Determining ESCC progression facilitates a profound understanding of the disease. Previous scRNA-seq studies unraveled the unidentified cell types, states, and potential functions involved in ESCC progression. Yao et al. deciphered the cell transition states in a multistep ESCC progression model in mice and identified a set of pivotal transitional signatures associated with the oncogenic evolution of epithelial cells [65]. Chen et al. revealed that epithelial cells activated fibroblasts in the microenvironment, thereby promoting ESCC development [16]. Nevertheless, single-cell profiling through tissue dissociation engenders the loss of information concerning the spatial localization of cells and their intercellular communications. Cancer growth and progression unfold intricately as spatial processes, entailing the breakdown of normal tissue organization, invasion, and metastasis. Tumor cells exist within an ecosystem comprising tumor cells and their surrounding TME. Therefore, dissecting and understanding this intricate ecosystem is crucial for unraveling the mechanisms driving tumor progression and developing effective novel treatments.

In this study, multi-region transcriptomic profiling was conducted using DSP to comprehensively characterize the cancer (immune) microenvironment throughout ESCC progression. DSP utilization represents a novel avenue for investigating cancer evolution, allowing comprehensive elucidation of molecular and spatial details and alterations occurring during the early stages of tumorigenesis, which is unattainable through bulk or single-cell transcriptomic studies alone. Herein, the distinct processes and microenvironment compositions linked to ESCC progression were elucidated, and potential biomarkers and therapeutic targets were identified. These findings were orthogonally validated using single-cell spatial transcriptomics, IHC, and publicly available bulk RNA-seq and scRNA-seq datasets. Functional studies performed on DTX3L and BST2 revealed their inhibitory effects on ESCC cell proliferation and migration, as well as on M2 macrophage polarization.

In addition to re-capitulating the well-established biological pathways recognized as the hallmarks of cancer [10], our approach provides a deeper understanding of the specific biological compartments in which these alterations occur. Changes in oncogenic, metabolic, and immune-related processes throughout ESCC tumorigenesis were delineated. Substantial shifts occurred in the cellular composition from normal tissue to cancerous lesions within the stromal segment. Stromal components within the TME play critical roles in tumor initiation, progression, and metastasis [66]. The current study sheds light on the extensive remodeling of the TME during ESCC progression. Increased proportions of immunosuppressive Tregs and macrophages were detected, implying the establishment of a progressively immunosuppressive microenvironment that favors tumor cell adaptation. Tumor-infiltrating macrophages within the TME are the key drivers of cancer progression and support tumor growth [42,43,44]. The cell communication network analysis revealed enhanced crosstalk between tumor epithelial cells and macrophages, mainly through MIF-(CD74 + CXCR4) interactions. The interaction between CXCR4 and CD74, which leads to MIF binding, triggers downstream signaling pathways implicated in tumor progression [67, 68]. Strategies aimed at inhibiting these pathways may offer novel therapeutic interventions for ESCC.

TAMs can be broadly categorized into two phenotypes: the M1 phenotype with anti-tumor properties and the M2 phenotype with pro-tumor characteristics [69]. Consequently, unraveling the crucial factors governing M2 macrophage polarization is pivotal in suppressing TAM-mediated tumor progression [70]. This study unveiled a compelling correlation between DTX3L and BST2 expression and M2 macrophage polarization. The functional study results indicated that decreased DTX3L or BST2 expression inhibited polarization of TAMs toward the M2 phenotype. Thus, increased DTX3L or BST2 expression during ESCC progression may foster an immune-suppressive TME. Similarly, Liu et al. highlighted the role of BST2 in promoting M2 macrophage polarization in cervical cancer [71], while He et al. reported significant involvement of BST2 in CRC progression and M2 TAM polarization [63].

The present study establishes a robust framework for the identification of novel spatial biomarkers by examining the dynamics of gene expression from normal tissue to cancer in pT1 ESCC, contributing to early disease interception and enhanced clinical management. Population screening programs strive to achieve early ESCC diagnosis and treatment, thereby reducing patient morbidity and mortality [72]. DTX3L, a member of the Deltex family [73], is involved in tumor progression in various cancer types [39, 60, 74]. As an E3 ligase, DTX3L regulates ESCRT-0 ubiquitination [75]. DTX3L also influences the DNA damage response pathway by modulating the mono-ubiquitination of histone H4 [76]. To the best of our knowledge, the current study is the first to demonstrate that DTX3L is overexpressed in ESCC tissues and its silencing impedes ESCC cell growth and migration.

BST-2 is a cell surface protein overexpressed in several solid tumors [77, 78]. Its tumorigenic role has been defined in gastric cancer [79] and CRC [63]. In the current study, treatment of Eca-109 or TE-1 cells with the CM collected from THP-1-derived macrophages with BST2 knockdown or treated with anti-BST2 caused a stark decrease in the proliferative and migratory abilities of ESCC cells. Encouragingly, immunotherapy with an anti-BST2 antibody has shown promise in patients with multiple myeloma in phase I studies. The BST2 antigen is also a novel immunological target for lung cancer treatment using anti-BST2 antibodies [80]. Additionally, promising results have emerged regarding treating glioblastoma with BST2-chimeric antigen receptor T-cell therapy [81]. Hence, the identification of BST2 as a promising immunological target for ESCC warrants further investigation.

Recently, Liu et al. (2023) conducted spatial transcriptomics analysis using DSP to investigate molecular markers indicative of the potential progression of esophageal squamous precancerous lesions (ESPL) into ESCC. They highlighted TAGLN2 and CRNN as candidate indicators for ESCC risk at ESPL stages. While our study and that of Liu et al. focus on ESCC progression using DSP technology, notable distinctions exist. Primarily, in this study, DSP’s segmentation feature was employed to segment ROIs into molecularly defined tissue compartments, including the epithelial (PanCK + /CD45 −), immune cell (PanCK − /CD45 +), and non-immune cell stromal (PanCK − /CD45 − /SYTO13 +) compartments. This approach facilitated the investigation of molecular alterations within each compartment and the identification of candidate biomarkers exhibiting specific aberrations during ESCC progression. Conversely, Liu et al. concentrated solely on the whole tissue without compartmentalization. Second, a comprehensive exploration of immune-related alterations during ESCC progression was performed. Additionally, although DSP technology offers highly multiplexed profiling capabilities, precisely identifying cell boundaries and capturing transcripts at the single-cell level remains challenging. These limitations were overcome in our study by using cutting-edge CosMx SMI technology to elucidate the spatial distribution of RNA within tissues at a subcellular resolution.

This study has certain limitations. First, the potential impact of genomic alterations on reprogramming was not explored due to a lack of genomic data. Additionally, despite pT1 ESCC representing an ideal model for spatial methodologies as different histology exits within the same lesion, expression levels in one compartment may be influenced by the presence of another compartment.

Conclusions

Through the implementation of spatial transcriptomic analysis, the dynamic biological processes and cell atlas associated with ESCC progression were elucidated. Furthermore, a robust framework for the identification of consistently altered biomarkers throughout ESCC tumorigenesis was established. This facilitated the evaluation of the underlying molecular characteristics driving ESCC progression, thereby paving the way for novel therapeutic strategies. Additionally, the role of TAMs in ESCC tumorigenesis was characterized. The findings gained from our research offer essential knowledge for comprehending the pathological progression of ESCC development and hold the potential to act as an initial alarm for ESCC, thereby aiding in the prevention and timely intervention of ESCC.