Deciphering the generation of heterogeneity in esophageal squamous cell carcinoma metastasis via single-cell multiomics analysis

Abstract

Background

Chromatin accessibility plays a crucial role in mediating transcriptional dysregulation and heterogeneity in Esophageal Squamous Cell Carcinoma (ESCC). Examining the chromatin accessibility of ESCC at single-cell level is imperative to understand how it activates oncogenes and contributes to the onset and metastasis of ESCC.

Methods

We performed single-cell assay for transposase-accessible chromatin sequencing (scATAC-seq) on cancerous and adjacent noncancerous tissues from four ESCC patients who were pathological staged as T1a, T2b, T3b, or T4a, to investigate whether regulatory elements are pivotal in instigating cellular heterogeneity during ESCC metastasis. In conjunction, we integrated these data with 55 scRNA-seq datasets, ChIP-seq or CUT&Tag sequencing data, Hi-C sequencing data, bulk RNA-seq data, and bulk ATAC-seq data from ESCC cell lines to dissect the mechanisms underlying the heterogeneity of ESCC and tumor microenvironment (TME).

Results

Our study identified enhancer-specific activation within epithelial cells orchestrated by the three-dimensional structure of chromatin that regulates SERPINH1 transcription, and promotes the epithelial-mesenchymal transition (EMT) and metastasis of ESCC. Additionally, chromatin element activation facilitated the expression of TNFSF4 in CD8 + exhausted T cells, thereby activating Tregs. Furthermore, we observed that chromatin accessibility promoted the differentiation of tumor-associated macrophages (TAMs) and cancer associated fibroblasts (CAFs).

Conclusions

In summary, utilizing multiomics analyses, we have revealed chromatin accessibility maps and illuminated the intricate molecular mechanisms that underlie cellular heterogeneity during ESCC metastasis, offering valuable insights to further advance research on tumor progression and deterioration.

Background

Esophageal Squamous Cell Carcinoma (ESCC) is a common malignant digestive tract tumor and is the sixth leading cause of cancer-related death, with a five-year survival rate of only 15–25% [1, 2]. Although endoscopic examination has become more popular in recent years, only a small number of patients with ESCC can achieve early diagnosis, and most individuals progress to advanced disease. Local lymphatic metastasis is the most common recurrence pattern and accounts for treatment failure. Most researches focused on the mechanisms of the progression and metastasis of ESCC [3, 4], but the activation process of key gene regulatory elements that promote the heterogeneous development and the tumor microenvironment (TME) of ESCC remains unclear. Therefore, a priority for oncologists is to explore ESCC heterogeneity and TME to achieve early diagnosis and prevent the metastasis of tumor.

ESCC exhibits significant tumor heterogeneity, especially evident in chromosome copy number, differential expression of cellular markers and functional disparities among its constituent cells throughout the occurrence and progression across various clinical stages [5,6,7,8]. This heterogeneity is initially observable in the peritumoral environment. Chen et al. have reported that some antitumor and inflammatory molecules, such as KRT13 and MUC21 were high expressed in non-malignant esophageal epithelial cells, whereas SERPINH1, a gene related epithelial-mesenchymal transition (EMT) was overexpressed in cancer cells [6, 9]. This result demonstrated the heterogenous gene expression between cancerous and noncancerous cells in ESCC tissues. In addition, the heterogeneity may not only exist between cancerous cells and noncancerous cells, the tumor cells even in one individual also exhibit a diverse array of cell subsets due to the continuous proliferation and distinct differentiation processes. These cell subsets exhibit unique characteristics related to immunology, cell cycle regulation, EMT-like properties and cell repair mechanisms, which further contribute to the overall heterogeneity of ESCC [6, 6,10,11,12]. Furthermore, the heterogeneity between metastatic and non-metastatic ESCC tissues also obviously exists. Yin, X et al. have revealed that the pathways related to keratosis, epidermal cell differentiation, and T-cell activation are significantly activated in ESCC without metastasis, while the tumor tissues in ESCC with lymphatic metastasis exhibit extremely high levels of immunoregulatory pathways through cNMF analysis [5]. Just as mentioned above, most studies have focused on the phenomenon of ESCC heterogeneity, the underlying mechanisms, especially the epigenetic mechanism which account for the heterogenous changes remains unclear.

Epigenetic disorders drive aberrant transcriptional programs, leading to cancer initiation and progression [13]. Chromatin accessibility, accompanied by histone modifications, can regulate the activity of oncogenes without changing the DNA sequence. Technologies such as ATAC-seq and ChIP-seq (such as H3K27ac ChIP-seq) enable the identification of the activation of regulatory elements [14]. By combining bulk ATAC-seq with TCGA analyses, researchers have shown that increased chromatin accessibly at EGFRe1, an enhancer near the EGFR gene locus, can promote the overexpression of EGFR and mediate the proliferation and migration of ESCC cells [15]. Furthermore, regions around transcriptionally active genes (NELL2 and PRIM2), exhibit more open chromatin structures, whereas those associated with genes with low-expression (such as HPGD, KAT2B, and RRAGD) present reduced chromatin accessibility [16]. However, these findings are primarily based on analyses of cell line and tissues, with a limited understanding of single-cell chromatin accessibility and its impact on the ESCC TME.

Single-cell assay for transposase-accessible chromatin sequencing (scATAC-seq) is a single-cell genomics technique that identifies open regions of chromatin to reflect gene expression and cellular function. Compared with bulk ATAC-seq, scATAC-seq can reveal regions of gene transcriptional activity at the single-cell level and can more precisely reveal the epigenetic mechanisms underlying tumor heterogeneity and TME evolution. Integrated scATAC and scRNA-seq analyses have revealed gene transcriptional activation mechanisms in various cancers, including endometrial cancer [17], colorectal cancer [18], clear cell renal cell carcinoma [19], bladder cancer [20], and glioma [21]. These studies demonstrated that chromatin accessibility changes activated gene transcription, thereby promoting tumorigenesis and metastasis at the single-cell level. Regner et al. reported that the increased accessibility of the LAPTM4B promoter and enhancers can promote the overexpression of LAPTM4B, increase chemotherapy resistance, and even correlate with the poor prognosis of ovarian cancer patients using the combination of scATAC-seq with scRNA-seq [17]. However, the role of chromatin accessibility in regulating gene transcription at the single-cell level in ESCC has not been reported.

The TME comprises an environment that encompasses tumor cells, adjacent tissues, and immune cells, fostering the survival of tumor cells [22, 23]. Variations in the TME cell composition and function significantly affect cancer initiation and progression [24]. Several researchers have proposed that CD8 + exhausted T cells are more abundant in tumor cells than in adjacent tissues [25]. Malignant epithelial cells can promote the transformation of normal fibroblasts (NFs) to cancer-associated fibroblasts (CAFs) by downregulating ANXA1, with CAFs contributing to the formation and migration of ESCC cells [26]. In addition, fibroblasts are significantly enriched in cancer regions, whereas immune cells (T cells, NK cells, and B cells) are more abundant in stromal regions in ESCC [8]. However, the epigenetic mechanisms underlying these TME changes remain unclear. Therefore, a deeper understanding of ESCC TME dynamics and their underlying epigenetic regulation is crucial for developing effective therapeutic strategies for ESCC.

In this study, we performed scATAC-seq on the tumors and adjacent tissues of 4 patients with ESCC at different disease stages, coupled with a comprehensive analysis of 55 groups of scRNA-seq datasets. The results revealed the epigenetic alterations in ESCC-related differentially expressed genes (DEGs). Furthermore, we explored the chromatin accessibility dynamics of ESCC at different stages and inferred the activity of transcription factors (TFs) that regulate the interaction of stage-specific elements. In addition, chromatin accessibility and transcription maps of CAFs, tumor associated macrophages (TAMs) and T cells in the TME were generated, providing a basis for the development of novel ESCC therapeutic strategies.

Methods

Patient selection and tissue collection

We retrospectively selected patients with ESCC who had undergone curative resection at Shandong Provincial Hospital Affiliated to Shandong First Medical University. These patients were pathologically staged as I-IV according to the 8th Edition of TNM Staging Manual published by the American Joint Committee on Cancer. From each stage, we selected 4 patients, resulting in a total of 16 patients included in our study. For each patient, a pair of ESCC and non-cancerous esophageal tissue (located more than 5 cm from the margin) were collected. To comprehensively analyze chromatin accessibility and gene expression, we employed a multi-faceted approach. 4 pairs of representative tissues (one per stage) were subjected to scACAC-seq (Supplementary Table S1). The remaining 12 ESCC tissues were analyzed using bulk ATAC-seq followed by qPCR to assess chromatin accessibility changes; Concurrently, IHC staining was performed on these 12 samples to evaluate the expression levels of the genes of interest (Supplementary Table S1). This study was approved by the Ethics Committee of Biomedical Research of Shandong Provincial Hospital Affiliated to Shandong First Medical University (NO.2022–160, March 1, 2022).

Cells culture

