A prognostic model based on 10 gene signatures associated with intestinal microbiota predicts survival prognosis of esophageal squamous cell adenocarcinoma

Abstract

This study aimed to construct a potential prognostic model for the clinical prediction of esophageal squamous cell adenocarcinoma (ESCA), based on intestinal microbiota-related genes. Microbiome, transcriptome, and clinical phenotype data were downloaded from the TCMA, UCSC, and GEO databases. Intestinal microbiota-related prognostic genes were screened using differential- and univariate Cox regression analyses. The ESCA subtypes were explored based on prognostic genes, followed by infiltration analysis of immune cells in the different subtypes. Feature genes were identified using LASSO and multivariate regression analyses, and a prognostic model was established and validated. Univariate and multivariate Cox regression analyses were used to identify independent prognostic factors. The amount of Bacteroidetes and Proteobacteria differed significantly between the tumor and paracancerous samples. A total of 17 intestinal microbiota-related prognostic genes were also screened. We then identified three ESCA subtypes that exhibited significant differences in microbiota distribution and immune cell infiltration. Finally, we identified 10 gene signatures that were significantly associated with the survival status of ESCA patients. Furthermore, patients with no immunotherapy response presented a significantly higher risk score and prognostic risk. We constructed a prognostic model based on 10 characteristic genes related to intestinal microbiota, which could potentially predict the prognosis of ESCA patients.

KEYWORDS:

• Esophageal squamous cell adenocarcinoma
• intestinal microbiota
• Bacteroidetes and Proteobacteria
• prognostic model
• immune infiltration

Introduction

Esophageal squamous cell adenocarcinoma (ESCA) is a major esophageal cancer with a high incidence and mortality rate (Liang et al. 2017). It accounts for approximately 90% of the 456,000 yearly esophageal cancer cases (Abnet et al. 2018). The outcomes of ESCA are generally poor with five-year survival rates of only 13–18% (Napier et al. 2014). Therefore, in-depth scientific investigations of the pathogenesis of ESCA to identify medical interventions are needed (Holmes and Vaughan 2007). The human intestinal tract hosts various types of microbiota, especially in the oral cavity and gastrointestinal tract, which varies greatly and is host specific (Hooper et al. 2012). Alterations in the microbiota is an important cause of human diseases, such as inflammation and cancer (Di Pilato et al. 2016). Integrated intestinal microbiome and metabolome analyses has laid the foundation for revealing the detailed molecular mechanisms of the microbiome in human cancer (Yang et al. 2019). Thus, emphasizing the need to investigate the integration of microbiome data into the prevention, diagnosis, and therapy of ESCA is crucial.

The close association between intestinal microbiota and esophageal cancer risk was previously reported (Chen et al. 2021). The intestinal microbiota is involved in host physiological processes, including metabolism and immune defense, and the monitoring and intervention of intestinal microbiota may aid esophageal cancer screening and therapy, potentially improving the prognosis (Zhou et al. 2021). A previous study indicated that the bacterium, Fusobacterium nucleatum, contributes to the prognosis of esophageal cancer in humans (Yamamura et al. 2016). Additionally, obligate anaerobes in the gut are predominant in patients with esophageal tumors (Tanaka et al. 2012). Another study also identified characteristic intestinal flora in the stool samples of patients with advanced esophageal cancer (Motoori et al. 2017). Alteration of the intestinal microbiota contributes to the development of Barrett's esophagus and ESCA (Lv et al. 2019). Notably, certain differentially expressed genes (DEGs) have been linked with the alteration of intestinal microbiota in the progression of cancer (Sui et al. 2020). Long et al. demonstrated that the up-regulation of SLC52A3, which is associated with intestinal microbiota, can be used as a prognostic biomarker in esophageal cancer (Long et al. 2018). The potential mechanisms of microbiota in the progression of human cancer, through complicated interactions with the tumor microenvironment (TME) and signaling pathways, were previously reported (Qiu et al. 2021). A recent study demonstrated that the subtle interaction between microbiota and tumor cells is key to improving the clinical diagnosis and prognosis in human cancers (Sheetal et al. 2021). Differences in microbial composition was proven to be an effective clinical diagnostic marker for patients with ESCA (Deng et al. 2021). However, comprehensive research on the intestinal microbiota-associated biomarkers in the prognosis of ESCA is lacking.

