A total of 23 patients with a definite diagnosis of ESCC were enrolled from March 2021 to December 2021 in this study. Clinicopathological information, including demographics, pathologic diagnoses, imaging examinations, and treatment history were collected from each patient. A total of 8 patients received radiotherapy alone and 15 patients received a combination of radio-chemotherapy, as detailed in Table 1 , the chemotherapy regimen was based on clinical guidelines, mainly with nedaplatin and Tegafur Gimer. Tumor tissue samples were collected from all participants who received radiotherapy to perform the whole exome (WES) and transcriptome (WTS) sequencing and matched peripheral blood samples which included three-time points, pre-treatment, on-treatment, and post-treatment were collected to perform the TCR-βsequencing. All procedures were conducted by the Declaration of Helsinki. This study was approved by the Ethics Committee of Shanxi Provincial Cancer Hospital (Taiyuan, China) (Approval No: KY2022011) and written informed consent was obtained from all participants.

DNA- and RNA-based NGS were performed in Geneplus-Beijing (Beijing, China) using WES and WTS as published in our previous study (17, 18). Genomic DNA (gDNA) from WBCs and tumor tissues were processed into indexed libraries using a DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). RNA was extracted from the tumor samples using an RNeasy FFPE Kit (Qiagen, Hilden, Germany). Sequencing libraries of genomic DNA and mRNA were prepared using the KAPA DNA Library Preparation Kit (Kapa Biosystems, MA, USA) and NEB Next Ultra™ RNA Library Prep Kit (Illumina, Inc., CA, USA), respectively. The DNA and RNA sequencing were performed and all experimental procedures followed the manufacturer’s instructions. The DNA and RNA indexed libraries were sequenced using a 100-bp paired-end configuration on a DNBSEQ-T7RS sequencer (MGI Tech, Shenzhen, China) or Gene+Seq-2000 sequencing system (GenePlus, Suzhou, China), producing 30G(150X), 12G (150X), 12G sequencing data for tumor tissues (WES), WBLs (WES), and tumor tissues (RNA), respectively. TCR sequencing and data analysis were performed in Geneplus-Beijing (Beijing, China) as previously described (14, 16). Multiplex PCR amplification of CDR3 of the TCR-β chain (TRB) was conducted including PCR1 and PCR2, inclusively and semi-quantitatively. Libraries were sequenced using the PE150 strategy on the Gene+Seq-2000 sequencing system (GenePlus, Suzhou, China), producing 2G/sample. Based on the ImMunoGeneTics (IMGT) V, D, and J gene references, the CDR3 sequence is characterized as the amino acids situated between the second cysteine in the V region and the conserved phenylalanine in the J region. The MiXCR software package is employed to identify and allocate CDR3 sequences (19).

