A novel fatty acid metabolism-related signature identifies MUC4 as a novel therapy target for esophageal squamous cell carcinoma

The clinical data and gene expression profile of 259 ESCC samples were obtained from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. Of all the ESCC samples, the RNA-seq data (FPKM format), regarding 80 ESCC samples, was retrieved for TCGA database. Meanwhile, GSE53625↗ incorporating gene expression profiles based on GPL18109↗ platform was downloaded from GEO database. The ‘caret’ package in R software randomly allocated all the ESCC samples (n = 259) into the training and testing cohorts in a 1:1 ratio. Subsequently, 130 ESCC samples were used as a training cohort, and 129 ESCC samples employing as the testing cohort.

We acquired the FA metabolism-related gene sets (KEGG fatty acid metabolism pathways, Reactome fatty acid metabolism genes, and Hallmark fatty acid metabolism genes^20–23), and further analyzed for FMGs. FMGs were extracted from the Molecular Signature Database v7.4 (MSigDB).

The consensus clustering arithmetic is used to assess the cancer features differences within a specified data set²⁴. To identify cluster subtypes of ESCC for optimization classification purposes in FMGs, consensus clustering was performed in all ESCC samples through the “ConsensusClusterPlus” package by the K-means algorithm, which comprised consensus matrix, cumulative distribution function (CDF), and delta area plots. Stable consensus clustering classification was presented with the specific parameters (re-samplings, 50; pltem, 0.8; pearson correlation distances, 1). The survival analysis between different cluster subtypes were performed through the Kaplan–Meier survival curves method, and chi-squared test was applied to explore the survival difference comparisons among different subtypes.

Based on the genes expression profiles and prognostic information intersected in training cohort, the FMGs that exhibited significance (p < 0.05) in univariate Cox regression analysis for survival were screen as the prognostic-related FMGs and then included into the least absolute shrinkage and selection operator (LASSO) analysis with tenfold cross-validation. Subsequently, LASSO analysis was applied for selecting best candidate FMGs into the risk model by using the “glmnet” R package. The specific risk score of each ESCC patients was analyzed based on the FMGs expression levels and corresponding regression coefficient, and the detailed formula was as follows: risk score = ∑1i(gene Expression × gene coefficient). Subgroups including high-risk and low-risk groups were determined based on their respective median risk score. The same algorithm analysis was used for the testing cohort and all ESCC cohorts. The relationship between the risk score and the clinicopathological features (age, gender, T stage, N stage, and TNM stage) was explored. Additionally, the "survminer" and "timeROC" R packages were utilize to estimate overall survival (OS) difference and the predictive accuracy of the risk model, respectively.

Clinicopathological variables and risk score were performed univariate and multivariate Cox regression analyses, and the “forestplot” R package was utilized to screen the independent prognostic factors. A nomogram was generated based on age, gender, T stage, N stage, TNM stage, and the risk score with “regplot”, “survival” and “rms” R packages. In this nomogram arithmetic system, each clinicopathological variable was serve as a score, and overall scores of all variables in each sample were calculated and exhibited. ROC curve was analyzed to evaluate the discrimination reliability of the nomogram.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp/kegg/kegg1.html↗) analyses of differentially expressed FMGs among different subtypes were conducted through the “org.Hs.eg.db” and “clusterProfiler” R packages employing all ESCC samples data. In addition, to further search inherent mechanisms associated with FMGs among the two risk groups, the gene set enrichment analysis (GSEA) was performed to distinguish hallmark pathways associated with FMGs using the “limma” R package, with c5.go.symbols.gmt and c2.cp.kegg.v7.0.symbols.gmt as reference gene sets. p < 0.05 and q false-discovery rate (FDR) value less than 0.25 were considered valid.

The composition of the 21 kinds of infiltrating immune cells and immune function were performed in ESCC samples using “CIBERSORT” R package, and the heatmap was visualized. The ESTIMATE algorithm identifies tumor microenvironment features associated with immune, stromal, and ESTIMATE score for each patient using “estimate” R package. The associations among the cluster subtypes, risk score, and immune checkpoint expression were detected using Spearman’s correlations. Kruskal Wallis test was applied to perform the differences analysis between cluster subtypes and risk groups.