In this study, the intestinal microbiota-related prognostic genes were identified based on data from the Cancer Microbiome Atlas (TCMA), University of California Santa Cruz (UCSC), and Gene Expression Omnibus (GEO) databases. This was followed by ESCA subtype identification and immune infiltration analysis on different subtypes. Subsequently, a microbiota-related feature gene-based prognostic model was established and the association of the prognostic model with survival prognosis was investigated. We hoped to develop a potential prognostic model based on prognostic feature genes of the intestinal microbiota for the clinical prognosis of ESCA.

Materials and methods

Microarray data and pre-processing

The RNA-seq data of 162 ESCA and 11 paracancerous samples were obtained from the UCSC Xene database (platform: Illumina HiSeq 2000 RNA Sequencing) (Tyner et al. 2017). The clinical information of these samples, including age, esophageal tumor location, neoplasm histologic grade, pathologic M, pathologic N, pathologic T, gender, tumor stage, and disease type, were also downloaded. The ESCA sample microbiome data (including 60 tumor and 22 paracancerous samples) were also acquired from the TCMA database (Dohlman et al. 2021). Among them, 50 samples were detected in both the Cancer Genome Atlas (TCGA) and TCMA databases. The microarray dataset, GSE53625, was acquired from the GEO database (Edgar et al. 2002) to externally validate the prognostic model. This dataset contains 179 ESCA tumor samples with clinical survival and prognostic information. Furthermore, GSE53625, which contains 132 ESCA samples, was used to validate genotyping based on gene clustering. Another immunotherapy dataset, GSE165252, which contains 71 esophageal adenocarcinoma samples (48 non-responders and 23 responders), was used to study the relationship between gene signatures and immunotherapy response. Dataset curation and study workflow are depicted in Figure 1.

Figure 1. Study flow chart.

Differential analysis of intestinal microbiota and genes

The relative abundance of intestinal microbiota at the phylum and genus levels of ESCA and normal samples was determined using TCMA. The t-test was used to compare the difference in intestinal microbiota abundance between patients with ESCA and normal subjects. The results were visualized using a box plot and the intestinal microbiota (P < 0.05) were considered differential intestinal microbiota involved in tumorigenesis. Moreover, the DEGs between tumor and normal samples were investigated with a false discovery rate of <0.05, and |log₂fold change| >1, using the limma package in R software (version: 3.6.1) (Ritchie et al. 2015).

Analysis of intestinal microbiota-related DEGs

The Pearson correlation coefficient (PCC) was calculated between differential microbiota and DEGs by the sample correspondence, followed by a significance analysis using the t-test. Subsequently, the Benjamini and Hochberg adjusted P value <0.05, and |PCC| >0.4 were considered as the cut-off values for the selection of microbiota-associated DEGs. Furthermore, the gene ontology (GO)-biological process (BP) and the Kyoto encyclopedia of genes and genomes (KEGG) pathway analyses were performed on intestinal microbiota-related DEGs using the DAVID (version: 6.8) database (Sherman and Lempicki 2009). A P < 0.05 was selected as the cut-off value for this enrichment analysis.

Investigation of prognostic DEGs associated with intestinal microbiota

A univariate Cox regression analysis (Wang et al. 2016) was performed using R software (version: 3.6.1) based on clinical survival information (overall survival [OS] and OS time) and the expression value of microbiota-associated DEGs to screen the prognostic DEGs related to intestinal microbiota. Statistical significance was set at P < 0.05.

Identification of ESCA subtypes

Based on the expression of the prognostic DEGs related to intestinal microbiota in each tumor sample, the ConsensusClusterPlus package (Zhang et al. 2020) in R software (version: 3.6.1) was used for the tumor subgroup investigation. Kaplan-Meier (KM) survival analysis was then carried out using the Survival package (version 2.41-1) to evaluate the survival correlations between different tumor subtypes, followed by statistical analysis and comparison. ESCA subtypes were further verified using an independent validation set, GSE53625, based on candidate gene expression.

TME analysis on different subtypes

Based on the gene expression patterns in ESCA tumor samples, the proportions of six immune cell types (including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and myeloid dendritic cells) in each subgroup were calculated using the online database Timer2.0 (Li et al. 2017). The Wilcoxon test in R software (version 3.6.1) was used to investigate the proportions of immune cells between different subtypes.

Prognostic model construction

