Metabolic and Immunological Subtypes of Esophageal Cancer Reveal Potential Therapeutic Opportunities

Abstract

Background: Esophageal cancer has the sixth highest rate of cancer-associated deaths worldwide, with many patients displaying metastases and chemotherapy resistance. We sought to find subtypes to see if precision medicine could play a role in finding new potential targets and predicting responses to therapy. Since metabolism not only drives cancers but also serves as a readout, metabolism was examined as a key reporter for differences.

Methods: Unsupervised and supervised classification methods, including hierarchical clustering, partial least squares discriminant analysis, k-nearest neighbors, and machine learning techniques, were used to discover and display two major subgroups. Genes, pathways, gene ontologies, survival, and immune differences between the groups were further examined, along with biomarkers between the groups and against normal tissue.

Results: Esophageal cancer had two major unique metabolic profiles observed between the histological subtypes esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (EAC). The metabolic differences suggest that ESCC depends on glycolysis, whereas EAC relies more on oxidative metabolism, catabolism of glycolipids, the tricarboxylic acid (TCA) cycle, and the electron transport chain. We also noted a robust prognostic risk associated with COQ3 expression. In addition to the metabolic alterations, we noted significant alterations in key pathways regulating immunity, including alterations in cytokines and predicted immune infiltration. ESCC appears to have increased signature associated with dendritic cells, Th17, and CD8 T cells, the latter of which correlate with survival in ESCC. We bioinformatically observed that ESCC may be more responsive to checkpoint inhibitor therapy than EAC and postulate targets to enhance therapy further. Lastly, we highlight correlations between differentially expressed enzymes and the potential immune status.

Conclusion: Overall, these results highlight the extreme differences observed between the histological subtypes and may lead to novel biomarkers, therapeutic strategies, and differences in therapeutic response for targeting each esophageal cancer subtype.

Keywords: biomarkers; cancer metabolism; esophageal cancer; immunological subtypes; metabolic subtypes.

FIGURE 1

Partial least-squares discriminant analysis (PLS-DA) reveals histological differences for metabolic enzymes. (A) Hierarchical and PLS-DA clustering of enzyme mRNA levels between esophageal adenocarcinoma (EAC), esophageal squamous cell carcinoma (ESCC), and normal tissue (left) and just the histological subtypes (right). (B) Heatmap row Z-scores of log₂ mRNA expression. (C,D) PLS-DA containing scaled log₂ enzyme mRNA for 80% of each histological subtype that was randomly sent to training (C), while the remainder went to the testing dataset (D). (E) The receiver operating characteristic (ROC) curve of the training (upper left), testing (right), and the full dataset (bottom). (F) The confusion matrix with incorrect predicted classifications shown in red after tuning the sparse partial least-squares discriminant analysis (sPLS-DA). (H) The optimal combination of one component and 40 enzyme expression was further examined for their eigenvalues (influence) on component 1. (G) The balanced error rate for the first six components with 1–10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, or 300 enzyme mRNA expression(s) included in the tuned sPLS-DA. Large diamonds indicate the lowest error rate for the given component and number of enzymes. Components 3–6 could not decrease the balanced error rate further and are hidden behind the line of the second component. (I) The 40 enzymes in component 1 were examined for stability.

FIGURE 2

Key enzymatic markers between histological subtypes. (A) Random forest error rate with increasing number of trees when differentiating subtypes based on enzymes. (B) Random forest variable importance from 500 trees. (C) The fold change log₂ mRNA expression and -log₁₀ p-values generated from a Mann–Whitney U test for the optimal enzymes discovered from the sparse partial least-squares discriminant analysis (sPLS-DA) in Figure 1. (D) Machine modules with hyperparameter tuning trained on enzyme expression with the best logloss model being reported for each category, when applicable. Abbreviations are as follows: AUC, area under the curve; AUCPR, area under the precision-recall curve; RMSE, root mean square error; MSE, mean square error. (E) All of the models resulted in the same confusion matrix with the same holdout dataset. (F–K) Variable importance is seen for each of the top tuned models, including deep learning based on a feedforward neural net (F), generalized linear model (G), gradient boosting machine (H,I), and distributed random forest (J,K). Enzymes reported multiple times in the top 10 for importance in panels (F–H) and (J) had the expression plotted between the subtypes, with the number in parentheses indicating the number of occurrences (L). Error bars represent the standard error of the mean. Statistics were calculated through Mann–Whitney U test. NS, not significant, *p < 0.05, **p < 0.01, ***p < 0.001.

FIGURE 3