Normal human esophageal epithelial cells (HEEC) and human ESCC cell lines ECA109 and TE1 were obtained by Scientific Research Center of the Fourth Hospital of Hebei Medical University. These cells were cultured with RPMI-1640 containing 10% Fetal Bovine Serum (FBS), at 37 °C with 5% carbon dioxide incubator. Total RNA of cells was isolated with TRIzol reagent (Thermo Fisher Scientific) and reverse-transcribed to cDNA with PrimeScript Strand cDNA Synthesis Kit (Thermo Fisher Scientific). qRT-PCR analysis was performed with MonAmp™ SYBR® Green qPCR Mix, and the risk model-relevant genes was analyzed with 2^−ΔΔCT method. The primer sequences of these gene and GAPDH are displayed in Table 1.

Human ESCC cell lines and HEEC cells, were obtained from the Research Center of the Fourth Hospital of Hebei Medical University (Shijiazhuang, China). For the current study, all cells were cultured in RPM1-1640 medium (Gibco, USA) containing 10% fetal bovine serum (FBS) (Invitrogen, USA) and 100 U/ml penicillin (Invitrogen, USA) and 100 μg/ml streptomycin (Invitrogen, USA) in the 37 °C and 5% CO₂ cell incubator.

The ECA109 and TE1 cells were transfected with human MUC4 small interfering RNA (si-MUC4) and negative control shRNA (si-NC) using Lipofectamine 2000 (Invitrogen, USA) in accordance with the instructions of manufacturer. After 24 h, ECA109 and TE1 cells were harvested for subsequent functional assays.

Cell proliferation was analyzed using Cell Counting Kit-8 (Med Chem Express Princeton, USA). ESCC cell lines were harvested and seeded into 96-well plates with five replicate wells for each group per experiment. Different groups of ESCC cells lines were incubated with 10 µl CCK-8 at 24 h, 48 h and 72 h. Then, absorbance was measured following the manufacturer’s instructions, the optical density (OD) value was determined at a wavelength of 450 nm by a Multiskan microplate (ELx800; Bio-Tek, Winooski, VT, USA). At least three times dependent experiments were repeated.

The ECA109 and TE1 cells were plated in 6-well plates. When these ESCC cells grew to confluence > 80%, a straight cell-free zone was scratched using a 200-µl pipette tip. These cells were cultured in 2 ml of serum-free medium. The scratched areas were examined at 0 h, 24 h by using the microscope and photographed.

Transwells assays was used for testing the invasive ability of ESCC cells. ECA109 and TE1 cells were diluted with serum-free medium and 200 µl RPMI-1640 medium were seeded into the Matrigel-coated (Corning Inc. USA) upper chamber. 600 µl RPMI-1640 medium containing 20% FBS was added into the lower chamber. After incubation for 24 h with 37 °C constant temperature, cells on the upper side of the chamber were wiped off with cotton swabs, and the cells were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet for 15 min. Representative images were captured under fluorescence microscope (Nikon Ti2, Japan).

In this study, we identified prognostic-related FMGs among 259 ESCC samples in the all cohorts. The univariate Cox regression analysis identified 128 FMGs associated with prognosis. Moreover, the interaction correlations and degree of risk for prognosis of top thirty prognosis-related FMGs were illustrated by a network diagram (Fig. 1). The results indicated that PTGS2, ALOX12, IDI1, PECR, EPHX2, ALOXE3, PLA2G4A, ALOX12B, CYP2C9, MORC2, ACAT1, and HSD17B8 were favorable factors for OS, others were risk factors for OS.

