We downloaded the RNA sequencing data from TCGA database (https://portal.gdc.cancer.gov/) containing 80 colorectal MAC tissues and 51 normal colorectal tissues with clinical data.

For the RNA sequencing data we obtained, the R package limma (version 3.40.6) was used to conduct differential expression analysis to identify significant DEGs between MAC and normal tissues. A false discovery rate (FDR) < 0.05 and |log(fold change)| ≥ 2 were used to select DEGs.

Enrichment analyses of Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) were conducted by the R software package clusterProfiler (version 3.14.3) to assess the results of the gene set enrichment (17). FDR < 0.01 was taken as statistically significant. The KEGG rest API (https://www.kegg.jp/kegg/rest/keggapi.html) was utilized to gain the latest gene annotations of the KEGG pathway. Genes from the R package org.Hs.eg.db (version 3.1.0) were utilized for GO annotation.

The TCGA gene expression profile was used to remove genes with a standard deviation of zero in each sample, and goodSamplesGenes of the R software package, WGCNA, was applied to remove outliers and samples (16). WGCNA was then performed to build a scale-free co-expression network. Gene significance was determined as the mediating p-value for each gene (gene significance = lgP) in the linear regression linking gene expression and clinical characteristics. Module eigengenes were determined as the first principal component of each gene module, and the expression of module eigengenes was regarded to represent all genes in a particular module. Module membership was determined as the association between module eigengene and a given gene expression profile. Hub genes bound to MAC were obtained by computing the module membership and gene significance values. The cut-off criteria were |module membership| > 0.8 and |gene significance| > 0.1.

The least absolute shrinkage and selection operator (LASSO) regression model provides a new variable-screening algorithm that can effectively solve the collinearity problem. In our study, LASSO was performed to eliminate redundant factors and to identify the most significant survival-associated hub genes. After excluding patients with a follow-up period of less than 30 days, the R software package glmnet was applied to integrate gene expression data, survival status, and survival time, and the LASSO-Cox method was used for regression analysis. Moreover, 5-fold cross validation was set up to obtain the optimal model.

Patients were classified into two groups depending on the risk scores obtained by the Cox proportional hazards model, and the prognostic difference of the two groups was analyzed by the survfit function of the R software package (18). The Kaplan-Meier survival curve and log-rank test was utilized to assess the statistical significance of the prognosis among the high-risk and low-risk groups. Then, we performed receiver operating characteristic (ROC) analysis at 1, 3, and 5 years using Proc (version 1.17.0.1) of the R software, and evaluated the area under curve (AUC) and confidence intervals using the ci function of Proc.

To visualize the prediction of prognosis in patients with MAC, we constructed a nomogram depending on a couple of clinicopathological factors (survival time, survival status, age, sex, T stage, N stage, and M stage) and the 10 genes-based signature using the rms package in R. Harrell’s concordance index (C-index) and calibration curves, which could evaluate the consistency from the predicted survival probability and real observed rates, were conducted to assess the performance of the nomogram.

We classified the patients into a high-expression group (≥ 50%) and low-expression group (< 50%) depending on the expression levels of the hub genes. The Kaplan-Meier survival curve and ROC analyses were conducted, as previously described. Wilcoxon and Kruskal-Wallis tests were used to analyze the association between the hub gene expression and clinical characteristics.

GSEA software (version 3.0) was utilized to performed GSEA. The c2.cp.kegg.v7.4.symbols.gmt subset of the Molecular Signatures Database (http://www.gsea-msigdb.org/gsea/downloads.jsp) was downloaded to assess the related pathways and molecular mechanisms. FDR < 0.05 was regarded as statistically significant.

The RNA based stemness scores (RNAss) were calculated using the mRNA signature for each sample according to the algorithm of Malta et al. (19). The relationship between the expression of ENTR1 and infiltration of tumor-infiltrating immune cells was calculated using the R package IOBR (20) with the CIBERSORT (21) and ESTIMATE (22) methods. Correlation analysis was performed using Pearson’s correlation coefficient. Wilcoxon test was used to explore the relationship between gene expression and tumor-infiltrating immune cells.

The tumor tissues and paired normal tissues from the same patient were surgically harvested at Peking University Shougang Hospital, in accordance with institution-approved protocols. Tissues were collected from 13 patients with CRC, confirmed pathologically as MAC, who underwent radical surgery between January 1, 2020 to December 31, 2020. Informed consent was obtained from all participants involved in the study.

Tissues were fixed in formalin before embedding in paraffin. Paraffin-embedded tissues were sectioned, dewaxed, and dehydrated. Sections were 4 µm in thickness and deparaffinized in xylene, rehydrated in graded ethanol, and washed with TBS (1:20 dilution of 20x TBS, Solarbio, No. T1080) containing 0.3% Triton X-100 (T8200; Solarbio) (TBST) thrice. Sections were pretreated for antigen retrieval using citrate buffer (pH 6.0), cooled to room temperature (RT), and rinsed thrice with TBST. After blocking with 10% goat serum (ZSGB-BIO, ZLI-9021) for 1 h at RT, the tissue slides were incubated at 4°C for 8 h with ENTR1 antibody (Atlas Antibodies, A104721, 1:200 dilution). The sections were washed by TBST thrice after being rewarmed to RT, incubated with 3% H₂O₂ for 15 min, and then incubated with goat anti-rabbit IgG (Abcam, ab6721, 1:1000 dilution) at RT for 2 h. A DAB Staining Solution Kit (Gene Tech, Shanghai, GK600705) was used to stain the sections. The sections were counterstained with hematoxylin. Finally, all tissue slides were imaged and assessed using the IHC Profiler plugin (23), based on ImageJ bundled with Java 1.8.0 172 software (24). IHC Profiler used the average gray value (staining intensity) and positive area percentage (staining area) of positive cells as IHC measurement indices, and finally gave the sections three grades of scores as positive (≥2+), low positive (1+), and negative (0). The sections were then evaluated by a pathologist blinded to the nature of the samples, and manual correction was performed for all assessment results. Each tissue sample was replicated thrice. Negative controls without ENTR1 antibody were set for each section.

SPSS 26.0 (SPSS, Inc., Chicago, IL, USA) and GraphPad Prism 8 (GraphPad, Inc., CA, USA) were used for analysis. The difference between the two groups was calculated using the paired two-tailed Student’s t-test or Mann–Whitney–Wilcoxon test. Statistical significance was set at p < 0.05.

The analysis outline followed in this study was displayed in Figure 1 . We extracted the RNA sequencing data of MAC and normal colorectal samples from the TCGA-COAD and TCGA-READ datasets, and identified DEGs between MAC and normal tissues. Characteristics of the MAC and normal tissues were shown in Supplementary Table S1 . There were no significant differences between the two groups in age, sex and tissue location. Compared with the normal colorectal tissues, a total of 6,876 genes were upregulated and 3,455 genes were downregulated in MAC tissues (|fold change| ≥ 2, FDR < 0.05). Volcano plot of these genes is shown in Figure 2A . Figure 2B displays the top 20 genes up- and down-regulated in MAC.

Subsequently, we conducted KEGG and GO enrichment analyses of DEGs. The top five pathways involved in the DEGs, revealed by KEGG enrichment analysis, were neuroactive ligand-receptor interaction (FDR = 4.59e-16), cytokine-cytokine receptor interaction (FDR = 4.35e-11), protein digestion and absorption (FDR = 2.36e-8), viral protein interaction with cytokines and cytokine receptors (FDR = 5.77e-8), and complement and coagulation cascades (FDR = 8.52e-6) ( Figure 2C ). We also explored three main categories of GO enrichment: biological process (BP), cellular component (CC), and molecular function (MF). The top five pathways in the BP category were ion transport (FDR = 2.04e-26), regulation of signaling receptor activity (FDR = 2.04e-26), system development (FDR = 6.42e-24), humoral immune response (FDR = 6.24e-23), and animal organ development (FDR = 1.23e-22) ( Figure 2D ). In the CC category, the top five pathways identified were plasma membrane (FDR = 3.53e-37), intrinsic component of plasma membrane (FDR = 6.29e-34), integral component of plasma membrane (FDR = 3.44e-32), extracellular region (FDR = 8.69e-32), and extracellular matrix (FDR = 8.19e-23) ( Figure 2E ). In the MF category, the top five pathways involved in the DEGs were receptor ligand activity (FDR = 2.81e-26), receptor regulator activity (FDR = 1.13e-23), antigen binding (FDR = 7.95e-22), inorganic molecular entity transmembrane transporter activity (FDR = 2.13e-20), and channel activity (FDR = 1.08e-17) ( Figure 2F ).

To determine the crucial modules most correlated with MAC, we performed WGCNA on the DEGs in MAC and normal tissues. Taking a soft-thresholding power of 5 (scale-free R² = 0.87), we identified 32 gene modules ( Figures 3A–D ). From the correlation heatmap of the module eigengenes and MAC, we noticed that the paleturquoise module had the most positively correlation with MAC (p = 5.2e-10; Figure 3E ). The scatter plot of module membership in the paleturquoise module and the gene significance for MAC indicated a high positive correlation (correlation coefficient = 0.90, p = 2.1e-145, Figure 3F ). By setting up the cut-off criteria (|module membership| > 0.8 and |gene significance| > 0.1), we identified 71 hub genes from the 391 genes in the paleturquoise module ( Supplementary Table S2 ).

The functions and signaling pathways of the hub genes were analyzed using KEGG pathway analysis and GO functional enrichment analysis. According to KEGG pathway analysis, these hub genes were enriched in the cell cycle (FDR = 1.35e-4) and DNA replication (FDR = 1.46e-4) pathways, when FDR < 0.01 was considered statistically significant ( Figure 4A ). GO functional enrichment analysis showed that the highest enriched GO terms in the BP, CC, and MF categories were cell cycle (FDR = 8.15e-8, Figure 4B ), nuclear lumen (FDR = 1.67e-8, Figure 4C ), and anaphase-promoting complex binding (FDR = 5.27e-4, Figure 4D ), respectively. In brief, these results indicated that hub genes were mainly involved in the cell cycle process.

To construct a risk score assessment for predicting the OS of patients with MAC, LASSO regression analysis was used based on the aforementioned 71 hub genes. Eventually, a prognostic signature was constructed containing 10 hub genes: TROAP, C19orf48, CCNF, ZMYND19, RUVBL1, PAFAH1B3, ENTR1, NTMT1, RANGAP1, and IFRD2 ( Figures 5A, B ). These ten hub genes, shown in Figure 5C , were closely related to each other in the expression of mRNA, which indicated that they may be functionally related. We used Cox proportional hazards regression analysis of the hub genes to develop a prognostic signature. After excluding patients with follow-up periods of less than 30 days, patients with MAC were divided into low-and high-risk groups, according to the risk scores. Clinical characteristics of the high-risk and low-risk groups were shown in Supplementary Table S3 . There were no significant differences in clinical characteristics between the two groups, except for more patients in the high-risk group with stage N2 (p = 0.025). The distribution of risk scores was analyzed, and the relationship between the prognostic signature and the expression of the 10 hub genes was observed. A marked decrease in the OS of patients with MAC was observed as the risk score increased ( Figure 5D ). The Kaplan-Meier survival curve revealed that patients in the high-risk group had significantly worse OS than those in the low-risk group (p < 0.0001, hazard ratio = 21.40, 95% CI 2.82–162.61, Figure 5E ). Hereafter, we explored the performance of the ten hub genes-based signature for prognosis prediction. The AUCs for 1-, 3-, and 5-year OS were 0.87 (95% CI, 1.00–0.71), 0.92 (95% CI, 1.00–0.81), and 0.89 (95% CI, 1.00–0.72), respectively ( Figure 5F ). ROC curve analysis showed that the signature of the 10 genes had good prognostic prediction ability for patients with MAC.

To further visualize the OS probability of patients with MAC, we constructed a nomogram with age, sex, T stage, N stage, M stage, and the 10 hub genes-based prognostic signature. Depending on the nomogram of total points-to-outcome, patients with higher total points were estimated to have worse 1-, 3-, and 5-year OS probabilities ( Figure 5G ). The C-index of the nomogram was 0.80 (95% CI 0.72–0.89, p < 0.001). Moreover, the nomogram calibration curves showed promising performance in predicting 1- and 3-year OS probabilities ( Figure 5H ).

When we investigated the prognostic value of the 10 hub genes, only ENTR1 was closely associated with OS. Depending on the expression level of ENTR1, the patients were classified into either a low-expression group (< 50%), or a high-expression group (≥ 50%). Clinical characteristics of the high-risk and low-risk groups were shown in Supplementary Table S4 . The high-expression group had worse OS than the low-expression group (p = 0.016, hazard ratio (HR) = 3.74, 95% confidence interval (CI), 1.19–11.75, Figure 6A ). The AUCs of ENTR1 for 1-, 3- and 5-year OS were 0.69 (95% CI, 0.90–0.48), 0.64 (95% CI, 0.87–0.41) and 0.71 (95% CI, 0.98–0.45), respectively ( Figure 6B ). We further explored the association between ENTR1 expression levels and clinical features ( Figures 6C-H ). ENTR1 was upregulated in patients with MAC (p < 0.0001). However, except that ENTR1 was more highly expressed in stage N0 than in N1 (p = 0.01), no differences were found in ENTR1 expression at different T stages, N stages, M stages, TNM stages, or tumor locations.

We conducted GSEA to investigate the potential functions of ENTR1 in MAC. GSEA revealed that the ENTR1 expression was associated with pyrimidine metabolism (FDR = 0.026, Figure 6I ).

To evaluate the correlation between ENTR1 expression and cell stemness of MAC, the RNAss of the MAC samples were calculated using the mRNA signature. The results indicated that ENTR1 expression was significantly positively correlated with the cell stemness of MAC (correlation coefficient = 0.44, p < 0.0001, Figure 7A ).

Subsequently, ESTIMATE was used to assess immune infiltration in the MAC samples. ESTIMATE analysis suggested that ENTR1 expression was negatively correlated with stromal scores (correlation coefficient = -0.24, p = 0.03, Figure 7B ). Immune and ESTIMATE scores displayed a downward trend with increased ENTR1 expression. However, these differences were not statistically significant ( Figures 7C, D ).

The association of ENTR1 with tumor-infiltrating immune cells in the high- and low-expression groups was analyzed using the CIBERSORT algorithm. Analysis of the profile of 22 types of tumor-infiltrating immune cells demonstrated that the number of CD8+ T cells and T follicular helper cells was significantly higher in the high-expression group (p = 0.02, Figure 7E ). Furthermore, we explored the relationship of ENTR1 expression with CD8+ T and T follicular helper cells. The results revealed that ENTR1 expression was positively correlated with CD8+ T cells (correlation coefficient = 0.27, p = 0.01, Figure 7F ), whereas no significant correlation was found with T follicular helper cells (correlation coefficient = 0.11, p = 0.31, Figure 7G ).

IHC was used to determine whether ENTR1 expression was higher in the MAC group (n = 13) than in the normal group (n = 13). Ultimately, both groups could be divided into negative (-), low positive (+), and positive (++/+++) groups, based on ENTR1 expression levels. Figure 8A showed the extracellular mucus of MAC, and ENTR1 was stained brown and expressed on the cytoplasm and membrane. The overall positive rate of ENTR1 expression in the MAC group was 92.3% (12/13), compared to 69.2% (9/13) in the normal group. The results of the Mann-Whitney-Wilcoxon test showed that there were significant differences in the expression levels of ENTR1 between the two groups (Z = 2.75, p = 0.01). The positive area of ENTR1 in the MAC group was significantly higher than that in the normal group (p = 0.038, Figure 8B ).