Based on the least absolute shrinkage and selection operator (LASSO) Cox regression analysis in R software (version 3.6.1), the optimal gene set was screened for the presence of intestinal microbiota-related prognostic DEGs in the TCGA training dataset, in which genes were analyzed with 50-fold (nfold = 50 s) cross-validation analysis using the glment package (version: 4.0-2) (Friedman et al. 2009). Genes in the optimal gene set were identified as gene signatures. Then, based on their expression and sample prognostic information, multivariate Cox regression analysis was performed to calculate the prognostic coefficient of each gene. Finally, the following prognostic model was constructed: $\text{Prognostic score }(\text{PS})=\sum \text{Coe}{\text{f}}_{\text{genes}}\times \text{Ex}{\text{p}}_{\text{genes}}$where Coef_genes refers to the multivariate Cox regression coefficient of the target gene and Exp_genes refers to the expression levels of genes in the TCGA training dataset. The KM curve in R software (version 3.6.1) was used to evaluate the relationship between the expression of feature genes and the survival status of samples in the training set.

Prognostic model validation

Internal and external validation of the prognostic model was carried out based on the TCGA training set and the GSE53625 dataset, respectively. The median PS was calculated according to the prognostic coefficient of gene signatures and the samples were divided into high-risk and low-risk groups. The KM curve was used to analyze the survival prognosis of patients in the two groups. Furthermore, the GSE165252 dataset was used to explore the correlation between prognostic risk score and immunotherapy response in patients with esophageal adenocarcinoma. According to the above formula, the risk scores of samples in the GSE165252 dataset were calculated and compared between non-responders and responders, and the statistical significance was set at P < 0.05. Finally, to further observe the relationship between intestinal microbiota and different risk groups, a Sankey diagram was constructed to show the correlations between microbiota, gene signatures, and their prognostic effects.

Independent analysis of the prognostic model

Univariate and multivariate Cox regression analyses in the survival package (version 2.41-1) of R3.6.1 were used to screen independent prognostic factors. The log-rank P < 0.05, was selected as the threshold of significant correlation.

Results

Analyses of differential intestinal microbiota and DEGs between ESCA samples and controls

A total of three statistically significant differential intestinal microbiota were found at the phylum and genus levels each (Figure 2A-B). After eliminating the intestinal microbiota with low abundance, the Bacteroidetes and Proteobacteria with significant differences in abundance at phylum level were enrolled for the subsequent analysis to include more species. The abundance of Bacteroidetes was significantly higher and that of Proteobacteria was significantly lower in cancer tissues than in paracancerous tissues (Figure 2A). Furthermore, a total of 1438 DEGs, including 1063 upregulated and 375 downregulated genes, were identified between ESCA and paracancerous samples.

Figure 2. Box plot indicating the relative abundance of differential intestinal microbiota. A, relative abundance of differential intestinal microbiota between tumor samples and paracancerous samples at the phylum level. B, relative abundance of differential intestinal microbiota between tumor samples and paracancerous samples at the genus level. X- axis and Y-axis represent the different types of intestinal microbiota and expression levels, respectively. Blue and red bars represent the tumor and paracancerous samples, respectively.

Investigation of intestinal microbiota-related DEGs

A network of interactions between intestinal microbiota and DEGs was constructed (Figure 3A), with 522 phylum-gene interactions, 174 Bacteroidetes-associated genes, and 348 Proteobacteria-associated genes. Enrichment analysis was then performed on DEGs associated with Bacteroidetes and Proteobacteria. These DEGs were mainly enriched in pathways including the cell cycle (hsa04110) and mismatch repair (hsa03430) (Figure 3B). Additionally, these DEGs were mainly assembled in GO functions (Figure 3C), such as cell division (GO:0051301) and mitotic nuclear division (GO:0007088).

Figure 3. Identification and enrichment analysis of intestinal microbiota-related DEGs: A, interaction network of intestinal microbiota and DEGs. Circles represent DEGs with darker colors indicating more significant P values. The hexagon represents differential intestinal microbiota. The red and blue lines between the two nodes represent positive and negative correlation, respectively. The red node represents up-regulation, while the green node represents down-regulation. B, significantly enriched KEGG pathways. C, assembled GO-BP. Darker colors indicate more significant P values. The larger the node, higher was the number of genes enriched.

Investigation of prognostic genes

The univariate Cox regression analysis was used to screen prognostic genes associated with intestinal microbiota. The results indicated that 17 genes were identified as prognostic genes. Detailed information for all these prognostic genes is shown in Supplementary Figure 1.

Investigation of the ESCA subtypes