Based on the screened prognosis-related FMGs expression data, the consensus matrix k value of 2, 259 ESCC patients were clustered into two subtypes (cluster 1 and cluster 2) with optimal consistency and stability (Fig. 2A,B). Kaplan–Meier survival analysis revealed that cluster 1 patients had a significantly poorer prognosis that in in cluster 2 patients (P = 0.0016; Fig. 2C). A heatmap depicted the prominent difference of prognosis-related FMGs expression and clinicopathological variables in these two clusters (Fig. 2D). These findings indicated that cluster subtypes could distinguish survival risk stratification of ESCC.

GO and KEGG functional analysis were conducted on the differentially expressed genes, with the enriched GO and KEGG terms in each category visualizing in Circos plots (Fig. 3A,B). The findings revealed that differentially expressed genes were remarkably associated with the FA metabolism related-regulation, such as FA degradation, arachidonic acid metabolism, FA elongation, beta-Alanine metabolism, and tryptophan metabolism (Table 2). The ssGSEA analysis revealed that the activated dendritic cells, CD56dim natural killer cells, gramma delta T cells, natural killer T cells, natural killer cells, and regulatory T cells had significantly higher in cluster 1 than in cluster 2 (Fig. 4A). Then, we analyzed the association with immune functions and found that APC co-inhibition, CCR, parainflammation, T cell co-inhibition, and type II IFN response had significantly different immune function scores in two subtypes (Fig. 4B). As shown in Fig. 4C–E, patients in cluster 1 had the higher stromal, immune score and ESTIMATE score than in cluster 2.

The prognostic-related FMGs (p < 0.05) were screened using univariate Cox regression analysis and these genes were included into the LASSO regression model for screening the candidate risk model-related genes. A eight- FMGs risk model was finally established based on the appropriate λ value. The risk score calculating formula in the training cohort, validation cohort and all cohorts was as follows: risk score = (0.0422) × FZD10 + (0.0034) × TACSTD2 + (− 0.0941) × MUC4 + (0.0313) × PDLIM1 + (− 0.1334) × PRSS12 + (0.1114) × BAALC + (0.2159) × DNAJA2 + (− 0.0464) × ALOX12B.

ESCC patients of training cohort, testing cohort and all cohorts were grouped into two risk groups depending on median value of risk scores. As shown in Fig. 5A–C, high-risk scores patients accounted for the higher percentage of patient deaths. A heatmap depicted the differential expression levels of genes among the different risk scores (Fig. 5D–F). The results of risk scores and survival rates in the training cohort were consistent with the testing cohort and all cohorts. The results revealed that high risk scores exhibited the unfavorable effect for prognosis. Kaplan–Meier survival analysis showed that patients in the low-risk survived longer than those in the high-risk group in the training cohort (p < 0.001), testing cohort (p = 0.001) and all cohorts (p < 0.001). The ROC curve results suggested that the AUCs for prognostic prediction were 0.748, 0.660, and 0.701, respectively (Fig. 6A–C).

We performed the univariate and multivariate prognostic analyses to investigate the predictive ability of FA metabolism risk score for survival. Combining FA metabolism risk score with clinicopathological features in all cohorts, the results showed that the age (p < 0.001), TNM stage (p < 0.001), and risk score (p < 0.001) were proved as independent prognostic factors (Fig. 7A,B). Figure 7C showed calibration curves for predicting ESCC patients at OS of 1, 3, and 5 years, indicating that nomogram had reliably prediction capability. Considering the clear and visual characteristics of nomogram, so we visualized the risk signature using a clinicopathologic nomogram. A nomogram integrating age, gender, T stage, N stage, and risk score was analyzed to determine the relationship between these factors and prognosis (Fig. 7D). Further analysis indicated that the nomogram (AUC = 0.769) has a better predictive performance than a single prognostic indicator such age (AUC = 0.564), TNM stage (AUC = 0.622), or the prognostic risk scoring model (AUC = 0.725; Fig. 7E).

