Gene expression profiles and corresponding clinical features of ESCC patients were obtained from the GSE53625 (179 tumor tissues and 179 normal tissues) and GSE121931 (125 tumor tissues) cohorts from the Gene Expression Omnibus (GEO) database. Raw proteomics data and corresponding clinical features from PXD021701 (124 ESCC tissues and 124 normal tissues) were obtained from the PRoteomics ID Entifications Database (PRIDE Archive). ESCC scRNA-seq data were acquired from GSE197677 (10 solid tumor tissues and four normal tissues). Microbial abundance profiles at different taxonomic levels were obtained from the Cancer Microbiome Atlas (TCMA, https://tcma.pratt.duke.edu/) database for 40 ESCC tumor samples. A total of 66 tryptophan metabolism-related genes (TrpGs) were obtained from the Human Gene Set “KEGG_TRYPTOPHAN_METABOLISM” in the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb), and from the “superpathway of tryptophan utilization” in PathCards (https://pathcards.genecards.org/Card/superpathway_of_tryptophan_utilization). The TrpGs are listed in Supplementary Table S1.

Differentially expressed genes (DEGs) between ESCC tumor tissues and normal tissues (GSE53625 cohort) were identified using the “Limma” R package with the criteria of |log(Fold change, FC)| > 0.5 and P-value < 0.01. Then, the differentially expressed TrpGs (DE-TrpGs) were selected by overlapping the DEGs in ESCC and TrpGs from the public databases.

The scRNA-seq data were processed and assessed using the “Seurat” R package. Quality control, normalization, and clustering were finished before data acquisition (Okuda et al., 2023). Therefore, a total of 29,473 genes and 67,626 cells, and an expression matrix were included for subsequent analyses in this study. Cells were annotated using known cell marker genes. The “COSG” R package was used to implement Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for each cell type, based on the top 2,000 marker genes in each cell type. We also detected the expression of DE-TrpGs in each cell population. The TrpG score for each cell population was calculated using the “UCell” R package based on 66 TrpGs. The cell population with the highest TrpG score was selected for further analysis.

The “CopyCAT” R package was used to analyze the copy number variant (CNV) characteristics of the epithelial cell populations and distinguish tumor cells from normal ones. The “COSG” R package was used for GO enrichment analysis based on the marker genes. Then, the “CellChat” R package was used to investigate the potential interactions between tumor cells and high-TrpG score cells. The CellChat analysis was complementarily used for cell communication analysis based on the CellChat DB, which contains 2,021 validated interactions, including 60% of paracrine/autocrine signaling interactions, 21% of extracellular matrix (ECM)–receptor interactions, and 19% of cell–cell contact interactions (Jin et al., 2021).

The least absolute shrinkage and selection operator machine learning algorithm (LASSO) regression was used to select the prognostic TrpGs in the GSE121931 cohort using the “glmnet” R package. Univariate Cox regression analysis was used to screen protective factors for ESCC among the above LASSO-generated risk factors in the GSE121931 and PXD021701 cohorts. The hazard ratio (HR) and 95% confidence intervals (CI) were calculated. TrpGs with HR > 1 were identified as risk factors, and those with HR < 1 were identified as protective factors. A bilateral p < 0.05 was considered statistically significant. Based on these protective factors, the “ggrisk” R package was used to calculate a risk score and construct a TrpG risk model in the GSE121931 cohort. The ESCC patients were distributed into high-risk score and low-risk score groups based on the median value of the risk score. Kaplan–Meier overall survival (OS) curves were drawn using the “survival” R package for the training (GSE121931) and validation (PXD021701) cohorts. Receiver operating characteristic (ROC) curves for 1-, 2-, and 3-year survival probabilities for ESCC patients were drawn using the “timeROC” R package. Then, a nomogram was established using the “rms” R package according to independent prognostic factors. The calibration curves and decision curves were drawn to verify the accuracy of the nomogram.

The marker gene set of 28 immune cells was obtained and referenced from a previous article (Charoentong et al., 2017). The “IOBR” R package was used to estimate the proportion of immune cells in the high- and low-risk groups. The correlation between infiltrated cells and risk scores was also assessed by the “IOBR” R package. The immune score, stromal score, and ESTIMATE score were calculated using the ESTIMATE algorithm in the “IOBR” R package. The differences in each score between the high- and low-risk groups were determined by Student’s t-test, ***p < 0.001, **p < 0.01, and *p < 0.05. The differences in the expression of immune checkpoint genes between the high- and low-risk groups were determined by Student’s t-test (***p < 0.001, **p < 0.01, and *p < 0.05), and the correlation between the expression of immune checkpoint genes and risk scores was detected using a Spearman’s correlation analysis.

The NCI-60 compound activity data and RNA-seq expression profiles from the CallMiner™ (https://discover.nci.nih.gov/cellminer/home.do) were obtained to investigate the drug sensitivity. Drugs approved by the Food and Drug Administration (FDA) or in clinical trials were selected for analysis. The “pRRophetic” R package and ridge regression analysis were used to predict and evaluate the half maximal inhibitory concentration (IC50) value for each sample. The differences in IC50 value between the high- and low-risk groups were detected by the Student’s t-test (***p < 0.001, **p < 0.01, and *p < 0.05).

The microbial abundance profiles at different taxonomic levels for 40 ESCC tumor samples were obtained from the TCMA database. The differences in microbiota were assessed using a Student’s t-test (***p < 0.001, **p < 0.01, and *p < 0.05).

A total of 10 ESCC tumor samples and 10 normal samples were obtained from patients who had been diagnosed with ESCC and had not received any preoperative treatment at the Shandong Public Health Clinical Center. The protocol for collecting clinical samples was approved by the Ethics Committee of the Shandong Public Health Clinical Center (approval number: 2021 XKYYEC-37), and the patients provided informed consent before samples were collected. The patients underwent surgical resection, and the samples were frozen in liquid nitrogen immediately after collection and stored at −80°C for further RNA extraction.

TRIzol reagent (Takara, Dalian, China) was used to isolate the total RNA from tissues, following the standard protocol. RNA quantification was performed with a Nanodrop (Thermo Fisher Scientific, MA, United States). To compare expression levels between two groups, quantitative real-time PCR (qRT-PCR) was performed using the RNA-direct SYBR Realtime PCR Master Mix (Takara) and a Stratagene Mx3000P Real-Time PCR System (Agilent Technologies, CA, United States). The expression of each gene was determined by using GAPDH as the reference gene. Each experiment was performed in triplicate to obtain the average value. The qRT-PCR results were analyzed with the 2^−ΔΔCt method. The primers used are: GAPDH: CCTCAACTACATGGCTGAGAAC and CAAGGGGTCTACATGGCAACT, MAOA: CCAGCGGTAGAAATCACCCA and TCTGATGAGCACATACACGTTACTT, AKR1A1: GGGTACCTGGAAGAGTGAGC and GATCTGAGCTGGAGATCGGC, ALDH9A1: CAACCGGCCGAGTGATAGC and TGTGGTCGGTTGATGAGTGG, HAAO: AACAAGCTCATGCACCAGGA and CATGGTGTCGCCCACATAGT, ALDH2: TAATCCAGGTTGCTGCTGGG and AGTGCTCACCTCCTCCTTG.

All experimental results are expressed as the mean ± standard deviation (SD). Statistical analyses were performed using GraphPad Prism software, version 8.3.0 (GraphPad Software, CA, United States). The significant difference was determined using a two-tailed Student’s t-test (two-group comparisons). P < 0.05 indicated statistical significance.

The workflow of this study is shown in Figure 1. A total of 6,029 DEGs (3,074 upregulated and 2,955 downregulated) between ESCC samples and normal samples in the GSE53625 cohort were identified with the criterion of |log(FC)| > 2 and p < 0.01 (Figure 2A; Supplementary Table S2). Then, 28 DE-TrpGs were obtained by overlapping the 6,029 DEGs and 66 TrpGs (Figures 2B, C), which contained eight upregulated (STAT1, IL4I1, TDO2, NMNAT3, KMO, IDO1, OGDHL, and ACAT2) and 20 downregulated DE-TrpGs (ALDH9A1, ALDH3A2, ADH1B, AOX1, MAOA, ECHS1, ALDH2, ACAT1, HAAO, ALDH7A1, HADHA, MAOB, AKR1A1, CYP4X1, INMT, CYP2U1, CAT, HADH, NMNAT1, and UGT2B11) between ESCC and normal samples (Figure 2B). These results indicate that tryptophan metabolism is potentially involved in the ESCC tumorigenesis and progression.

We further investigated the influence of tryptophan metabolism in the tumor microenvironment (TME) of ESCC based on single-cell gene expression profiles. Eight primary cell populations were identified with marker genes in ESCC (Figures 3A, B), including 25,627 T cells, 1,477 mast cells (Mast), 3,324 natural killer (NK) cells, 7,998 myeloid cells (Myeloid), 14,425 fibroblasts (Fib), 1,574 endothelial cells (Endo), 8,344 B cells, and 4,857 epithelial cells (Epi). We also found an increased abundance of immune cells (T cells and B cells) and epithelial cells, and a decreased abundance of fibroblasts in ESCC tumor samples compared with normal samples (Figure 3C). GO and KEGG enrichment analysis illustrated the different biological functions in the different cell populations (Figures 3D, E). Especially, extracellular matrix structural constituents, collagen-containing extracellular matrix, focal adhesion, and ECM-receptor interaction were significantly enriched in fibroblasts. The T cell receptor signaling pathway was significantly enriched in the T cell population. We also found the fibroblast population with the highest TrpG score in the tumor microenvironment (TME) of ESCC (Figure 3F), and the TrpG score in ESCC tumor samples was lower than in normal samples (Figure 3G). We next explored the expression of 28 DE-TrpGs, from which bulk RNA-seq data analysis in different cell types showed that 18 DE-TrpGs (ALDH9A1, ALDH3A2, ADH1B, AOX1, MAOA, ECHS1, ALDH2, ACAT1, HAAO, ALDH7A1, HADHA, AKR1A1, CYP4X1, INMT, CYP2U1, CAT, HADH, and NMNAT1) were significantly expressed at the single-cell level, especially the downregulated expression of fibroblasts in ESCC tumor samples compared with normal samples (Figure 3H). The single-cell RNA-seq results were consistent with the bulk RNA-seq results: 18 DE-TrpGs were downregulated in ESCC tumor samples compared with normal samples and were clearly observed in fibroblasts.

The TrpG scores of fibroblasts were found to be significantly different between the normal groups and the ESCC group. Therefore, we focused on the fibroblast subtypes in ESCC. As shown in Figure 4A, four subtypes of fibroblasts were identified, including normal fibroblasts (normal_Fib), myofibroblasts (myCAF), antigen_CAFs, and perivascular‐like fibroblasts (PVL). The subtypes of normal_Fib were mainly enriched in normal samples, and the others were significantly enriched in ESCC samples (Figure 4B). Furthermore, 564 malignant epithelial cells (aneuploid) and 83 benign epithelial cells (diploid) were identified by calculating the CNV using the “copyCAT” R package (Figures 4C, D). Biological functional enrichment analysis indicated that malignant epithelial cells were significantly enriched in MYC targets, MTORC1 signaling, interferon alpha/gamma response, and PI3K/AKT/MTOR signaling pathways (Figure 4E).

Then, we created a communication network to illustrate the communication between fibroblasts and epithelial cells. ECM–receptor analysis results indicated that malignant epithelial cells mainly interacted with normal fibroblasts and myofibroblasts, but a weaker association was found between benign epithelial cells and myofibroblasts (Figure 4F). Heatmaps displayed the normal fibroblasts and myCAF, and malignant epithelial cells represented the main signaling providers and the receptors, respectively (Figure 4G). The intercellular interactions were mainly active in signaling pathways, including pathways involving LAMNIN, HSPG, and ARGN (Figure 4G). Furthermore, secreted signaling analysis results indicated a strong interaction between malignant epithelial cells and normal fibroblasts and myofibroblasts and a weaker interaction between benign epithelial cells and myofibroblasts (Figure 4H). Heatmaps also represented normal fibroblasts or myofibroblasts as the main signaling providers, and the malignant epithelial cells were represented as the main receptors (Figure 4I). The intercellular interactions were mainly involved in activation of several signaling pathways, including MK, PTN, EGF, and IGF (Figure 4I). Taken together, a strong interaction between malignant epithelial cells and normal fibroblasts or myofibroblasts and a weak interaction between benign epithelial cells and fibroblasts or myofibroblasts indicate that communication between malignant epithelial cells and normal fibroblasts or myofibroblasts is transformed in the TEM of ESCC.

Based on bulk RNA-seq data and scRNA-seq data analyses, 18 DE-TrpGs (ALDH9A1, ALDH3A2, ADH1B, AOX1, MAOA, ECHS1, ALDH2, ACAT1, HAAO, ALDH7A1, HADHA, AKR1A1, CYP4X1, INMT, CYP2U1, CAT, HADH, and NMNAT1) were screened in ESCC, then incorporated into a LASSO regression model to select the diagnostic TrpGs in ESCC. As a result, 13 TrpGs (ALDH9A1, ADH1B, AOX1, MAOA, ECHS1, ALDH2, ACAT1, HAAO, ALDH7A1, HADHA, AKR1A1, CAT, and HADH) were identified in the GSE121931 cohort (Figures 5A, B). In addition, univariate Cox analysis was performed to identify the protective factors for ESCC: 10 protective factors (MAOA, INMT, ALDH7A1, AOX1, AKR1A1, ALDH9A1, ACAT1, HAAO, ALDH2, and HADH) were identified in the GSE121931 cohort (Figure 5C). Ten protective factors (HAAO, KYAT3, ADH1B, ALDH2, ALDH9A1, ALDH9A2, MAOA, CAT AKR1A1, ECHS1) were identified in the PXD021701 cohort (Figure 5C). Finally, five protective factors (MAOA, AKR1A1, ALDH9A1, HAAO, and ALDH2) were obtained by overlapping the above results (Figure 5D).

A TrpG risk score was calculated, and the patients were distributed into high-risk and low-risk groups based on the median risk score in both the GSE121931 and PXD021701 cohorts (Figures 5E, H). Kaplan–Meier OS curves indicated that patients with high-risk scores have a poor survival time in both the GSE121931 and PXD021701 cohorts (Figures 5F, I). AUC values for 1-, 2-, and 3-years were 0.630, 0.651, and 0.633 in the GSE121931 cohort, respectively (Figure 5G). AUC values for 1-, 2-, and 3-years were 0.629, 0.620, and 0.600 in the PXD021701 cohort, respectively (Figure 5J). These findings suggested that the five-TrpG signature had good prognostic value in both cohorts. Multivariate Cox regression analysis was used to identify the independent prognostic factor for ESCC patients by incorporating the clinical clinicopathological characteristics (tumor grade and lymph node stages), and the five-TrpG signature, HAAO, ALDH2, and lymph node stages were identified as the independent factors for ESCC to generate a nomogram for 1-, 2-, 3-year survival prediction (Figure 5K). The calibration curve indicated the good performance of the nomogram for 2-year survival prediction (Figure 5L). The decision curves indicated the sensitivity of the factors (Figure 5M). Additionally, we collected ESCC samples to examine the expression levels of five protective factors (MAOA, AKR1A1, ALDH9A1, HAAO, and ALDH2). The results revealed that the expression of MAOA, AKR1A1, ALDH9A1, HAAO, and ALDH2 was significantly downregulated in ESCC tumor tissues compared to non-tumor tissues (Figure 6). These findings suggest that these five genes may serve as potential protective factors and hold promise as diagnostic biomarkers for ESCC.

We next evaluated the correlation between the TrpG risk score, immune characteristics, and drug sensitivity. The results indicated that the high enrichment scores of regulatory T cells (Treg), CD56 dim natural killer (NK) cells, eosinophils, type 17 T helper (Th17) cells, activated CD8 T cells, neutrophils, and activated B cells were significantly enriched in the high-TrpG risk score group compared with the low-TrpG risk score group. However, CD56 bright (NK) cells, activated dendritic cells (DCs), effector memory CD4 T cells, gamma delta (γδ) T cells, central memory CD4 T cells, and immature DCs were significantly decreased in the high-TrpG risk score group compared with the low-TrpG risk score group (Figures 7A, B). We also found a high immune score in the high-TrpG risk score group compared with the low-TrpG risk score group (Figure 7C). Moreover, we found high expression of immune checkpoint genes (KIR2DL1, LGALS9, TNFRSF18, and TNFRSF4) in the high-TrpG risk score group compared with the low-TrpG risk score group (Figure 7D). The above findings suggest that the TrpG risk score is positively associated with immune cell infiltration; a high TrpG risk score possibly indicates more immune activation. We also estimated the sensitivity to drug therapy, resulting in ESCC patients with high-TrpG risk scores being less sensitive to drugs, such as lapatinib, gemcitabine, etoposide, erlotinib, elesclomol, doxorubicin, docetaxel, dasatinib, cisplatin, and bleomycin; however, they were more sensitive to imatinib and axitinib (Figure 7E).

Furthermore, to explore whether there were specific candidate microbial abundances associated with the TrpG risk score. As shown in Figure 8, we discovered that the relative abundance of Fusobacteria and Tenericutes was decreased in the high-TrpG risk score group compared to the low-TrpG risk score group. These findings indicate that microbes may be involved in the progression of ESCC.