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 µm2. 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% CO2 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 × 106 cells). Tumor growth was monitored over a period of 28 days. Tumor volume was calculated using the standard formula a2 × 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.