To further explore the potential molecular mechanism of the risk model, the transcript message in ESCC patients in the all cohorts were analyzed using GSEA. The GO analysis further revealed that glycoside metabolic process, reverse cholesterol transport, positive regulation of mitochondrial translation, mitochondrial ATP synthesis coupled proton transport, and antigen processing and presentation of endogenous peptide antigen were enriched in the high-risk group (Fig. 8A), and immune-related activities in low-risk group were shown in Fig. 8B. Interestingly, the activity of metabolic pathways such as phenylamine metabolism, glutathione metabolism, and tyrosine metabolism, were fundamentally enriched in the high-risk group, according to KEGG analysis (Fig. 8C). Additionally, JAK STAT signaling pathway, FC gamma R mediated phagocytosis, T cell receptor signaling pathway, cytokine-cytokine receptor interactions were enriched in the low-risk group (Fig. 8D). Taken together, our results revealed that the FA metabolism reprogramming condition can act on the immune activity through internal regulatory mechanisms.

To examine the ESCC tumor microenvironment among the two risk groups patients, we explored the distribution of 21 immune cell types in different risk groups in all cohorts using CIBERSORT. The proportions of immune cells in the high-risk and low-risk groups were shown in Fig. 9A. The results indicated that the high-risk group had lower levels of naïve B cell, activated mast cells, plasma cells, and CD4 T cell memory resting, compared to the low-risk group. The StromalScore, ImmuneScore, and ESTIMATEScore were computed using the “ESTIMATE” package, revealing that ImmuneScore, StromalScore, and ESTIMATEScore were significantly higher in high-risk group than in low-risk group (p < 0:05) (Fig. 9B–D). Additionally, we evaluated the relationship between the tumor-infiltrating immune cells and prognosis, observing that increased plasma cells, monocytes, mast cells activated, and dendritic cells activated were closely related to good prognosis (Fig. 9E–H). High infiltration of T cells memory activated, M0 macrophages, M1 macrophages, and M2 macrophages showed unfavourable effects for survival (Fig. 9I–L).

The expression of 39 immune checkpoints in the high-risk and low-risk groups were also compared. The results revealed that most immune checkpoints were significantly expressed in different risk groups (Fig. 10A). Accordingly, the relationship between the immune checkpoints and risk score was performed, which suggested that the expression level of CD200, KDR, KIR2DL1, LAG3, NRP1, TGFBR1 and TNFSF15 were significantly associated with risk score (Fig. 10B–H).

The eight-FAMs risk model included the following genes: FZD10, TACSTD2, MUC4, PDLIM1, PRSS12, BAALC, DNAJA2 and ALOX12B. To further verify the expression level of the above mentioned FMGs, we detected the expression levels of these genes in HEEC, EAC109 and TE1 cells (Fig. 11A–H). The results of qRT-PCR revealed that the MUC4 and BAALC was upregulated, whereas the expression of TACSTD2, FZD10, PDLIM1, and ALOX12B was downregulated in ESCC cells than that in HEEC. However, we found no significant differences in expression levels of PRSS12 or DNAJA2 between HEEC and ESCC cell lines.

To explore the oncogenic roles of MUC4 in ESCC, we selected ECA109 and TE1 cells with low MUC4 expression by using si-MUC4 transfection subsequent experiments. MUC4 expression was knocked down in ECA109 and TE1 cells, and was demonstrated by qRT-PCR analysis (Fig. 12A). To evaluate the role of MUC4 on the proliferation, CCK-8 assay was used to validate and results showed the viability of MUC4 groups with knockdown were obviously lower than that of NC groups at 24, 48 and 72 h, respectively (Fig. 12B,C). Moreover, wound-healing assay suggested that after culture for 24 h, the knockdown of MUC4 expression significantly reduced the migration ability of ECA109 and TE1 cell lines (Fig. 12D,E). In addition, knockdown of MUC4 markedly attenuated the invasion ability in ECA109 and TE1 cell lines via transwell assay (Fig. 12F,G). Our findings showed that knockdown of MUC4 can suppresses ESCC cells proliferation, migration and invasion.

Data is provided within the manuscript or available from the corresponding author on reasonable request.