Based on the 17 prognostic genes identified, subtype analysis was performed on 151 tumor samples with prognostic information. The results identified three different subtypes, including cluster 1 (41 tumor samples), cluster 2 (57 tumor samples), and cluster 3 (53 tumor samples) (Figure 4A). Principal component analysis (PCA) (Figure 4B) and heatmap (Figure 4C) showed that all genes could be separated by different subtypes. The correlation analysis of prognosis among subtypes revealed a significant difference in survival status among the three subgroups (P = 0.0077). Cluster 3 had a better clinical survival prognosis, while clusters 2 and 1 had poorer clinical survival prognosis (Figure 4D). By validating the gene clusters in the GSE53625 dataset, we identified three subtypes (Figure 4E), of which, cluster 3 had a significantly favorable prognosis compared to clusters 1 and 2 (P < 0.0001, Figure 4F). Moreover, the distribution analysis of different intestinal microbiota in the three clusters showed that, compared with cluster 3, the relative abundance of Bacteroidetes was higher in clusters 1 and 2, suggesting that Bacteroidetes was a predictor of poor prognosis (Figure 4G). In addition, compared with clusters 1 and 2, the relative abundance of Proteobacteria was higher in cluster 3, indicating that Proteobacteria was a predictor of a favorable prognosis (Figure 4H).

Figure 4. Investigation of the three ESCA subtypes. A, consensus heatmap shows the three ESCA subtypes identified, including clusters 1, 2, and 3. B, PCA plot shows the distribution of samples in three subgroups. C, heat map shows the expressions of these 17 prognostic DEGs in different subgroups and stratified by different clinical information. D, the KM curve shows the difference in survival prognosis between different subtypes. X-axis and Y-axis represents survival time (months) and survival probability, respectively. E, genotyping based on gene clustering was validated in an independent validation dataset (GSE53625). F, survival status was compared between three clusters identified from GSE53625, using a KM curve. G, box plot shows Bacteroidetes distribution in the three subtypes. H, box plot shows Proteobacteria distribution in the three subtypes. I, box plots show the differences in the infiltrations of six immune cells between the three clusters.

TME analysis on different subtypes

The infiltration abundance of six types of immune cells in 151 cancer samples was calculated using the timer database. The results showed that, except for B cells, the other five types of immune cells, including CD4+ T cells, CD8+ T cells, neutrophil cells, macrophages, and myeloid dendritic cells, were significantly different between clusters 1 and 3, and their infiltration abundance was higher in cluster 3 than in cluster 1 and cluster 2 (Figure 4I).

Establishment and validation of the prognostic model

Based on the 17 prognostic genes obtained from the TCGA training set, 10 gene signatures were screened using LASSO Cox regression analysis (Table 1). Survival analysis of these 10 genes showed that high expression of ABRACL, ANGPTL3, DUSP2, HIST1H1E, HIST1H2BG, HSPA1B, HSPH1, and PTGES3 was associated with decreased survival, whereas high expression of MT1M and SCARA5 was associated with favorable survival (Figure 5A). A prognostic model was established based on the regression coefficients and expression values of these 10 genes. Samples were assigned to two different risk groups, based on the median risk score, in which patients in the high-risk group showed adverse survival compared to those in the low-risk group in both the training and validation sets (Figure 5B-E). By considering the risk score in an immunotherapy dataset, we also found significant higher prognostic risk scores in non-responder patients with ESCA, compared with responders in the GSE165252 dataset (Figure 5F), indicating that patients without immunotherapy response may present increased prognostic risks. Moreover, heatmap results showed that the expression of these 10 characteristic genes differed significantly between the high- and low-risk groups (Figure 6A). To observe the relationship between the different risk groups and ESCA subtypes, the distribution proportions of different clusters in the high-risk and low-risk groups were investigated using a bar chart (Figure 6B). The samples in clusters 1 and 2 were mainly assembled in the high-risk group, which indicated a worse prognosis. Furthermore, the positive correlations between microbiota and feature genes, as well as their prognostic roles in ESCA, were visualized using a Sankey diagram (Figure 6C). Seven gene signatures were positively correlated with Bacteroidetes, while the protective factors, including MT1M and SCAR5, were positively correlated with Proteobacteria, which further indicated that Bacteroidetes suggested a poor prognosis, while Proteobacteria suggested a favorable prognosis.

Figure 5. Construction and validation of the prognostic model. A, KM survival curves show the relationships between the expression levels of the 10 gene signatures and survival status of ESCA patients. B and D, distributions of the PS in the TCGA training dataset (B) and GSE53625 validation set (D). C and E, KM survival curves show the survival differences between high-risk and low-risk groups in the TCGA training dataset (C) and GSE53625 validation set (E). F, risk scores of samples in the immunotherapy dataset (GSE165252) were calculated and compared between non-responder and responder ESCA patients.