Four established human ESCC cell lines (KYSE30, KYSE150, KYSE510 and KYSE180) were purchased from the ATCC cell bank, and one normal human esophageal epithelial cell line (HEEC; BNCC359279) was a gift from Dr. Zhaohua Xiao (The Second Hospital Affiliated to Shandong University). All the cells were cultured in Roswell Park Memorial Institute medium 1640 (RPMI 1640). The media were supplemented with 10% fetal bovine serum (FBS), 100 U/mL penicillin, and 100 U/mL streptomycin (all of these reagents were obtained from Gibco, Invitrogen, Inc, USA). All the cell lines were maintained in a humidified atmosphere with 5% CO₂ at 37℃.

Plasmids and cell transfection

SERPINH1 overexpression plasmids were constructed by cloning the human SERPINH1 gene into the pcDNA3.1-3xFlag-C vector (pcDNA3.1-SERPINH1-3xFlag). The cells were transfected with pcDNA3.1-SERPINH1-3xFlag (OE_SERPINH1) or pcDNA3.1-3xFlag-C (Empty Vector) using Lipofectamine 3000 (Invitrogen, USA).

scATAC library construction and sequencing

Approximately 0.1 g of tissue was obtained, ground, digested with a micropestle, and passed through a 40-micron cell strainer, after which all the cells were collected by resuspension in precooled PBS. The cells were pelleted via centrifugation at 500 ×g for 5–10 min at 4 °C. After washes with PBS and centrifugation, the cells were resuspended in 50 μl of cold nuclei lysis buffer and incubated on ice. Nuclei were isolated by centrifugation at 500 ×g for 10 min at 4 °C. Following Trypan blue staining and cell counting the nuclei concentration was adjusted to the desired capture number. scATAC libraries were prepared using a 10 × Genomics single-cell ATAC kit and sequenced on the Illumina NextSeq 500 platform according to the manufacturer's recommendations.

scATAC-seq data analysis

The raw reads from each patient sample were aligned to the hg38 reference genome, and a list of unique ATAC-seq fragments with associated barcodes was generated using 10 × Genomics Cell Ranger ATAC. The list of unique fragments per barcode for each patient’s sample was read into the R package ArchR (version 1.0.2) [27]. Cells with > 1000 unique nuclear fragments and > 4 TSSs were preserved in accordance with the strict quality control (QC) criteria for scATAC-seq data. Doublet removal was performed via “addDoubletScores” and “filterDoublets” functions within ArchR. The filterRatio was configured to be 1 during this process. Then, further quality control of the sequencing depth of each dataset was performed according to the TSS enrichment and log10(nFrags) criteria [27].

For dimensionality reduction of the scATAC-seq data, iterative latent semantic indexing (LSI) in ArchR [27] was performed using the "addIterativeLSI" function, with 25,000 variable features and 30 dimensions. A batch effect correction tool called “Harmony” was applied to remove strong batch effect differences. Subsequently, clusters were identified by utilizing the 'addClusters' function with a resolution parameter set at 0.8. Thereafter, a two-dimensional visualization of the data was generated through the employment of the 'addUMAP' function, with the parameters configured as nNeighbors = 30, minDist = 0.5, and metric = cosine, to enhance the clarity and interpretability of the data representation. Pseudo bulk replicates were generated using the “addGroupCoverages” function to obtain measurements of statistical significance. Peak calling was conducted using the “findMacs2” function, followed by reproducible peak set creation and peak matrix generation with the “addReproduciblePeakSet” and “addPeakMatrix” function. Furthermore, the “getMarkerFeatures” function was used to identify marker peaks based on TSS enrichment and the log10(nFrags) value, and a t-test for statistical significance was performed. Coaccessibility was used to determine correlations in accessibility between two peaks, and peak-to-gene linkage leakage was used to integrate scRNA-seq data and determine correlations between peak accessibility and gene expression [27].

scRNA-seq data analysis