After the removal of terminal adaptor sequences and low-quality reads (>50% N rate, >50% bases with Q<5) by FASTP (v0.12.6) (20), the remaining reads were aligned to the reference human genome (hg19) and aligned using BWA (version 0.7.10) (21) and HISAT (22) for DNA and RNA sequencing, respectively. Duplicated reads were removed using the MarkDuplicates tool in Picard (version 4.0; Broad Institute). Genomic single nucleotide variants (SNVs), small insertions and deletions (InDels), copy number variants (CNVs), and structural variants (SVs) were called with default parameters by MuTect (version 1.1.4) (23)/NChot, GATK (v3.6-0-g89b7209; Broad Institute) and CONTRA (version 2.0.8) (24), respectively. Transcript assembly was performed using StringTie (version 1.2.3) (25). Chromosomal CNV burden represented the total level of amplifications or deletions at the chromosome level. Significantly recurrent regions with amplification or deletion were detected using Gistic 2.0 with a noise threshold of 0.3, a broad length cutoff of 0.5 chromosome arms, a confidence level of 95%, and a copy-ratio cap of 1.5 (26). The mutational landscape was portrayed using the R package ‘maftools’ (version 2.14.0). The CDR3 sequences were identified and assigned using the MiXCR software package (version 3.0.3) (19). The relative abundance or distribution of each clonotype. Shannon’s entropy was calculated on the clonal abundance of all productive TCR sequences. Clonality score is defined as 1- (Shannon index)/ln(# of productive unique sequences) (16).

Non-synonymous SNVs and Indels with a mutant allele frequency greater than 5% per megabase in the coding region were included in the calculation of Tumor mutation burden (TMB). R package ‘yapsa’ (version 1.24.0) was performed to infer the composition of known Catalogue of Somatic Mutations in Cancer (COSMIC) mutational signatures in EC using the COSMIC mutational signatures version 2 (27). Normalization, estimation of dispersion, and statistical testing of differential expression were performed using the DESeq function in DESeq2 (version 1.38.3) with default parameters. Genes with an adjusted p-value (q-value) < 0.05 and an absolute log2 fold-change (log2FC) > 1 were considered as significantly differentially expressed. R package ‘clusterProfiler’ (version 4.7.1.3) was used to perform gene set enrichment analysis (GSEA) enrichment analysis (28) using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (29) and the Gene Ontology Biological Process (GO BP) database (30). Cold and hot tumors were distinguished based on the immune infiltration scores of samples. The approach involved evaluating immune cell scores using the Gaussian calculation of the ssGSEA (28) method from the standardized TPM matrix using GSVA (31). Immune infiltrating cells were activated, suppressed, and classified into other cell types based on their cellular characteristics. The TME deconvolution method from the R package “IOBR” was utilized to assess the immune microenvironment, and the function “iobr_cor_plot” was employed to compute the expression of features or genes and generate plots. The “iobr_cor_plot” function dynamically generates statistical results by processing the calculated score values through the scale function (32).

Statistical analysis and visualization were performed using the software R 4.3.2 (R Foundation for Statistical Computing, Vienna, Austria). Sample clustering to distinguish immune-cold and immune-hot tumors was achieved through the Euclidean method within the ConsensusClusterPlus function in software R. For the comparison of TCR dynamics, immune microenvironment, genomic characteristics, and PFS (progression-free survival), as well as the NLR (neutrophil-to-lymphocyte ratio) between groups, a t-test will be performed for variables following a normal distribution, while the Wilcoxon rank-sum test will be applied for non-normally distributed variables.Correlation analysis of different indicators was performed using the Pearson correlation method. A two-tailed P < 0.05 was considered statistically significant.

In this study, a total of 23 patients with pathologic diagnoses of locally advanced ESCC were enrolled, and their clinicopathologic features and demographics were summarized in Table 1 . The median age at diagnosis was 69 years, ranging from 53 to 81 years, with 47.8% (11/23) males. Among them, 43.5% (10/23) and 26.1% (6/23) patients had a smoking and drinking history respectively, and none of the patients had a family history. Of all, 30.4% (7/23) patients had stage II disease, 39.2% (9/23) had stage III and 30.4% (7/23) with stage IV. And 82.6% (19/23) of patients were identified as positive for lymph node metastasis. The lesions were located in the cervical segment in 3 cases, the thoracic segment in 19 cases, and both the stiff segment and thoracic segment in 1 case. The median neutrophil-to-lymphocyte ratio (NLR) was 1.79, ranging from 0.93 to 8.47. Radiotherapy and radio-chemotherapy were conducted in 34.8% (8/23) and 65.2% (15/23) of patients, respectively. Based on a rigorous evaluation of radiological evidence (MRI and radiography) by two independent radiologists, 19 patients achieved complete response (CR) or partial response (PR), and 4 patients were diagnosed with stable disease (SD). The median progression-free survival (PFS) was 14 months, ranging from 1 to 27 months.

Using WES data, a total of 3517 somatic mutations were detected in 23 patients with ESCC, including 3520 SNVs and 197 Indels. The median number of mutations was 158 (range 22 - 438). TP53, TTN, PCLO, FAT1, MUC16 and SYNE1 were the top 6 commonly SNVs mutated genes, and mutated in 87%, 57%, 30%, 26%, 26%, and 26% of ESCC, respectively ( Figure 1A ). We detected 143 gene amplifications in these samples, including 89 gene gain alterations and 54 gene loss alterations. The most frequent gene amplifications were of PLD1, TMEM212, FNDC3B, GHSR, TNFSF10, NCEH1, ECT2, SPATA16 and NLGN1, and the amplification of these genes were all located on chromosome 3q26.31 (Data not shown).

A total of 11 significantly amplified and 6 deleted regions were identified by Gistic 2.0 (26) with the q value < 0.1 ( Figure 1B ). We also examined the composition of six possible base-pair substitutions by the Catalogue of Somatic Mutations in Cancer (COSMIC) mutational signature analysis, we found that nearly 50% of mutations are C>T ( Supplementary Figure S1A ), and 75% (41% + 34%) of mutational signatures are attributed to either signature 1 (spontaneous deamination of 5−methylcytosine) or signature 2 (APOBEC Cytidine Deaminase), both of which are associated with C>T mutations ( Supplementary Figure S1B ).

The median TMB of the 23 ESCC patients was 3.47 mutations/Mb (range 1-52.86 mutations/Mb). We compared the TMB levels of the 23 ESCC patients in this study with TMB from other cancers derived from TCGA and found that ESCC had a relatively high level of TMB ( Supplementary Figure S2 ). Subsequently, we allocated mutant genes to 10 hallmark oncogenic pathways. The TP53, RTK/RAS and NOTCH pathways were consistently prevalent in ESCC (33, 34). Besides, some less prevalent pathways, including WNT and HIPPO pathways also exhibited superior frequencies in ESCC (35). Among distinct pathways, multiple RTK/RAS alterations and WNT alterations tended to be concurrent in one patient ( Figure 1C ).

All ESCC biopsies were subjected to RNA sequencing (RNA-Seq). Unsupervised clustering by the euclidean method within the ConsensusClusterPlus function based on immune infiltration was used to classify the 23 ESCC patients into two clusters, the immune-hot group and the immune-cold group ( Figure 2A ). By comparing the differences in the immune microenvironment between immune-hot and immune-cold tumors, we found that tumors from the immune-hot group were highly infiltrated with B-cells, macrophages, CD45, CD8+ T-cells, cytotoxic cells, neutrophils, NK cells, T-cells, and Th1 cells ( Figure 2B , Supplementary Figure S3 ). To validate the stability of the results, we extracted RNA-seq data from 96 esophagus cancer (ESCA) samples in the TCGA database and performed a comparison of the immune microenvironment between immune-hot and immune-cold tumors, which yielded similar results. In the ESCA analysis, hot tumors also exhibited high infiltration of CD8+ T-cells, macrophages, and NK cells ( Supplementary Figure S4 ).

In our study, Differential gene expression analysis was performed on the raw expression counts of all genes between immune-cold group and immune-hot group. Thresholds for the adjust p-value (padj) and |log2FC| were set at <0.05 and ≥1, respectively. As a result, a volcano plot revealed 382 genes significantly upregulated in immune-hot group tumors, and 446 genes significantly upregulated in immune-cold tumors, with adjusted p-values less than 0.05 ( Figure 2C ). We next compared the mRNA expression profiles and signal pathway enrichment between the two groups. Olfactory transduction, Cytokine−cytokine receptor interaction, Neuroactive ligand−receptor interaction and JAK-STAT signals were significantly enriched in the immune cold group ( Figure 2D ). GSEA enrichment analysis was conducted to identify pathways upregulated in cold and hot tumors. Pathways related to the detection of stimulus (NES = 1.48, padj = 2.2 x 10-2) were upregulated in the immune-hot tumors ( Supplementary Table S1 ).

To explore the relationship between genomic and transcription characteristics and prognosis, we divided patients into two groups: PFS-low group (PFS ≤ 14m, n=12) and PFS-high group (PFS > 14m, n=11). Fisher’s exact test was used to find the differentially mutated genes between the PFS-low group and the PFS-high group. The result was not significant. We found the median TMB in the PFS-low group was significantly higher (p=0.049). However, no significant difference in TNB (p=0.15), chromosomal CNV burden (p=0.065), and MSI score (p=0.16) between two groups ( Figure 3A ). The composition of COSMIC mutational signatures was also compared between the PFS-high group and the PFS-low group ( Figure 3B ). Mutational signature 1 was present in both the PFS-high group and the PFS-low group and it was the highest in both groups. The PFS-high group was mainly associated with signature 1 (40%) and signature 4 [36%, exposure to tobacco (smoking) mutagens], and the PFS-low group were a signature 1 (43%) and signature 13 [34%, APOBEC Cytidine Deaminase (C>G)] ( Figure 3B ). The analysis suggests that patients with a better prognosis tend to have a lower TMB.

Based on RNA-seq transcriptome data, we next compared the pattern of gene expression between PFS-high group and PFS-low group using the Expression data (ESTIMATE) algorithm to calculate the ImmuneScore, StromalScore and ESTIMATEScore. We found no differences in these areas between the two groups ( Supplementary Figure S5A ). We then assessed differences in immune cell composition between two groups and found no difference in the ratio of 22 immune cells ( Supplementary Figure S5B ). This suggests that the immune composition may not differ significantly between patients with different prognosis. Different immune microenvironment may have limited influence on prognosis of ESCC.

We next used the systemic inflammation biomarker, NLR, to stratify patients and further analyze differences in genomic and transcription characteristics between different NLRs. We divided patients into two groups: the NLR-low group (0.93≤NLR ≤ 1.79, n=12) and the NLR-high group (1.93≤NLR ≤ 8.47, n=11). KMT2D was mutated in 5 samples (5/12) of the NLR-low group and no mutation was detected in the NLR-high Group ( Supplementary Figure S6A ). Meantimes, we found KMT2D mutant samples tended to have higher TMB but no significant difference ( Supplementary Figure S6B ). No significant difference in TMB, TNB, and chromosomal CNV burden between NLR-low group and NLR-high group (p = 0.93, p = 0.15 and p = 0.065, respectively) ( Supplementary Figure S6C ). The composition of COSMIC mutational signatures was also compared between the NLR-high group and the NLR-low group. The NLR-high group was mainly associated with signature 1 (62%), and NLR- low group was signature 1 (43%) and 2 (30%) ( Supplementary Figure S6D ).

In addition, we calculated the ImmuneScore, StromalScore, and ESTIMATEScore comparing the NLR-high group and the NLR-low group. StromalScore in the NLR-low group was significantly higher than that in the NLR-high group (p = 0.032) ( Figure 3C ). There was no significant difference in the ESTIMATEScore and ImmuneScore between two groups (p = 0.079, p = 0.13) ( Figure 3C ). The proportion of 22 kinds of immune cells in the two groups was analyzed. It was found that the proportion of B_cells_naive and T_cells_CD4_memory_activated in the NLR-low group was significantly higher than that in the NLR-high group, and there was no significant difference in other immune cells ( Supplementary Figure S7 ). This suggests that naïve B cells and activated CD4+ memory T cells and NLR may have some correlation with NLR levels ( Supplementary Table S2 ).

We compared TCR clonality, Simpson index (a type of diversity index, the probability that two clones randomly sampled belong to the same population, and the lower the Simpson index value, the higher the diversity), Shannon index (a type of diversity index, and the lower the Shannon index value, the higher the clonal diversity), and richness among the pre-treatment (pre-treat), on-treatment (on-treat) and post-treatment (post-treat) of ESCC patients by t test using matched samples. Clonality was significantly increased from pre-treat to post-treat (p = 0.016) and from on-treat to post-treat (p = 0.03) ( Figure 4A ). Shannon index was significantly decreased from pre-treat to post-treat (p = 0.0009) and from on-treat to post-treat (p = 0.029) ( Figure 4B ). Richness was significantly decreased from pre-treat to post-treat (p = 0.019) ( Figure 4C ). There was no significant difference in the Simpson index among the pre-treat, on-treat, and post-treat ( Figure 4D ). Two patients exhibited a continuous increase of TCR CDR3 (CHCLPAED_AGGGELFF) frequency in post-treatment samples compared to on-treatment samples. Three patients had elevated levels of CDR3 amino acid sequence (CASSLDSNQPQHF) after treatment. The lower the Shannon Index, the higher the diversity of TCR, and the diversity of TCR decreased after treatment. This suggests that radiotherapy in ESCC patients may result in reduced TCR diversity. To further validate our findings, we performed a similar analysis in the GSE120101 dataset, which showed comparable results. Specifically, Clonality and Simpson index exhibited an upward trend post-treatment, while Shannon index and Richness both decreased from pre-treatment to post-treatment, consistent with our study’s findings ( Supplementary Figure S8 ).

We assessed the relationship between CDR3 diversity and clinical and molecular characteristics in ESCC patients. As immune status can be reflected by TCR diversity, we first evaluated TCR Clonality, Simpson index, Shannon index, and Richness differences between patients with different NLRs during the course of treatment to further investigate how TCR diversity in peripheral blood reflects immune status in ESCC. Unfortunately, there was no difference between the NLR-high group and the NLR-low group in the pre-treat, on-treat, or post-treat samples ( Supplementary Figure S9 ). We continued to explore the prognostic value of the TCR repertoire for ESCC patient outcomes post-radiotherapy. To further demonstrate the similarity between the PFS-high and the PFS-low group, we focused on diversity in the Clonality, Simpson Index, Shannon index, and richness analysis, for each patient were analyzed during treatment. We found that the richness of PFS-low groups was significantly higher than that of PFS-high groups in the pre-treat sample (p = 0.018) ( Figure 4E ), while there was no significant difference in the on-treat and post-treat samples ( Supplementary Figure S10 ). There was no significant difference in the Clonality, Simpson index, and Shannon index between the PFS-high and the PFS-low groups ( Supplementary Figure S10 ). Subsequently, the changes in the TCR Diversity of 23 ESCC patients were analyzed. The results showed that patients with a significant increase in TCR-Diversity had better PFS and sustained clinical benefit in the PFS-high group (p = 0.024) ( Figure 4F ). In the correlation analysis, we found that clone type was significantly negatively correlated with NLR (p < 0.01), and the Shannon index was negatively correlated with NLR ( Figure 4G ).

To examine the effect of gene mutations on the TCR clonal pattern, we assessed the clonal differences between the wild-type and mutant variants, employing an unpaired t-test for statistical comparison. We found that the baseline clonality in mutant FBXW7 (MT) was significantly lower than that in wild-type (WT) patients (p = 0.0076), while they were not significant after treatment between FBXW7-MT and WT patients (p = 0.059 in on-treatment samples; p = 0.18 in post-treat samples) ( Supplementary Figure S11 ). Intriguingly, in FBXW7-MT patients, the Clonality that underwent therapy was higher than that at baseline. During radio-chemotherapy, the clonality in RYR1-MT and UNC79-MT patients had a similar phenomenon to the FBXW7-MT patients. The baseline clonality in RYR1-MT and UNC79-MT were significantly lower than that in WT patients (p = 0.0073 for RYR1 and p = 0.033 for UNC79), while they were not significant after treatment between MT and WT patients ( Supplementary Figure S11 ). The results suggested that FBXW7, RYR1, and UNC79 mutant patients might have better treatment outcomes upon radio-chemotherapy.