Figure 6. Correlation between gene signatures and the three clusters. A, heatmap shows the expression of the 10 gene signatures in the high- and low-risk groups. B, bar chart shows the distribution of samples of the three subtypes in high- and low-risk groups. C, Sankey diagram shows the relationships between microbiota, gene signature, and their prognostic effects.

Table 1. Information of the 10-gene-signature.

Gene name	HR (95% CI)	P-value	Coefficient
ABRACL	1.897 (1.27−2.835)	0.002	0.432
ANGPTL3	2.231 (1.091−4.563)	0.028	1.302
DUSP2	1.394 (1.013−1.918)	0.041	0.224
HIST1H1E	1.56 (1.162−2.094)	0.003	0.308
HIST1H2BG	1.409 (1.085−1.831)	0.01	0.133
HSPA1B	1.474 (1.13−1.922)	0.004	0.159
HSPH1	1.664 (1.185−2.338)	0.003	0.139
MT1M	0.774 (0.605−0.991)	0.042	−0.208
PTGES3	2.343 (1.269−4.324)	0.006	0.110
SCARA5	0.702 (0.493−1)	0.05	−0.221

Notes: HR, hazard ratio; CI, confidence interval; P < 0.05 was the cut-off value for significant difference.

Independent analysis for the prognostic model

Results of the univariate Cox regression analysis on the clinical factors and risk groups showed that pathologic M, pathologic N, tumor stage, and risk group factors were significantly associated with ESCA prognosis (Figure 7A). Multivariate Cox regression analysis was used to screen for independent prognostic factors. The results showed that the risk group was a unique independent prognostic factor for ESCA (Figure 7B).

Figure 7. Analysis of independent prognostic factors. Univariate (A) and multivariate (B) Cox regression analyses were used to identify independent prognostic factors according to different clinical characteristics.

Discussion

ESCA is one of the most aggressive squamous cell carcinomas and is highly prevalent in Asia (Ohashi et al. 2015). In recent years, tumor immunotherapy research has been extensively carried out in the clinic and the intestinal microbiota has also been found to be related to immune checkpoint inhibitors in the treatment of tumors (Gopalakrishnan et al. 2018). However, the exact clinical value of intestinal microbiota in the pathogenesis of ESCA remains unclear. In this study, based on the microbiome data in TCMA, the abundance of Bacteroidetes and Proteobacteria differed significantly between tumor samples and paracancerous tissues. A total of 522 microbiota-associated DEGs were investigated based on the dataset in UCSC Xene. These DEGs were mainly enriched in the mismatch repair pathway. The microbial environments of ESCA tissues are primarily composed of Firmicutes, Bacteroidetes, and Proteobacteria (Shao et al. 2019). Recent studies indicated that the intestinal microbiota of Bacteroidetes is related to the occurrence and progression of human cancer (Flemer et al. 2017). Moreover, Proteobacteria levels are decreased in ESCA patients with positive lymphatic metastasis (He et al. 2020). Liu et al. showed that Bacteroidetes had significantly higher relative abundances among N+ (positive lymph node) esophageal squamous cell carcinoma patients than in N- (negative lymph node) controls, whereas Proteobacteria were less abundant in N+ patients (Liu et al. 2018b). However, in the current study, all risk factors were positively correlated with Bacteroidetes, indicating that Bacteroidetes strongly associated with poor prognosis in ESCA. Additionally, changes in the intestinal microbiota are realized through the participation of related genes in specific functions or pathways. Deficiency in mismatch repair is closely associated with tumor progression (Cahill et al. 2007). Mismatch repair is commonly targeted with homologues of Dam and MutH, which are related to Proteobacteria (Fishel 2021). A genomic analysis of DNA repair genes in the alpha-Proteobacterium indicated that the change of genes in Proteobacterium leads to DNA repair dysfunction (Martins-Pinheiro et al. 2007). These results suggest the importance of Bacteroidetes and Proteobacteria in human cancers. Thus, based on the results of this study, we speculated that Bacteroidetes and Proteobacteria might participate in the progression of ESCA via association with the mismatch repair pathway.

