Data were obtained from COAD level 3 transcriptome sequencing data in the TCGA database. The protein assembly data and protein relative abundance details were obtained from the CPTAC database. The GMT file construction of specific pathway-related genes was obtained from the KEGG database. Bioinformatic analysis of the obtained data was performed using the R package (Limma, WGCNAC, ConsensusClusterPlus, and GSVA GSEAbase), GSEA 4.1 (http://www.gsea-msigdb.org/gsea/index.jsp), String online tool (https://cn.string-db.org/), and Cytoscape (https://cytoscape.org/index.html).

Human CRC tissue microarrays consisting of 79 pairs of tumors and adjacent normal tissues were obtained from Shanghai Outdo Biotech, China. The anti-EPHX2 antibody was purchased from ProteinTech (Ag1283). Intensity score was in accordance with the following criteria: 0, no appreciable staining; 1, weak staining; 2, moderate staining; and 3, strong staining. The percentage score was based on the percentage of CRC positive cells (0–100). The immunohistochemical (IHC) staining was scored by two independent pathologists, and the final IHC score was calculated by multiplying the staining intensity score and positive staining percentage score. IHC scores ranged from 0 to 12 and were graded (low, 0–4 points; medium, 5–8 points; and high, 9–12 points) as described previously (12). Paired t-test was used to compare H score between tumor tissues and adjacent normal tissues.

The normal human intestinal epithelial cell line (NCM460) and human colonic carcinoma cell lines (HCT116, HT29, HCT15, and LOVO) were purchased from the Culture Collection of Chinese Academy of Sciences (Shanghai, China). In short, the HCT116 and HT29 cell lines were maintained in McCoy’s 5A medium (iCell, China). NCM460 was maintained in DMEM (Gibco, USA), and HCT15 was maintained in RPMI 1640 medium (Gibco). LOVO was maintained in the F12k medium (iCell, China). All media were supplemented with 10% fetal bovine serum (Gibco), supplemented with 1% penicillin and streptomycin (Gibco). All cell lines were cultured in a humidified incubator containing 5% CO₂ at 37°C.

Total RNA was isolated using the TRIzol reagent (Invitrogen, USA), following the manufacturer’s protocol, and RNA purity was detected using a NanoDrop 2000 spectrometer (Thermo Fisher Scientific, USA). The quantitative real-time polymerase chain reaction (PCR) was performed using SuperReal PreMix Plus (Invitrogen) in a StepOnePlus Real-time PCR Detection System (Applied Biosystems, CA, USA). Relative gene expression was calculated using the 2^-△△Ct method. Human GAPDH was used as an endogenous control for mRNA expression in the analysis.

Two groups, the EPHX2-overexpressing lentivirus vector (EPHX2-OE group) group and negative control vector (EPHX2-NC group) group, were established for the experiments. Tissues with EPHX2 overexpression were co-transfected with siRNA-EPHX2 (EPHX2-OE-siRNA group) or siRNA-negative control (EPHX2-OE-siRNA_NC group). The EPHX2-overexpressing lentivirus vectors were designed and constructed by the OBiO Company (Shanghai, China). CRC cells (HTC116), infected with the lentivirus, were screened continuously for 2 weeks using 2 µg/mL puromycin. After the efficacy verification of EPHX2 overexpression, the lentivirus-infected colonic carcinoma cells were used in subsequent steps of the experiments.

The proteins were extracted from cells, boiled with RIPA buffer (Beyotime, Shanghai, China), loaded and separated on a 10% SDS-PAGE gel (EpiZyme, Shanghai, China), transferred onto polyvinylidene fluoride membranes, and then incubated with the following primary antibodies: anti-EPHX2 (ProteinTech, Catalog number:10833-1-AP, USA) and anti-β-actin (abclonal, Catalog number: AC026, China). Protein expression was detected using the ECL reagent. β-actin was used as control. Three biological replicates were performed in each experiment.

Annexin V–IF647/propidium iodide (PI) apoptosis assay was performed using an Annexin V–IF647/PI kit (Servicebio, China). The cells in each group were subsequently incubated with 5 μL of Annexin V–IF647 and 5 μL of PI solutions at 37°C for 10 min in the dark. The number of apoptotic cells was measured using a flow cytometer (BD Biosciences Co. Ltd., USA). Three biological replicates were performed in each experiment.

After lentiviral transfection, colonic carcinoma cells were seeded onto the basement membrane matrix (EC matrix, Chemicon, CA, USA) present in the insert of a 24-well culture plate. Thereafter, a complete culture medium was added to the lower chamber as a chemoattractant. After 48 h, the non-invading cells and EC matrix were gently removed with a cotton swab. The invasive cells located on the lower side of the chamber were fixed with 4% paraformaldehyde, stained with crystal violet, air-dried, and photographed with a microscope (Mshot, China). Three biological replicates were performed in each experiment.

Total RNA of HC116-GL120 and HCT116-H21064 were extracted using the TRIzol reagent (Invitrogen, CA, USA), following the manufacturer’s instructions. RNA sequencing was performed on an Illumina NovaSeq 6000 system of LC Sciences (USA). FastQC software (https://github.com/OpenGene/fastqc) were used to remove the reads that contained adaptor contamination, low quality bases and undetermined bases with default parameter. Then sequence quality was also verified using FastQC. We used HISAT2 (https://ccb.jhu.edu/software/hisat2) to map reads to the reference genome of Homo sapiens GRCh38.The mapped reads of each sample were assembled using StringTie (https://ccb.jhu.edu/software/stringtie) with default parameters. Then, all transcriptomes from all samples were merged to reconstruct a comprehensive transcriptome using gffcompare (https://github.com/gpertea/gffcompare/). After the final transcriptome was generated, StringTie and was used to estimate the expression levels of all transcripts. StringTie was used to perform expression level for mRNAs by calculating FPKM (FPKM = [total_exon_fragments/mapped_reads(millions) × exon_length(kB)]). The KEGG signal pathway enrichment analysis of FPKM data obtained by using GSEA4.1 software.

The cell coverslips were first washed with phosphate-buffered saline (PBS), subsequently fixed with 4% paraformaldehyde for 20 min, washed three times with PBS, and incubated with 60% isopropyl alcohol for 10 s, followed by staining with freshly prepared 60% Oil Red O solution (100% solution: 0.5 g of Oil Red O dissolved in 100 mL of isopropylene) for 30 min and washing with PBS. The samples were counterstained with hematoxylin. Three biological replicates were performed in each experiment. The cell climbing pieces were scanned using a Motic Easy Scan (Motic, China) and analyzed using the image-Pro plus 6.0 software (https://www.xrayscan.com/software-image-pro-plus/).

A cellular ROS assay kit (cat. no. ab186027; Abcam) was used to determine ROS levels. In a short, the transfected cells were plated (4 × 10⁴/100 µL per well) in a 96-well plate overnight. Subsequently, the cells were stained with ROS red stock solution for 30 min at room temperature. The absorbance was measured using a microplate reader, and ROS levels were determined by fluorescence increase at Ex/Em = 520/605 nm. Three biological replicates were performed in each experiment.

Weighted correlation network analysis (WGCNA) was used to analyze RNA sequencing data (level 3) from 388 tumor samples from TCGA-COAD to assess the correlation between gene modules and patient survival, tumor TNM stage, and distant metastasis ( Figure 1 ). The dark green module was mostly associated with the clinical characteristics of CRC through WGCNA, and the genes of this module may have a role in inhibiting the progression of CRC.

In bioinformatics prediction, the combination of mRNA and protein expression analyses of filtered data is beneficial to obtain more robust results. Therefore, the Colon Cancer Confirmatory Study (early release) protein assembly and protein relative abundance data (downloaded from the CPTAC database) were used to further analyze the genes in the dark green module. Genes with protein expression levels that were negatively correlated with TNM stage were used as the target gene set ( Figure 2 ). A gene set of nine genes, which may inhibit the progression of CRC, were obtained.

Cluster and survival analyses were used to further verify the role of target gene sets in CRC. Based on the mRNA expression levels of target genes set, the cluster analysis of 388 TCGA CRC samples was performed using the R package ConsensusClusterPlus. Consensus Clustering took a subsample from a set of microarray data and determined a cluster with a specified number of clusters (k). For each k, the paired consensus values (i.e., the proportion of the number of occurrences of two samples in the same subsample in the same cluster), were calculated and stored in a symmetric consensus matrix. Consensus matrices were summarized in several graphical presentations to enable the users to determine a reasonable number and membership of clusters. Cluster analysis was performed on 388 patients in the TCGA database according to the expression levels of target genes (Consensus Cumulative Distribution Function, CDF). The results showed that when k = 2, the consistency and clustering confidence was maximal. Furthermore, the survival analysis of patients with different cluster samples was performed. The Kaplan–Meier method was used to estimate the survival time distribution of different clusters, which was presented in the form of a survival curve to analyze the survival characteristics. The survival curve between both groups was tested using the log-rank test. The results showed that the survival rate of the cluster1 group was higher than that of the cluster2 group (p <0.05) ( Figure 3 ). This further illustrates the correlation between the target gene set and inhibition of CRC progression.

The core genes in the target gene set were predicted. First, GSEA enrichment analysis was performed on tumor samples based on cluster grouping. According to the normalized enrichment score (NES) and the false discovery rate (FDR), three optimal enrichment KEGG pathways were obtained (NES >2 and FDR <0.001). Thereafter, protein–protein interaction (PPI) analysis of target genes and genes of the optimal enrichment pathway was performed using the String online database and Minimum Required Interaction Score set to the highest confidence (0.9), and disconnected nodes were hidden in the network. The maximal clique centrality (MCC) method was used to calculate the key proteins in the PPI network (using CytoHubba plug-in of Cytoscape software). Of the nine target genes, EPHX2 was the first key protein to appear (top 12) ( Figure 4 ). First, the likely KEGG pathways of the target gene set were analyzed using GSEA. Subsequently, by PPI analysis, the first gene in the target gene set was selected and calculated as the core node. EPHX2 was considered to be the core gene involved in the inhibition of CRC progression.

To investigate the expression of EPHX2 in patients with CRC, the protein expression profiles of EPHX2 in human CRC specimens and adjacent normal tissues were examined. Consistent with previous bioinformatic analysis, IHC analysis showed that EPHX2 expression was significantly downregulated in human CRC specimens, compared with adjacent normal tissues ( Figure 5 ).

To verify the effect of EPHX2 on the progression of CRC at the cellular level, a CRC cell line (HCT116) stably transfected with EPHX2 overexpression (GL120) was constructed, and the transfection empty vector (H21064) was used as the negative control group (NC). The results showed that, compared with the NC group, overexpression of EPHX2 could inhibit the invasion of CRC cells and promote apoptosis. Moreover, EPHX2siRNA co-transfection could partially reverse the phenotype of cells ( Figure 6 ).

We further analyzed the molecular mechanism of EPHX2’s inhibition of the progression of CRC at the cellular level and performed RNA-seq analysis on HCT116 cells stably transfected with EPHX2, with the NC group as control and three cell lines in each group. GSEA enrichment analysis was performed for both groups. The results were screened using FDR <0.05, and 30 candidate KEGG pathways were obtained. The optimal enrichment signal pathway of the previously obtained target gene set in the clinical samples of CRC was used for intersection analysis ( Figure 7A ). Finally, peroxisome and fatty acid degradation signal pathways, which may be the most likely KEGG pathways for EPHX2 to act as the core gene in the target gene sets in the inhibition of the progression of CRC, were labelled ( Figure 7B ).

Hypoxia is an internal environment in almost all solid tumors. Hypoxia interacts with fatty acid metabolism and can regulate tumor immunity (17, 18). Therefore, the differences in the immune cell content and hypoxia signal pathway between the cluster1 and cluster2 groups were evaluated. The CIBERSORT R package was used to analyze TCGA-COAD RNA-seq data and perform differential analysis for the cluster1 and cluster2 groups. The KEGG pathway database was used to obtain genes related to the hypoxia signaling pathway and create GMT files ( Supplementary Material 1 ). The difference in hypoxia signal between the cluter1 and cluster2 groups was evaluated using the GSVA method. The expression of CD8T lymphocytes was upregulated in the cluster1 group, whereas the downregulation of M0 macrophage expression ( Figure 8A ) suggested that the target gene set may play a role in promoting the immunity of CRC. However, whether EPHX2 is directly involved in this biological effect remains to be further evaluated. No significant difference in hypoxia signal between cluster1 and cluster2 groups ( Figure 8B ) existed, which suggests that EPHX2 as a core gene inhibits colon cancer through peroxisome involvement in fatty acid degradation. Because it is not disturbed by the anoxic microenvironment of colon cancer, it cannot evaluate whether the hypoxic microenvironment of CRC is related to the expression of EPHX2.

Enhancing fatty acid degradation to generate ATP helps cancer cells to resist glucose deprivation or hypoxia, but the production of large amounts of ROS can lead to chronic oxidative stress. Under long-term oxidative stress, cancer cells are more likely to initiate the cell death process than normal cells (19). Lipid droplets (LDS) are highly organized spherical organelles and a repository of cellular fatty acids (20). The anoxic environment of solid tumors can promote the formation of lipid droplets by inducing PLIN2, which helps cells to use them for antioxidant defense after reoxygenation (21). Compared with normal colon tissue, CRC cells have several lipid droplets (22, 23). In previous studies, minimal difference was present in hypoxia signal in CRC samples grouped according to the target gene set. Therefore, the mechanism of EPHX2’s inhibition of the progression of CRC was preliminarily verified by detecting ROS content and oil red staining in both groups.

In HCT116 cells, the content of lipid droplets in the OE group was lower than that in the NC group (p <0.01) ( Figures 9A–C ). Additionally, higher levels of ROS were present in cells (p <0.01) ( Figure 9D ). This observation preliminarily confirmed that EPHX2 could promote the degradation of fatty acids and increase the release of ROS in HCT116 cells.