Metabolic pathways are altered between esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (EAC). Gene set enrichment analysis (GSEA) compared each patient against the average of normal tissues. Normalized enrichment scores (NESs) were plotted and hierarchically clustered for metabolic pathways, defined as ≥10 minimum enzymes with a minimum of 90% of enzymes, for (A) all Gene Ontology (GO) classifications and (B) metabolic pathways within the Kyoto Encyclopedia of Genes and Genomes (KEGG). (C) Each patient’s NES (top) and group average (bottom) for the tricarboxylic acid (TCA) cycle from GO biological process (BP) (left) and KEGG (right) (D) along with the heatmaps of the comprising genes within the pathway of GO (top) and KEGG (bottom). (E) Patient survival was plotted by separating the upper and lower quartiles of TCA NES through KEGG’s pathway for ESCC (left), EAC (middle), and both (right). (F) Patient survival was plotted by separating the upper and lower quartiles of COQ3 expression in ESCC (left), EAC (middle), and both (right). (G) Each patient’s NES (top) and group average (bottom) for glycerolipid enrichment from GO BP (left) and KEGG (right) (H) with the heatmaps of the comprising genes within the pathway of GO (top) and KEGG (bottom). (I) Patient survival was plotted by separating the upper and lower quartiles of glycerolipid NES through KEGG’s signature for ESCC (left), EAC (middle), and both (right). (J) Gene expression within each subtype. Pathway enrichment significance was calculated by a Student’s t-test, Mantel–Cox log-rank for survival, and gene expression was compared with a one-way ANOVA with a Bonferroni’s multiple comparison test. **p < 0.01, ***p < 0.001.

FIGURE 4

Esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (EAC) histological subtypes show distinct inflammatory responses. Normalized enrichment score (NES) was calculated for each patient’s tumor by comparing each patient’s expression to the normal adjacent tissue average expression (left) and as a group against the average of each subtype (right) for Gene Ontology biological processes (GO BPs). (A) The largest fold change of ESCC compared to EA was discovered to be acute phase response, (B) while the third was regulation of acute inflammatory response. (C) The same calculation was made for the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to show the NES Z-score (left) and average group enrichment (right). (D) Heatmap mRNA Z-scores of differentially expressed genes (Student’s t-test, p < 0.05) for GO term “cytokine activity” (GO: 0005125). The color bar on top indicates normal (gray), EAC (orange), and ESCC (blue). (E) mRNA expression levels for cytokines of interest relating to a Th17 signature. Differences in enrichment scores were calculated by the Student’s t-test and corrected with a Benjamini–Hochberg correction when examining multiple pathways in GO. Gene expression was compared with Kruskal–Wallis H test with Dunn’s multiple comparison test between all columns. **p < 0.01, ***p < 0.001.

FIGURE 5

Differences in the immune environment and predicted response to immune therapy between esophageal adenocarcinoma (EAC) and esophageal squamous cell carcinoma (ESCC) histological subtypes. Heatmap Z-scores, ranging from –2 (blue) to 2 (red), of differentially regulated (q < 0.05) immune infiltrates according to (A) CIBERSORT, (B) xCell, (C) and TIminer utilizing the dataset of Angelova et al. (2015). (D) mRNA Z-scores of differentially expressed (>1.5-fold, q < 0.05) genes of potential immunotherapeutic interest between EAC and ESCC. The color bars on top indicate esophageal tissues from normal (gray), adenocarcinomas (orange), and squamous cell carcinoma (blue). (E) Kaplan–Meier survival plots for patient survival based on CIBERSORT’s immune prediction. (F–I) Shown here are predicted therapy benefit (F), predicted therapy response (G), immune exclusion (H), and immune dysfunction prediction (I) according to the computational method of Tumor Immune Dysfunction and Exclusion (TIDE). Survival statistics utilized Mantel–Cox log-rank tests. Therapy prediction utilized Fisher’s exact test and a Student’s t-test for immune dysfunction, immune exclusion, and heatmaps, with the latter corrected by Benjamini–Hochberg correction for q-values.

FIGURE 6

Enzyme expression correlates with immune function prediction. (A–C) Legends (A) for Z-score of differentially expressed enzymes between subtypes (B) and Spearman’s rho for correlation between expression and immune function prediction from CIBERSORT and Tumor Immune Dysfunction and Exclusion (TIDE) utilizing both subtypes (C). (D) Kaplan–Meier survival curves for ALDH3B1. (E) Ranked Spearman’s rho correlation between enzyme mRNA and TIDE score for tumor immune dysfunction and exclusion utilizing both subtypes. Survival statistics utilized Mantel–Cox log-rank. *p < 0.05, **p < 0.01.