A total of 17 prognostic DEGs were identified in the current study, followed by three ESCA subtypes. A poor prognosis was observed in both clusters 1 and 2, with a high abundance of Bacteroidetes. In contrast, a favorable prognosis was observed in cluster 3, with a high abundance of Proteobacteria. Moreover, compared with clusters 1 and 2, the relative infiltration abundance of immune cells, including CD4+ T cells and CD8+ T cells, was higher in cluster 3, possibly explaining the differences among subtypes. Tumor-infiltrating lymphocytes are associated with a survival benefit in several tumor types and in response to immunotherapy (Deschoolmeester et al. 2010). Patients with tumor cells with higher CD8+ T cell infiltration tend to have a better prognosis (Schumacher et al. 2001). The relationship between high CD8+ T cell abundance and favorable prognosis was explained by Jansen et al.: the immune system established a mechanism of ‘forward operating bases’ in some tumors, and patients with well supported immune cells in tumors were more likely to control the growth of tumor cells over a longer period (Jansen et al. 2019). The different recruitment results of CD8+ T lymphocytes into ESCA tissue lead to different survival outcomes in patients with different ESCA subtypes (Liu et al. 2015). Notably, Cho et al. indicated that CD4+ and CD8+ T cells strongly correlate with ESCA prognosis (Cho et al. 2003). A previous study indicated that CD4+ lymphocytes can be used to predict the survival rate of ESCA patients after clinical surgery (Chen et al. 2017). Interestingly, polysaccharide A from the capsule of Bacteroidetes induces clonal CD4+ T cell expansion, indicating a close relationship between CD4+ T cell infiltration and the biological function of intestinal microbiota (Johnson et al. 2015). Thus, based on results of the current study, the difference in infiltration abundance of CD4+ and CD8+ T cells might explain the survival difference among the three subtypes. Additionally, Bacteroidetes and Proteobacteria may affect ESCA prognosis by mediating T cell infiltration.

A prognostic model was constructed based on 10 genes, including ABRACL, ANGPTL3, and MT1M. The distribution analysis of the three clusters in the high-risk and low-risk groups showed that samples from clusters 1 and 2 were mainly assembled in the high-risk group, which led to a poor outcome. It has been shown that miR-145-5p functions in the proliferation, migration, and invasion of ESCA cells by stimulating ABRACL (Fan et al. 2021). ANGPTL3 is an indicator of prognosis in esophageal cancer patients (Zhu et al. 2015; Bao et al. 2018; Chen et al. 2019). This study found that ESCA patients with high ANGPTL3 expression had a better prognosis. Thus, the overexpression of ANGPTL3 might be an important clinical indicator for favorable ESCA outcomes. In terms of intestinal flora, we found that Proteobacteria may also be a potential biomarker for favorable ESCA prognosis. Although Mennini et al. found an increase in the abundance of Proteobacteria in esophagitis, enterocolitis syndrome, and allergic proctocolitis (Mennini et al. 2020), the relationship between Proteobacteria and the prognosis of esophageal cancer has not been confirmed. This study initially found that increased Proteobacteria and ANGPTL3 overexpression may be related to the prognosis of ESCA, but their underlying regulatory mechanisms still require further investigation. Regarding MT1M, Li et al. found that the inhibition of MT1M expression was a negative impression factor for esophageal cancer, further supporting the prognosis value of MT1M in esophageal cancer (Liu et al. 2018a; Li et al. 2021). All these studies suggest that ABRACL, ANGPTL3, and MT1M are valuable predictors of the prognosis of ESCA. Thus, we speculated that the prognosis model constructed in this study could predict the overall survival of ESCA patients.

However, the current study had some limitations, such as small sample size and lack of experimental analysis to verify the relationship between gene signatures and intestinal microbiota, owing to difficulty in collecting stool and solid tumor tissue samples from ESCA patients. Therefore, we will consider including a larger sample size into the next research plan to study the relationship between gene signatures and intestinal microbiota in improving ESCA prognosis. We will also investigate the molecular regulatory mechanisms through experimentation.

In conclusion, Bacteroidetes and Proteobacteria are not only associated with the prognosis of ESCA, but are also involved in the progression of ESCC via the mismatch repair pathway. Three ESCA subtypes were investigated and the difference in ESCA prognosis might be associated with CD4+ and CD8+ T cells. Furthermore, a prognostic model was constructed based on 10 gene signatures, which could be used to predict the overall survival of ESCA patients.

Acknowledgements

None.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available in UCSC Xene (http://xena.ucsc.edu/), TCMA (https://tcma.pratt.duke.edu/), and NCBI-GEO (https://www.ncbi.nlm.nih.gov/geo/).