scRNA-seq data were obtained from the GSE160269 dataset of GEO (https://www.ncbi.nlm.nih.gov/geo/) and the PRJNA777911 and PRJNA971344 (in house) datasets of SRA (https://www.ncbi.nlm.nih.gov/sra) (Supplementary Table S2). The detailed clinical information is summarized in Supplementary Table S2. The raw reads of the latter two datasets were aligned to the hg38 reference genome and were counted using 10 × Genomics Cell Ranger with the default parameters. In addition, Seurat (version 4.2.0) [28] was used to convert the matrix to a Seurat object and for QC. Samples whose number of cells was less than 1000 were removed. To eliminate potentially low-quality cells, those with a high proportion of mitochondrial gene expression were excluded. Additionally, genes that were expressed in fewer than three cells were also removed from the analysis. Furthermore, the data were integrated via CCA to eliminate the influence of data sources on the data. Meanwhile, the 'FindVariableFeatures' function was employed to detect the top 2000 highly variable genes. Then, Principal component analysis (PCA) was subsequently performed with the “RunPCA” function in Seurat. The “FindNeighbors” and “FindCluster” functions in Seurat were used to sort the cells. The clustering results were visualized in UMAP plot generated with the “RunUMAP” function with the top ten principal components [29].

A cell annotation approach based on well-characterized gene markers was adopted, as detailed in a previous publication [9]. Briefly, EPCAM, KRT19, and KRT5 in epithelial cells; FN1, COL1A1, and COL3A1 in fibroblasts; VWF, PECAM1, and SELE in endothelial cells; CD3A and CD3D in T cells; CD79A and CD79B in B cells; CD163 in macrophages; TPSAB1 in mast cells; and CLEC9A in dendritic cells were analyzed.

scATAC integration with scRNA-seq and scATAC-seq annotations

Gene activity scores were inferred via scATAC-seq using the “addGeneScoreMatrix” function from ArchR before the marker was transferred from the scRNA-seq data to the scATAC-seq data. The scATAC-seq gene scoring matrix (GeneScoreMatrix) was compared with the scRNA-seq gene expression matrix (RNAAssay) using the “addGeneIntegrationMatrix” function. The scATAC-seq data were subsequently integrated with the scRNA-seq data. This integration allowed us to transfer cell type annotations from the well-characterized scRNA-seq data to the scATAC-seq data. Briefly, cells from scATAC-seq data were directly aligned with corresponding cells from the scRNA-seq data by comparing the scATAC-seq gene score matrix with the scRNA-seq gene expression matrix, using the “FindTransferAnchors” function from the Seurat R package and the “addGeneIntegrationMatrix” in the ArchR R package.

Copy number variation (CNV) analysis

The CNV assessment of epithelial cells was performed with the infercnv (version 1.10.1) R package [30]. Firstly, the RNA microarray matrices were extracted from scRNA-seq data using the “GetAssayData” function and from scATAC-seq data using the “GeneScoreMatrix” function, with the average gene copy number of T cells used as a reference. Subsequently, AnnoProbe (version 0.1.7, GitHub—jmzeng1314/AnnoProbe) was used to construct the gene map reference files. Finally, by comparing the gene copy numbers and references of epithelial cells, cells with CNV scores above 1000 were considered malignant. This threshold was chosen based on a reference standard, where a CNV score above 1000 is typically associated with genomic instability, a characteristic commonly observed in malignant cells.

Pathway enrichment analysis

First, “FindAllMarkers” from Seurat was used to calculate tumor and normal DEGs via scRNA-seq. Moreover, "getMarkerFeatures" from ArchR were applied to calculate DEGs in scATAC data according to TSS enrichment and log10 (nFrags) based on GeneScoreMatrix. Three pathway enrichment analyses were applied, and the details are provided below.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis [31] was performed using the Visualization and Ensemble Discovery Database (DAVID, https://david.ncifcrf.gov/home.jsp) [32].

For Gene Set Enrichment Analysis (GSEA), clusterProfiler (version 4.2.2) [33] was used to compare the differentially expressed genes obtained from RNA-seq data from SERPINH1 overexpressing and control cells with the gene set of the MSigDB (https://www.gsea-msigdb.org/gsea/msigdb) hallmark gene sets, after which the GSEA gene set was scored according to the alignment.

For Gene Set Variation Analysis (GSVA), after the samples were grouped, the gene set of the Homo sapiens “H” biological pathway was selected, and pathway enrichment scores for each sample were calculated using GSVA (version = 1.42.0) [34]. The matrix employed in this process was the normalized expression matrix, and the parameter of kcdf was set to Gaussian. Pheatmap (version = 1.0.12, https://cran.r-project.org/web/packages/pheatmap/index.html) was applied to visualize the pathway enrichment score based on the group information in a heatmap.

SERPINH1 regulation and enhancer prediction

Hi-C contact matrices were generated and normalized by HiCPro (version 2.11.1) [35] at 10 kb and 40 kb, and then chromatin loops with contact counts > 10 were identified by Fit-Hi-C (version = 2.0.8) [36] based on contact matrices at a resolution of 10 KB. Interaction heatmaps around SERPINH1 were visualized via Juicebox (version = 2.16.00) at 10 kb and 25 kb. Integrated Genomics Viewer (IGV) (version = 2.12.3) [37] was employed to visualize scATAC-seq, ATAC-seq, ChIP-seq and RNA-seq data. Loops were visualized with the R package Sushi (version = 1.10.0) [38], with contact counts marked on the vertical coordinates. See Supplementary Table S3 for specific information [39,40,41,42].

CUT&Tag and ChIP-Seq data analysis

The sequencing data were mapped to hg38 using Bowtie2 (version = 2.2.5) [43] and then sorted and indexed with SAMtools (version = 1.6) [44]. MACS2 (version = 2.2.7.1) [45] was subsequently employed to perform peak calling. The peak annotations were obtained and displayed using the ChIPseeker R package (version = 1.30.3) [46]. The SAMtools processing output files were transformed into bw format using the bamCoverage tool from deepTools (version = 3.5.1) software [47], and the enrichment of peaks near specific regions such as the TSS was assessed using computeMatrix from deepTools software. IGV [48] was employed for data visualization.

Analysis of bulk RNA-seq data

The sequencing data were assessed in FASTQC (version 0.11.9, https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) format, and based on the quality evaluation, Trimmomatic (version 0.39) [49] was used to filter the data, remove adapters and remove low-quality reads. The filtered result files were subsequently aligned to the hg38 reference genome using HISAT2 (version 2.2.1) [50], and HTSeq (version 2.0.2) [51] was then utilized to quantify the gene expression levels in each sample. Afterward, the R package DEseq2 (version 1.34.0) [52] was subsequently employed to identify and statistically analyze the DEGs. Genes with an absolute log2FoldChange ≥ 1.5 and a P value < 0.05 were considered significantly differentially expressed.

Analysis of transcription factors

For the analysis of scATAC-seq data, the peaks were annotated with "addMotifAnnotations" and the motifs were sorted in each cell type with "peakAnnoEnrichment" from ArchR [27]. For the scRNA-seq data, the python package pySCENIC (version 0.12.0) [31, 53] was used based on the hg19-tss-centered-10kb-7species database to evaluate the activities of TFs [31]. The input matrix was the normalized expression matrix from Seurat. Coexpression modules were inferred with GRNboost (https://github.com/aertslab/GRNBoost), pruned modules for targets with RcisTarget (version 1.24.0, https://github.com/aertslab/RcisTarget), and enriched cellular regulons with AUCell (version 1.22.0) [31]. Finally, the top 50 motifs were intersected in the scATAC-seq data with the upregulated motif in the scRNA-seq and NetworkAnalyst (https://www.networkanalyst.ca) was used to map the TF interaction network.

Trajectory analysis

ArchR was employed to create cellular trajectories that order cells across a lower N-dimensional subspace in pseudotime within an ArchRProject of the scATAC-seq data [27]. The “addTrajectory” function was used to create a trajectory and incorporate this trajectory into the data, which were visualized using the “plotTrajectory” function. Moreover, Monocle3 (version 1.3.4) was used to perform the trajectory analysis to explore the pseudotime trajectory of the scRNA-seq data [54]. The “estimateSizeFactors” and “estimateDispersions” functions were used to construct statistical models to characterize the data. In addition, the “reduceDimension” and “orderCells” functions were used to reduce the dimensions of the cells and to place the cells on a pseudotime trajectory.

Cellchat analysis

The R package CellChat (version 1.6.1) was used to explore cell–cell interactions among epithelial cells, CAFs, TAMs and T cells. Moreover, CellChat was utilized to identify tumor-specific pathways and ligand-receptor interactions between different cell subsets [55, 56].

Cell migration assay

Both the SERPINH1-overexpressing KYSE150 cells (OE_SERPINH1) and the control cells (Vector) were precultured in serum-free medium for 24 h. Then, the starved cells were resuspended in serum-free culture medium, and 5 × 10⁵ cells were seeded into the upper chambers of 24-well transwells (8-μm pore size polycarbonate membrane; Merck Millipore Bioscience, Germany) with 200 μL of serum-free medium. Six hundred microliters of basal culture medium containing 20% FBS was added to the lower chamber. The plate was then placed in a cell culture incubator for 24 h. Afterward, the medium in the transwell was discarded, and the cells were washed with PBS three times. Methanol was used to fix the cells for 30 min, followed by three additional washes with PBS. The cells were subsequently stained with 0.1% crystal violet for 30 min. The upper layer of cells was removed using a cotton swab, and the cells were washed three times before being observed under a microscope.

Immunohistochemistry (IHC)

The streptavidin-peroxidase method was used to detect DPP9 and TPM1 expression in the samples. Formalin-fixed and paraffin-embedded surgical samples were sequentially cut into 4 μm sections. The sections were then dewaxed, and antigen retrieval and hydrogen peroxide incubation were performed in sequence. Rabbit anti-DPP9 (Abcam, Cambridge, UK) and rabbit anti-TPM1 antibodies (Proteintech, IL, USA) were used at dilutions of 1:100 and incubated with the sections at 4 °C overnight. The subsequent experimental steps were performed using a secondary biotinylated antibody kit purchased from ZSGB Biotech (Beijing, China) according to the manufacturer’s instructions.

Bulk ATAC library construction and real-time qPCR quantification

Nuclei were extracted from tissues using the protocol detailed in the “scATAC Library Construction and Sequencing” section. Briefly, a 0.1 g tissue sample was ground and filtered through a 40-micron cell strainer, followed by washes with cooled PBS. Cold cracking buffer was added to the pellet, and the mixture was agitated, and centrifuged for 10 min at 4 °C. The supernatant was removed, and the pellet was chilled on ice. Nuclei were counted after staining with 0.04% Trypan blue. The ATAC library was generated using the Hyperactive ATAC-Seq Library Prep Kit for Illumina (Vazyme, China) according to the user manual with minor adjustments. The isolated nuclei were subjected to chromatin fragmentation for 30 min at 37 °C with fragmentation buffer containing Tn5 transposase and PCR adapters, during which PCR adapters were ligated to fragmented DNA ends. The reaction was terminated by an incubation with stop buffer for 5 min. The resulting DNA fragments were purified using ATAC DNA extraction beads and amplified by PCR. Finally, library size selection was performed with 0.55 × and 1.2 × volume ATAC DNA Clean Beads. The purified ATAC-Seq library was enriched for open chromatin fragments. Their relative abundance was quantified by real-time qPCR using the library as a template. The following real-time qPCR primers were used:

Target	Forward sequence (5′-3′)	Reverse sequence (5′-3′)
GAPDH	GAGCAGTCCGGTGTCACTAC	TCGAACAGGAGGAGCAGAGA
SERPINH1	TGTGCAGATGATCCAAGGGG	CACACCCAACCCCATCTCTG
TPM1	CAGGGCTCCATGGTATTGCA	CTGTTGAGAGGCTCACGCAT
DPP9	ATCAGGGCGTTGGTGGTATAG	AACTCGCAAAAGTGGTACCG

Results

Single-cell chromatin accessibility and the transcriptome landscape in ESCC

scATAC and scRNA sequencing were performed on ESCC tumors at different clinical stages to study the effect of chromatin accessibility on the occurrence and development of ESCC (Fig. 1a). In our previous study, we revealed heterogeneity in ESCC using scRNA-seq data from many clinical samples (Fig. 1b) [9]. In this study, we performed scATAC-seq on 8 groups of tumoral and peritumoral tissue samples from 4 patients with ESCC at various TNM stages (pT1N0M0 (stage IA), pT3N0M0 (stage IIB), pT3N1M0 (stage IIIB), and pT4N2M0 (stage IVA)) (Supplementary Table S1) and constructed a chromatin accessibility map for ESCC and explored its association with gene transcription to investigate the relationship between chromatin accessibility and ESCC heterogeneity (Fig. 1a). After quality control and dimensionality reduction of the scATAC-seq data and annotation with reference to the scRNA-seq cell population, we ultimately identified 8 cell populations within the scATAC-seq data: epithelial cells (n = 12,061), fibroblasts (n = 20,111), macrophages (n = 5008), T cells (n = 9472), endothelial cells (n = 2252), B cells (n = 4331), dendritic cells (DCs, n = 38), and mast cells (n = 2046) (Fig. 1b). Among these cell populations, epithelial cells and fibroblasts were the two most abundant cell populations, while T cells proportions varied significantly across samples (Fig. 1c). Moreover, we annotated the scATAC-seq peak (Fig. 1d) and identified markers of different cell types: EPCAM, KRT19, and KRT5 for epithelial cells; FN1, COL1A1, and COL3A1 for fibroblasts; VWF, PECAM1, and SELE for endothelial cells; CD3A and CD3D for T cells; CD79A and CD79B for B cells; CD163 for macrophages; TPSAB1 for mast cells; and CLEC9A for dendritic cells (Fig. 1e, f). As these genes are well-established markers for different cell types, the above results indicate that the degree of chromatin accessibility is closely related to cellular identity.

Fig. 1

Single-cell chromatin accessibility and the transcriptome landscape in ESCC. a Schematic overview of scATAC and scRNA-seq ESCC sampling and subsequent analysis. Tumor and adjacent nontumor tissue samples were collected from patients with esophageal squamous cell carcinoma (ESCC) and subjected to scATAC-seq and scRNA-seq. The annotations from scRNA-seq were mapped onto scATAC-seq cells, and a combined analysis of gene expression levels and chromatin accessibility within the same cell subpopulations was conducted to explore the impact of chromatin accessibility on gene expression levels. b UMAP plot of all the scRNA-seq cells colored by cell type (left); UMAP plot of all the scATAC-seq cells colored by inferred cell type (right). c Stacked bar charts showing the contribution of each sample to each cell type according to scATAC-seq. d Row-scaled heatmaps of statistically significant distal peak-to-gene links. Each row represents the expression of a gene (right) correlated with the accessibility of a distal peak (left). e UMAP plots illustrating key marker genes for identifying cell types in scATAC, including epithelial cells (KRT5), fibroblasts (FN1), endothelial cells (VWF), T cells (CD3D), B cells (CD79A), and macrophages (CD163). f Heatmap of Z scores for significant differential elements in scATAC-seq clusters. Single-cell cluster identities are indicated at the top of the plot to identify the different clusters

Heterogeneity in the chromatin accessibility of ESCC epithelial cells

Given the predominance of epithelial cells, we subdivided them into five distinct subpopulations (Fig. 2a). The InferCNV analysis revealed that one of the subsets contained more nonmalignant epithelial cells than malignant cells; thus, we named these cells “epithelial cells”. The remaining subsets contained more malignant cells; thus, we termed these cells “tumor cells 1–4” (Fig. 2a, b). We integrated scRNA-seq and scATAC-seq data to explore the atypia between malignant tumor cells and nonmalignant normal epithelial cells in ESCC. Compared with their normal counterparts, malignant tumor cells presented upregulation of 3709 peaks and downregulation of 6030 peaks, which in turn altered gene expression (Fig. 2c). We found that 471 DEGs were consistent between the scATAC and scRNA-seq data (Supplementary Fig. S1a), among which 243 genes, including CCND1, ANO1, and CTTN, were upregulated in ESCC compared with peritumor samples in both the scATAC and scRNA-seq data (Fig. 2d). Compared with normal esophageal tissues, increased chromatin accessibility at these gene loci correlated with increased gene expression in malignant epithelial cells (Fig. 2e, g; Supplementary Fig. S1b, c). We mapped the chromatin accessibility map of CCND1 to observe the chromatin openness of CCND1 in detail and observed a strong chromatin accessibility signal at the CCND1 promoter in malignant epithelial cells compared with that in nonmalignant epithelial cells (Fig. 2f), suggesting that increased chromatin accessibility at the CCND1 cis-regulatory element promotes its elevated expression in malignant tumor cells. Conversely, 149 genes, such as KRT13 and RHCG, were downregulated in nonmalignant epithelial cells according to both the scATAC and scRNA data (Fig. 2d). These genes exhibited greater chromatin accessibility, active transcription, and gene expression in adjacent nonmalignant epithelial cells than those in malignant tumor cells (Fig. 2e, g; Supplementary Fig. S1b, d).

Fig. 2

Heterogeneity in the chromatin accessibility of ESCC epithelial cells. A UMAP plot of all epithelial cells colored by the cell subtype according to the scATAC-seq data. B UMAP plot of all epithelial cells colored according to inferCNV type in the scATAC-seq data. C Volcano plots showing differential epithelial peaks in malignant and nonmalignant cells. A total of 3709 upregulated peaks and 6030 downregulated peaks were identified. The log2mean is the logarithm of the arithmetic mean of the absolute value of the deviation between genescore and the mean of genescore. D Scatter plots showing DEGs between scATAC and scRNA-seq data from malignant and nonmalignant epithelial tissues in ESCC. Genes with log2FC > 0.25 in both the scATAC-seq and the scRNA-seq data were considered upregulated genes; those whose log2FC < -0.25 in both the scATAC-seq and the scRNA-seq data were considered downregulated genes. E The UMAP plots are colored according to log-normalized gene scores, indicating the accessibility of cis-regulatory elements linked to the indicated gene (CCND1 and KRT13) from the scATAC-seq data (upper). The UMAP plots show all epithelial cells colored according to gene expression (CCND1 and KRT13) from the scRNA-seq data (lower). F Browser track showing the accessibility profile at the CCND1 locus across malignant and nonmalignant cells. G UMAP plots of all epithelial cells colored by inferCNV type based on the scRNA-seq dataset. H KEGG signaling pathway diagram of DEGs associated with regulatory elements. The horizontal axis represents -log10(P-value), and the color indicates the number of genes involved. I Vein plot showing the intersection of DEGs in scATAC-seq (Log2FC > 0.25), scRNA-seq (Log2FC > 0.25), RNA-seq (KYSE150 vs. HEEC, Log2FC > 1), and RNA-seq (ESCC vs. peritumor clinical samples, Log2FC > 1) data from malignant cells and nonmalignant cells (left). Browser track showing the accessibility profile at the SERPINH1 locus across malignant and nonmalignant cells. The putative cancer-specific dRE of SERPINH1 is highlighted by a light yellow shadow (middle). Matching of the scRNA-seq expression of SERPINH1 is shown for all subclusters (right). J The UMAP plot is colored according to log-normalized gene scores, showing the accessibility of cis-regulatory elements linked to SERPINH1 based on the scATAC-seq GeneScoreMatrix

We performed a KEGG pathway enrichment analysis on the DEGs identified above to clarify the functions of genes regulated by cis-regulated elements in malignant tumors and normal cells. In malignant tumor cells, genes related to the hippo signaling pathway, focal adhesion, proteoglycans in cancer and other cancer-related signaling pathways were enriched (Fig. 2h). Notably, CCND1 is a key gene in these pathways, and CTTN is involved in the function of proteoglycans in cancer. Conversely, genes associated with decreased chromatin openness, such as KRT13, were enriched primarily in the p53 pathway, axon guidance, the estrogen signaling pathway and other signaling pathways (Fig. 2h).

We intersected the above DEGs with the DEGs from the KYSE150 and HEEC cell lines, as well as DEGs from clinical bulk RNA-seq data, to further identify the key genes that promote the progression of ESCC. We identified 6 genes that were specifically expressed in tumors, including ALCAM, FAT1. FSCN1, SERPINH1, ECT2 and HAS3, in all the datasets listed above (Fig. 2i, Supplementary Fig. S1e-i). Gene Ontology (GO) functional enrichment analysis of the above six DEGs revealed that they were closely related to cell migration (Supplementary Fig. S1j).

Based on our previous finding that SERPINH1 is a pivotal gene involved in the EMT [9], we hypothesized that high chromatin accessibility within in SERPINH1 region might promote gene transcription and facilitate ESCC metastasis. A chromatin accessibility map of SERPINH1 revealed that cis-regulatory elements proximal to the SERPINH1 locus were more specifically activated in ESCC tissue than in peritumor tissue (Fig. 2i). Intriguingly, despite the specific overexpression of SERPINH1 in tumor tissues, we observed notable variations in the chromatin accessibility and transcript levels of SERPINH1 across different epithelial cell subtypes. Notably, while the SERPINH1 promoter region presented the highest degree of chromatin accessibility in the tumor cells 2 subset, the tumor cells 1 subset presented the highest level of SERPINH1 transcription. This difference was attributed primarily to the more open chromatin state of the five cis-regulatory elements in the tumor cells 1 subset than in the other subsets, which facilitated gene transcription (Fig. 2i). Notably, SERPINH1 is predominantly highly expressed in early stage (I and II) ESCC (Fig. 2i, j), consistent with the findings of our previous study [9], and indicates that SERPINH1 plays an early role in activating the EMT signaling pathway. This comprehensive understanding of SERPINH1's role in ESCC pathogenesis offers valuable insights for the development of targeted therapeutic strategies.

Enhancers regulate the expression of SERPINH1 through the three-dimensional structure of chromatin

Topologically associated domains (TADs) are highly conserved in evolution, and rearrangements of TAD boundaries may lead to dysregulated of gene expression by altering the interactions between genes and enhancers. We compared the three-dimensional chromatin structures surrounding the SERPINH1 locus in KYSE150 and HEEC cells at a resolution of 25 kb and found that the intra-TAD interactions in KYSE150 cells were significantly greater than those in HEEC cells (Fig. 3a). At a resolution of 10 kb, we observed that the intensity of chromatin loop interactions in KYSE150 cells was twice as high as that in HEEC cells (Fig. 3a, c). Through an integrated analysis of single-cell ATAC-seq data from clinical samples, and ChIP-seq, CUT&Tag, ATAC-seq and bulk RNA-seq data from different ESCC cell lines, we found that the regulatory elements near SERPINH1 were significantly increased in ESCC cells compared to normal esophageal cells. Focusing on H3K27ac (transcription initiation marker) and H3K4me1 (active enhancer marker) histone modifications, we observed increased chromatin accessibility and H3K27ac enrichment at the Chr11:75,561,998–75,562,390 promoter region of SERPINH1 in all the ESCC cell lines and tissues compared with the normal controls (Fig. 3b). In the tissue samples, 10 regulatory elements were activated. By integrating H3K27ac and H3K4me1 data from cell lines, we identified Enhancer 1 (Chr11:75,492,270–75492469), Enhancer 2 (Chr11:75,499,416–75,499,784), and Enhancer 3 (Chr11:75,537,149–75,537,568) as active enhancers in both tissues and the KYSE30 and KYSE150 cell lines, as highlighted by the first three pink shaded regions in Fig. 3b. These enhancer regions were enriched with H3K4me1 modifications, indicating their potential to enhance gene transcription. Notably, these enhancer regions physically interacted with the H3K27ac-modified promoter region, collaboratively promoting the transcriptional initiation of SERPINH1. However, the enhancers were not activated in KYSE70, KYSE510, KYSE200, or TE5 cells, which may be attributed to the heterogeneity among different cell lines. We constructed bulk ATAC libraries for the KYSE150, KYSE510 and KYSE180 cell lines to further validate the chromatin accessibility of the SERPINH1 enhancer in different ESCC cell lines. We designed qPCR primers targeting Enhancer 2 to quantify the relative abundance of open chromatin within this region. Our results demonstrated that the accessibility of the SERPINH1 enhancer in KYSE150 cells was significantly greater than that in KYSE180 and KYSE510 cells (Fig. 3d). This increased accessibility resulted in significantly higher levels of SERPINH1 expression in KYSE150 than in KYSE180 and KYSE510 cells (Fig. 3b), suggesting that the accessibility of Enhancer 2 plays a key role in SERPINH1 transcriptional regulation. Additionally, our examination of clinical samples revealed greater chromatin openness around the SERPINH1 locus in early-stage ESCC samples than in late-stage samples (Fig. 3e). These data suggest that the decreased transcription level of SERPINH1 in the late clinical stage of ESCC may be attributed to the reduced enhancer accessibility within the SERPINH1 locus. Notably, the regulatory elements we identified exhibited significant chromatin interactions with the SERPINH1 promoter (Fig. 3c). These interactions may be achieved through loops and bridging structures in the three-dimensional chromatin structure, ensuring effective communication between enhancers and promoters. Therefore, we identified cancer-specific regulatory elements activated near SERPINH1, that interact with the SERPINH1 promoter through chromatin interactions to jointly regulate the expression of SERPINH1 in ESCC.

Fig. 3

Enhancers regulate the expression of SERPINH1 through the three-dimensional structure of chromatin. a Heatmap of contact domains around SERPINH1 based on hg38. The left shows the difference in chromatin interactions between HEEC (top) and KYSE150 (bottom) cells at 25 kb resolution (red: interactions with high frequency; blue: interactions with low frequency). The contact domain is marked by a blue rectangle. The right panel shows a magnified view of the domains around SERPINH1, with differences in the interactions highlighted in green circles. b Integrated Genomics Viewer (IGV) visualization of the differences in chromatin accessibility, histone modifications and RNA expression between normal and tumor cells and tissues. The promoter of SERPINH1 and possible enhancers are highlighted by red rectangles. c Significant intrachromosomal loops around SERPINH1. Only loops at 10 kb resolution with contact counts greater than 10 are shown. The vertical coordinates represent the interaction frequency, which varies greatly between tumor and normal cell lines. d Real-time qPCR quantification of SERPINH1 Enhancer 2 abundance in bulk ATAC libraries from the KYSE150, KYSE510 and KYSE180 cell lines (n = 3). e Real-time qPCR quantification of the SERPINH1 Enhancer 2 abundance in bulk ATAC libraries generated from ESCC samples at early (stages I and II) and advanced (stages III and IV) stages (n = 6). f, g GSEA showing the enrichment of EMT-related gene sets. h Transwell migration assays were used to assess the effect of SERPINH1 overexpression in KYSE150 cells (OE-SERPINH1) compared with control cells (Vector). i GSVA of RNA-seq data from the KYSE150 and KYSE30 cell lines, with OE and WT representing overexpression and control cells, respectively, and KD and NC representing the knockdown and negative control groups, respectively

Previous studies have showed that SERPINH1 is an early marker of the EMT in ESCC. We manipulated SERPINH1 expression levels through SERPINH1 knockdown and overexpression in the KYSE30 and KYSE150 cell lines, followed by bulk RNA sequencing and analysis to further explore the mechanism by which SERPINH1 is involved in the ESCC EMT and its role in the occurrence and development of ESCC. The functional enrichment analysis revealed that SERPINH1 overexpression (OE-SERPINH1) upregulated the expression of EMT-related genes, such as TIMP3, suggesting that SERPINH1 activated the EMT pathway (Fig. 3f). Consistent with these findings, the transwell cell migration assay confirmed that the migration ability of KYSE150 cells was significantly increased after the overexpression of SERPINH1 (Fig. 3h). Moreover, SERPINH1 overexpression activated multiple tumor-promoting pathways, including the p53 pathway, hedgehog signaling pathway, and inflammatory response pathway (Fig. 3i). Conversely, SERPINH1 knockdown (KD-SERPINH1) suppressed the expression of EMT-related genes, such as SEMA3C, ITGAV, and ITGA2, revealing an inhibitory state of the above pathways (Fig. 3g, i). In summary, our results indicated that the activation of SERPINH1 enhancers facilitates SERPINH1 transcription, not only regulates the EMT but also promotes tumor development and metastasis through multiple pathways, providing a new potential therapeutic target for ESCC.

Heterogeneity in ESCC epithelial cells at different stages

The effects of chromatin openness on ESCC heterogeneity and the stage-specific molecular characteristics were further investigated. When we examined epithelial cells according to the sample source, tumor cell subsets 1–4 represented malignant tumor cells from different patients with TNM1-4 tumors, respectively (Fig. 4a). This phenomenon persisted despite normalization and debatching, indicating significant differences in chromatin accessibility among ESCC tumor cells at different stages. We further explored this phenomenon from four aspects: gene copy number, chromatin openness, gene expression transcription, and TFs.

Fig. 4

Heterogeneity in ESCC cells at different stages. a UMAP plot of all epithelial cells colored by sample source and tumor stage according to the scATAC-seq data. b Heatmap visualizing differential peaks in epithelial clusters according to scATAC-seq data. c Heatmap illustrating large‐scale CNVs in epithelial cells (rows along the y‐axis) from tumor and peritumoral tissues. Red: amplification; blue: deletion. Epithelial cells from different cell subtypes and from different chromosome ranges are labeled with different color bars on the left and top of the heatmap, respectively. d Dot plot showing the expression levels of marker genes for epithelial cells from different stages for the scATAC-seq and scRNA-seq data. The dot size corresponds to the percentage of cells expressing the marker gene, and the dot color indicates the average expression level. e Browser track showing the accessibility profile at the TPM1 locus across malignant and nonmalignant cells. The putative cancer-specific dRE of TPM1 is highlighted in light yellow shadows. f Real-time qPCR quantification of TPM1 and DPP9 promoter abundance in bulk ATAC libraries generated from ESCC samples at various stages (n = 3). g Immunohistochemical staining results for the TPM1 (top) and DPP1 (bottom) proteins in stage I-IV ESCC tissues. h Scatter plot representing ATAC-seq chromVAR bias-corrected deviations in the 50 most highly expressed TFs across all epithelial cell types in the scATAC-seq data. Different colors indicate different single-cell cluster identifications. i Regulatory network of TFs identified by integrating the scATAC-seq and scRNA-seq data from patients with stage III ESCC. j Representative TF footprints of FOSL2 with motifs in the indicated scATAC-seq epithelial cells. The Tn5 insertion bias track is shown

We identified stage-specific peaks, with 3350, 2095, 2827, and 6652 peaks specifically upregulated in the tumor cells 1, tumor cells 2, tumor cells 3, and tumor cells 4 subsets, respectively (Fig. 4b). The copy number variation (CNV) analysis revealed more gene copies in tumor cells than in normal cells. However, significant differences in the degree of gene copy changes and chromosome numbers of gene variations in tumors of different stages were observed. The CNV density of TNM2-4 ESCC cells was significantly greater than that of TNM1 ESCC and normal esophageal epithelial cells, indicating that the CNV density of ESCC cells gradually increased with increasing tumor malignancy (Fig. 4c).

The genes expressed at different stages exhibited specificity at the cis-regulatory element activation and transcriptional levels. In normal esophageal epithelial cells, genes such as SPRR3 and RHCG are associated with increased chromatin openness and overexpression. MYL9 and TPM1 were activated in TNM1 ESCC, MT1X and CALM1 were activated in TNM2, LY6D and LY6K in TNM3, and UHRF1 and PDD9 in TNM4 (Fig. 4d). We mapped the chromatin openness of TPM1, a representative gene in the TNM1 stage to illustrate the differences in the chromatin accessibility of specific genes at different stages. The results showed that its four cis-regulatory elements were significantly activated in the TNM1 stage, driving the transcription of TMP1 and significantly upregulating its expression (Fig. 4e). In addition, real-time qPCR quantification of open chromatin in bulk ATAC libraries and IHC staining of ESCC tissues of different clinical stages revealed that TPM1 and DPP9 had higher chromatin openness and protein expression levels in TNM1 and TNM4, respectively, which was consistent with the results of the single-cell data analysis. (Fig. 4f, g).

We predicted motifs in different epithelial cell subsets and identified TFs that play key roles in different stages to identify the key TFs involved in the progression of ESCC (Fig. 4h). We found that the SRF, SRY and others was prominent in TNM1, the FOS family in TNM2, the KLF family in TNM3, and TP family in TNM4 (Fig. 4h). Moreover, we performed an analysis of the scRNA-seq data with pyscenic to identify TFs at the transcriptional level (Supplementary Fig. S2g). By intersecting the top 100 TFs predicted by scATAC with those upregulated in the pyscenic analysis, we obtained the key TFs in each tumor stage, constructed corresponding regulatory networks (Fig. 4i; Supplementary Fig. S2a, c, e), and blotted key TFs at different stages (Fig. 4j; Supplementary Fig. S2b, d, f). MEF2C, MAFF, FOSL2 and TP73 emerged as the core TFs for TNM stages 1–4, respectively (Fig. 4j; Supplementary Fig. S2b, d, f).

Chromatin accessibility influences cancer-associated fibroblast evolution

Cancer-associated fibroblasts (CAFs) are important components of the TME in ESCC. Previous studies have documented the progression of fibroblasts from NFs (normal activated fibroblasts (NAFs) / normal mucosa fibroblasts (NMFs)) to inflammatory cancer-associated fibroblasts (iCAFs) and then to myofibroblastic cancer-associated fibroblasts (myCAFs) in ESCC [26]. In accordance with our previous studies (Fig. 5a), we performed a quasi temporal analysis of fibroblasts, which is consistent with the identification of the evolutionary trajectory of ESCC (Fig. 5b) 9. However, whether this evolutionary trajectory is affected by chromatin accessibility is currently unclear. In this study, we annotated the fibroblast population within scATAC data with reference to the cell subsets of scRNA-seq data and identified six fibroblast subsets: NMF, NAF, iCAF1, iCAF3, myCAF1, and myCAF2 (Fig. 5c). Interestingly, a parallel evolutionary trajectory from NF (NAF, NMF) to iCAF and then to myCAF was observed at the chromatin accessibility level (Fig. 5d). A comparative analysis of the specific peaks revealed distinct patterns among the fibroblast subsets, with 2629, 7778 and 2137 peaks significantly upregulated in the NMF, iCAF1, and myCAF subsets, respectively. Moreover, we observed a progressive shift in peak profiles from NMF to iCAF to myCAF. For example, NAF exhibited upregulation of NMF peaks and a subset of iCAF peaks. Fifty-three peaks were specifically upregulated in the iCAF3 subset and were slightly upregulated in the iCAF1 and myCAF subsets (Fig. 5e).

Fig. 5

Chromatin accessibility influences cancer-associated fibroblast evolution. a UMAP plot of all fibroblasts colored by the cell subtype according to the scRNA-seq data. b Pseudotime analysis of fibroblast differentiation from NAFs, NMFs, and iCAFs to myCAFs using scRNA-seq data. c UMAP plot of all fibroblasts colored by the cell subtype according to the scATAC-seq data. d Pseudotime analysis of fibroblast differentiation from NAF, NMF, and iCAFs to myCAFs using scATAC-seq data. The line represents a double-spline fitted trajectory across pseudotime. e Heatmap visualizing differential peaks in fibroblast subtypes based on scATAC-seq data. f Pseudotime heatmap showing variable gene scores across fibroblast differentiation states based on scATAC-seq. g Pseudotime heatmap showing chromVAR TF motif bias-corrected deviations in the indicated lineage trajectory based on scATAC-seq. h Heatmap illustrating the dynamic changes in gene expression and GO pathway enrichment in fibroblast differentiation from NAF, NMF, and iCAFs to myCAFs based on scRNA-seq data. i Scatter plots showing differentially expressed genes between scATAC-seq and scRNA-seq data from fibroblasts in tumor and peritumoral tissues of ESCC. j, k Pseudotime representation of chromatin accessibility in scATAC-seq data (left) and gene expression in scRNA-seq data (right). The line represents a double-spline fitted trajectory across pseudotime

We hypothesized that the differentiation of fibroblasts is strongly correlated with the activation of regulatory elements during cancer development. Indeed, we observed dynamic changes in chromatin openness, TFs expression and gene expression during fibroblast differentiation (Fig. 5f-h). The results revealed that, in the NF stage, fibroblasts were mostly regulated by TFs such as TCF21 (Fig. 5g), and genes such as LAMA2 exhibited high chromatin accessibility (Fig. 5f, j). Consequently, LAMA2 was specifically transcribed, promoting ameboida-type cell migration and epithelial cell invasion (Fig. 5h). In the iCAF stage, fibroblasts were mostly regulated by TFs such as BATF (Fig. 5g), with genes such as KCNJ15 showing increased chromatin accessibility (Fig. 5f). The expression of these genes gradually increased, and the functions of RNA cleavage, respiration, and metabolism increased (Fig. 5h). In the myCAF stage, under the regulation of TFs such as NFYB (Fig. 5g), the levels of the iCAF markers gradually decreased, whereas the chromatin accessibility of genes such as CARMN increased (Fig. 5f, k). The subsequent upregulation of these genes enhanced cell adhesion and cell matrix adhesion (Fig. 5h) and ultimately contributed to tumor development and metastasis through increased TME interactions.

We compared chromatin accessibility and transcriptome profiles between adjacent NFs and CAFs to better understand the process of CAF differentiation. A total of 265 genes were specifically upregulated in tumors because of increased chromatin accessibility (Fig. 5i). Among these genes, 15 genes, including ENO1 and ETV1, are TFs, whereas 21 genes, such as SCARA5, ODF3B and COL1A2, are involved in fibrogenesis evolution and ESCC progression (Fig. 5i).

Chromatin accessibility influences the evolution of tumor-associated macrophages

Tumor-associated macrophages (TAMs) are key players in the TME, and promote tumor growth and metastasis by inducing inflammation, angiogenesis, cancer cell proliferation, and immunosuppression [57, 58]. Based on the marker gene expression in the scRNA-seq data [59] (Supplementary Fig. S3a), we annotated macrophages, and identified 6 macrophage subsets: the Angio-TAM, IFN-TAM, Inflam-TAM, LA-TAM, Prolif-TAM, and RTM-TAM subsets (Fig. 6a). These cell subsets were regulated by different TFs (Supplementary Fig. S3c), and exhibited distinct gene expression (Supplementary Fig. S3f), and functional characteristics (Supplementary Fig. S3b). For example, Angio-TAMs promoted tumor angiogenesis, IFN-TAMs upregulated the expression of IFN regulatory-associated genes, Inflam-TAMs recruited and regulated immune cells during the inflammatory response, LA-TAMs upregulated coagulation pathways, and Prolif-TAMs were closely related to the cell cycle and cell proliferation (Supplementary Fig. S3a, b). RTM-TAMs, characterized by high expression of FOLR2, LYVE1 and HES1, were associated with embryonic precursors and enriched in nearby normal tissues (Supplementary Fig. S3a, d). The annotations of scRNAs and scATAC data were then mapped to further investigate the epigenetic mechanisms of TAMs, and four TAM cell populations were identified in the scATAC data, namely, Angio-TAM, Inflam-TAM, LA-TAM, and RTM-TAM (Fig. 6b, Supplementary Fig. S3e).

Fig. 6

Chromatin accessibility influences the evolution of tumor-associated macrophages. a, b UMAP plot of all macrophages colored by cell subtype according to the scRNA-seq data (a) and scATAC-seq data (b). c Scatter plots showing DEGs between the scATAC and scRNA-seq data of ESCC macrophages in tumor and peritumor tissues. d, e Browser track visualization of the accessibility profile at the SPP1 (d) and PID1 (e) loci across malignant and nonmalignant cells. The putative cancer-specific dREs of SPP1 (d) and PID1 (e) are highlighted in light yellow. f Pseudotime representation of macrophage differentiation (upper) and CXCL3 expression (lower) from RTM-TAMs to inflam-TAMs using scRNA-seq data. The line represents a double-spline fitted trajectory across pseudotime. g Pseudotime representation of macrophage differentiation (upper) and CXCL3 chromatin accessibility (lower) from RTM-TAMs to inflam-TAMs using scATAC-seq data. The line represents a double-spline fitted trajectory across pseudotime. h, i Browser track showing the accessibility profile at the CXCL3 (h) and VEGFA (i) loci across malignant and nonmalignant cells. The putative cancer-specific dREs of CXCL3 (h) and VEGFA (i) are highlighted in light yellow. j Pseudotime analysis of macrophage differentiation (upper) and VEGFA expression (lower) from RTM-TAMs to inflam-TAMs using scRNA-seq. The line represents a double-spline fitted trajectory across pseudotime. k Pseudotime representation of macrophage differentiation (upper) and VEGFA chromatin accessibility (lower) from RTM-TAMs to inflam-TAMs using scATAC-seq data. Pseudotime representation from RTM-TAMs to inflam-TAMs in scATAC data. The line represents a double-spline fitted trajectory across pseudotime

We performed differential peak analysis to clarify the specific cis-regulatory elements of TAMs in ESCC and found that 1236 peaks were upregulated in tumors (Supplementary Fig. S3g). These peaks increased the chromatin accessibility of 1260 genes, of which 138 genes were upregulated in the scRNA-seq data (Fig. 6c). GO analysis revealed the enrichment of these 138 upregulated genes in lysosomal and antigen processing pathways (Supplementary Fig. S3i). Conversely, 29 genes were downregulated and enriched mainly in the adhesion pathway (Fig. 6c; Supplementary Fig. S3i). SPP1 and CCL3L3 exhibited increased chromatin accessibility and transcription in tumor tissues compared with adjacent normal tissues (Fig. 6c, d). However, SPRR3 and PID1 display the opposite pattern (Fig. 6c, e). We also performed a motif prediction of TAMs and mapped the TAM transcription factor regulatory network. We found that the FOS and JUN families were highly enriched in TAMs (Supplementary Fig. S3h).

RTM-TAMs are macrophage subsets resembling normal RTMs [59]. By performing a quasi temporal analysis of tumor-associated macrophages, we identified two differentiation pathways involved in ESCC: RTM-TAMs to Inflam-TAMs and RTM-TAMs to Angio-TAMs. The expression of CXCL3, a classic cell marker of Inflam-TAMs [59], gradually increased during the progression of RTM-TAMs to Inflam-TAMs (Fig. 6f, g). Moreover, 8 cis-regulatory elements within the CXCL3 gene were identified as promoting this cellular transformation (Fig. 6h). In addition, the evolutionary process from RTM-TAMs to Inflam-TAMs can be divided into four distinct stages according to the differences in TFs and cis-regulatory elements. The first stage was characterized by the upregulation of TFs such as EHF and NFYA and the increased chromatin accessibility of genes such as FKBP5. In the second stage, which was regulated by the FOS family and JUN family, the chromatin accessibility of FKBP5 decreased, whereas that of CTNNA3 and NRG3 gradually increased. The third stage involves the upregulation of chromatin accessibility in genes such as DPP6 and EYS under the regulation of TFs such as ZEB1 and ZBTB3. In the fourth stage, under the regulation of the SPI family and the CEBP family, the expression of inflammation-related genes such as CCL3L3 increased, indicating the cellular transition to Inflam-TAMs (Supplementary Fig. S3j, k).

During the differentiation of RMT-TAMs to Angio-TAMs, the chromatin openness score and expression of the VEGFA gene gradually increased in response to cis-regulatory elements (Fig. 6i–k). This process can be divided into three stages. The first stage, the RTM-TAM stage, is regulated primarily by TFs such as SPI, BCL, and ELF, with high expression of PDK4 and DAAM2. The second stage, the transition stage, is regulated by TFs such as TP63, and is characterized by increased chromatin openness of ROBO2 and CTNNA2. In the third stage, under the regulation of TFs such as CEBPA and CEBPD, the expression of RAPGEF1 and HDAC4 is increased, and the cells transform into Angio-TAMs (Supplementary Fig. S3l, m).

Chromatin accessibility landscape of tumor-associated T cells

We annotated T cell subsets based on scRNA-seq marker gene expression to explore the chromatin accessibility landscape of tumor-associated T cells (Supplementary Fig. S4a). Ten T-cell subsets were identified: CD4 + memory T cells, CD4 + naive T cells, CD8 + effector T cells, CD8 + exhausted T cells, CD8 + memory T cells, CD8 + mitotic T cells, CTL, NKT, Tfh, Treg (Fig. 7a). By mapping this annotation to the scATAC-seq data, a total of 9 T-cell subsets were identified (Fig. 7b). Notably, the proportion of CD8 + exhausted T-cell subsets was significantly higher in tumor tissues than in normal tissues according to both the scRNA-seq and scATAC-seq data, which may indicate that CD8 + exhausted T cells play an important role in ESCC (Fig. 7a–c).

Fig. 7

Chromatin accessibility landscape of tumor-associated T cells. a UMAP plot of all T cells colored by subtype based on scRNA-seq data (left). Stacked bar charts showing the contribution of each sample to each T-cell sub type using scRNA-seq data (top right). UMAP plot of all T cells colored by tissue using scRNA-seq data (bottom right). b UMAP plot of all T cells color-coded by subtype using scATAC-seq data (left). Stacked bar charts showing the contribution of each sample to each T-cell subtype based on scATAC-seq data (top right). UMAP plot of all T cells colored according to tissue type based on scATAC-seq data (bottom right). c Scatter plots showing the proportions of CD8 + exhausted T cells in the scRNA-seq (upper) and scATAC-seq (lower) cohorts. d Volcano plots showing differential peaks in CD8 + exhausted T cells between tumor and peritumoral tissues. e Scatter plots showing DEGs between scATAC-seq and scRNA-seq data in tumor and peritumoral tissues among ESCC CD8 + exhausted T cells in ESCC. f KEGG pathway analysis of CD8 + exhausted T cells in the tumor and peritumor scRNA-seq and scATAC-seq cohorts. g Dot plot of the expression levels of marker genes involved in antigen processing and presentation. The dot size corresponds to the percentage of cells expressing the marker gene, and the dot color indicates the average expression level. h Browser track visualization of the accessibility profile at the TNFSF4 locus across tumor and peritumoral cells. The putative cancer-specific dRE of TNFSF4 is highlighted in light yellow. i Heatmap showing the interaction of the OX40 signaling pathway between T-cell subtypes and macrophages in tumor tissues. The bar charts above and on the right represent the sum of the values in each column and row. j Heatmap showing the sender, receiver, mediator and influencer genes in the OX40 signaling pathway. k Violin plots showing the expression levels of OX40 pathway ligands and receptors in different macrophage and T-cell subsets in tumor and peritumoral tissues. l Circle plot showing the interaction between T cells and macrophage subsets of the classic ligand-receptor pair (TNFSF4-TNFRSF4) in the OX40 pathway in tumor tissues

We performed a variance analysis of CD8 + exhausted T cells in adjacent tissues and tumor tissues and detected 834 peaks that were upregulated in CD8 + exhausted T cells in tumors compared with adjacent tissues (Fig. 7d). HLA-DQA1 and TNFSF4 expression increased in tumor-associated CD8 + exhausted T cells, whereas the cis-regulatory elements of IL7R and CXCL8 were inactivated, leading to decreased gene expression (Fig. 7e). We performed a KEGG analysis of genes exhibiting concordant chromatin openness, as shown by scATAC and transcript levels based on scRNA analysis, to further clarify the function of CD8 + exhausted T cells regulated by cis-regulatory elements in ESCC. We found that CD8 + exhausted T cells upregulated the antigen processing and presentation signaling pathway (Fig. 7f) and expressed HLA family antigen-presentation-related cell markers in tumors (Fig. 7g). These findings indicate that while CD8 + exhausted T cells retain antigen presentation capabilities in tumor immunity, their tumor-killing function is compromised.

TNFRSF4, a marker of T-cell activation, is expressed primarily on activated effector T cells and regulatory T cells (Tregs). TNFSF4 is a key regulator of TNFRSF4 expression [60, 61]. In this study, we observed a strong chromatin accessibility signal at the TNFSF4 promoter, accompanied by three cancer-specific regulatory elements. These cis-regulatory elements enhanced TNFSF4 transcriptional activity in CD8 + exhausted T cells (Fig. 7h). We performed cell–cell communication analysis of TAMs and T-cell subsets to further explore the effect of CD8 + exhausted T cells on TNFSF4 upregulation in the TME. The results suggested that tumor cells communicated more specifically through OX40 and other signaling pathways than adjacent tissues did (Supplementary Fig. S4c, d). Notably, OX40 was enriched mainly in CD8 + exhausted T cells and Tregs (Fig. 7i). Under the influence of CD8 + mitotic T cells and CTLs, CD8 + exhausted T cells activate Tregs through TNFSF4-TNFRSF4 ligand-receptor binding, maintaining the tumor in an immunosuppressed state (Fig. 7j–l).

The impact of epigenetic alterations on cell‒cell communication in ESCC

We performed cell communication analysis on the four cell subsets analyzed above, epithelial cells, fibroblasts, macrophages, and T cells, to study their interactions in ESCC. The results revealed that the interactions among these cell types were significantly increased in ESCC tissues compared with normal tissues (Supplementary Fig. S5a). Notably, the enhanced communication in which fibroblasts act as senders communicating with macrophages in tumors was most obvious (Supplementary Fig. S5b).

Further analysis of the cell–cell communication pathways of different cell subsets revealed tumor-specific pathways involving LIGHT, RESISTIN and ANGPT, as well as significant differences in the TWEAK pathway between ESCC and adjacent tissues (Fig. 8a). Previous studies have demonstrated that TWEAK, which is expressed in tumor cells, induces proliferation and survival signaling in tumor cells [62]. In this study, we found that TWEAK was expressed mainly by IFN-TAMs, Angio-TAMs, Prolif-TAMs and inflam-TAMs in peritumoral tissue and acted on fibroblasts such as iCAF3 and myCAF1 cells. However, in tumor tissues, TWEAK was produced by IFN-TAMs and Prolif-TAMs and received by epi1-5 (the subcell type information was published in our previous assay [9]) and myCAFs (myCAFs1 and myCAFs2) (Fig. 8c; Supplementary Fig. S5c). The alteration of this pathway is closely related to the interaction between the TNFSF12 ligand and TNFRSF12A receptors. TNFSF12 is a ligand for both tumor and peritumor TAMs, while TNFRSF12A functions as a receptor that is specifically overexpressed in CAFs and malignant tumor epithelial cells (Fig. 8d; Supplementary Fig. S5d). Interestingly, the chromatin accessibility analysis revealed no significant differences in the degree of chromatin openness of TNFSF12 between tumor and normal tissues. However, compared with normal epithelial cells, TNFRSF12A displayed significantly increased chromatin accessibility in malignant epithelial cells, as well as in fibroblasts (iCAF3 and myCAF1) compared to NFs (Fig. 8d, e). These findings were consistent with the findings of the transcriptome cell communication analysis. These results suggest that altered chromatin accessibility in ESCC causes transcriptional dysregulation of ligands and receptors, ultimately influencing the cell-to-cell communication between ESCC cells and the TME.

Fig. 8

The impact of epigenetic alterations on cell‒cell communication in ESCC. a, b CellChat heatmap showing the initiators (a) and receivers (b) of epithelial, fibroblast, macrophage, and T-cell subsets in tumor (right) and normal (left) tissues, as well as the intensity of action. c Circle plot showing the interactions among the TWEAK signaling pathway in subsets of epithelial cells, fibroblasts, macrophages, and T-cells in peritumoral (upper) and tumor (lower) tissues. d Violin plots showing the expression levels of TWEAK pathway ligands and receptors in different epithelial, fibroblast, macrophage, and T-cell subsets in tumors and peritumor tissues. e UMAP plot colored according to log-normalized gene scores showing the accessibility of cis-regulatory elements linked to TNFRSF12A and TNFSF12 in macrophages (left), epithelial cells (middle) and fibroblasts (right)

Discussion

Previous studies utilizing bulk ATAC-seq and RNA-seq data have demonstrated a correlation between chromatin accessibility and the expression of tumor-specific cell markers in Esophageal Squamous Cell Carcinoma (ESCC). However, these investigations have primarily been conducted at the cellular or tissue level, limiting our understanding of chromatin accessibility's role in single-cell tumor heterogeneity. To address this gap, we performed scATAC-seq on ESCC tissue and integrated the results with scRNA-seq data to explore the impact of chromatin accessibility on tumor heterogeneity at the single-cell level. Our study revealed a strong association between regulatory elements and the transcriptional activation of classic cell markers, emphasizing a pivotal role for chromatin accessibility in defining cellular identity within ESCC. Specifically, a comparison of the cell annotation results of scATAC-seq with those of scRNA-seq revealed that the regulatory element serves as a crucial molecular switch, mediating the expression of defining cell markers. This observation not only deepens our understanding of the epigenetic landscape in ESCC but also points to the intricate interplay between chromatin state and gene expression that underpins cellular differentiation and identity.

In addition, according to the statistical analysis of cell counts in tumor and peritumor tissues, scATAC-seq revealed that some cell populations were absent or altered in population compared to those of scRNA-seq. A similar pattern was observed in a study of gynecologic malignancies [17]. These differences are likely attributable to inherent variations in the technical approaches, emphasizing the need for further optimization of scATAC-seq protocols, including increased sequencing depth and sample sizes, to ensure comprehensive and accurate characterization of the epigenetic landscape. SERPINH1 is overexpressed in a variety of cancers, such as invasive breast carcinoma, cholangiocarcinoma, colon adenocarcinoma, esophageal carcinoma, and glioblastoma [63], and can promote EMT in gastric cancer [64]. SERPINH1 can also promote the metastasis of clear cell renal carcinoma [65] and affect the prognosis of nasopharyngeal carcinoma [66] and adrenocortical carcinoma [63]. However, there have been no previous reports on the mechanism of chromatin regulation of SERPINH1. The openness of the cis regulatory element of SERPINH1 in early ESCC was observed in our study, suggesting a potential role for chromatin accessibility in driving its overexpression. By constructing knockdown and overexpression in cell lines, it was found that SERPINH1 may also participate in tumorigenesis and development through other pathways in addition to EMT. For example, it may affect cell proliferation, apoptosis, autophagy, and other biological processes through p53, regulating the growth and survival of tumor cells. In addition, SERPINH1 may also participate in the immune escape of tumor cells through inflammatory response pathways, reducing the possibility of tumor cells being recognized and cleared by the immune system through interactions with the immune system. The three enhancers of SERPINH1 identified by RNA-seq, ATAC-seq, Chip-seq, CUT-Tag, and Hi-C reveal a new mechanism for SERPINH1 to participate in tumorigenesis and development. Future research will further deepen our understanding of the role of SERPINH1 in tumorigenesis and development, and provide new ideas and methods for tumor treatment.

Differential expression of specific genes at different clinical stages leads to cancer heterogeneity [67]. In this study, there were differences in chromatin accessibility among malignant epithelial cells from different stages of cancer. By analyzing scRNA-seq and scATAC-seq data from clinical samples spanning multiple stages, we gained insights into the dynamic interplay between chromatin openness and gene transcription during tumor development. Our findings revealed that chromatin accessibility is closely related to the level of gene expression, giving rise to distinct chromatin architectures and biomarkers characteristic of each ESCC stage. This underscores the importance of considering chromatin dynamics in the context of tumor progression and highlights the potential for chromatin-based biomarkers in disease staging and therapeutic targeting. Moreover, the identification of stage-specific chromatin markers holds promise for refining the clinical staging of ESCC, enabling more precise patient stratification and personalized treatment strategies.

TME is a complex ecosystem that significantly impacts tumor progression. A notable finding of our investigation is the upregulation of TNFSF4 expression in CD8 + exhausted T cells, a phenomenon closely tied to the activation of regulatory elements. TNFSF4, as a key ligand for TNFRSF4, is well-established for its involvement in T-cell activation [60, 61]. Our data revealed that this upregulation is not merely a byproduct of cellular exhaustion but rather a strategic response within the TME, with implications for the overall immune balance. The overexpression of TNFSF4 in CD8 + exhausted T cells and its subsequent engagement with TNFRSF4 on Tregs suggests a complex cross-talk between these immune subsets. This interaction not only activates Tregs but also highlights the potential for CD8 + exhausted T cells to indirectly contribute to an immunosuppressive microenvironment through their TNFSF4-mediated effects. This finding challenges the traditional view of CD8 + exhausted T cells as purely dysfunctional and emphasizes their nuanced role in modulating the immune response. Importantly, our study highlighted the dual-edged sword nature of CD8 + exhausted T cells in the context of cancer immunology. While they retain some cytotoxic potential, their accumulated TNFSF4 activity seems to skew the immune response towards immunosuppression, potentially facilitating tumor immune escape. Our study provided a comprehensive view of the complex interplay between chromatin accessibility, immune cell function, and tumor progression within the TME, paving new ways for research and therapeutic intervention in the field of cancer immunology.

Additionally, our findings indicated that chromatin accessibility plays a crucial role in mediating CAF-TAM crosstalk, thereby driving ESCC metastasis. Specifically, we proposed that chromatin accessibility promotes RTM-TAM differentiation into Angio-TAMs and Inflam-TAMs, which subsequently secrete TWEAK to activate CAFs [68]. Activated CAFs then release growth factors and cytokines that stimulate distant cancer cell growth and invasion, ultimately enhancing metastasis [69]. This novel insight highlights their active and multifaceted role in promoting tumor progression. By differentiating into Angio-TAMs and Inflam-TAMs, these cells actively contribute to the creation of a permissive microenvironment that supports tumor growth and metastasis. Our findings emphasize the importance of understanding the epigenetic regulation of TAM differentiation and function in the context of cancer immunology.

Conclusions

This study provides novel insights into the role of chromatin accessibility in shaping ESCC tumor heterogeneity, progression, and the TME. By integrating scATAC-seq and scRNA-seq data, we generated a comprehensive resource for understanding the underlying molecular mechanisms driving ESCC pathogenesis. Our findings offer novel insight for developing targeted therapies aimed at modulating chromatin accessibility and disrupting key regulatory pathways.