A retrospective analysis was performed on samples from patients with colorectal tumors treated at the Department of Gastroenterology, The Seventh Medical Center of Chinese PLA General Hospital (Beijing, China) between October 1, 2016 and October 31, 2016. Patient selection was based on the following inclusion criteria: i) Pathological confirmation of colorectal adenoma or adenocarcinoma; ii) availability of fresh-frozen tumor tissue with paired adjacent normal tissue (within 5 cm of the tumor margin); and iii) complete clinical and pathological data. The exclusion criteria were: i) Prior neoadjuvant chemotherapy or radiotherapy; and ii) diagnosis of synchronous multiple cancers or hereditary CRC syndromes. Based on these criteria, a total of 10 patients were included. The cohort comprised 2 cases of tubular adenoma, 3 of villous adenoma, 3 with early-stage cancer and 2 with advanced tumors. In addition, 4 adjacent normal tissue samples from these 10 patients were collected as controls.

The normal human colon cell line NCM460 (cat. no. IM-H445; cell batch, IM-H445031603) was purchased from the Xiamen Yimo Biotechnology Co, Ltd. These cells were cultured in Dulbecco's Modified Eagle's Medium (DMEM; Procell Life Science & Technology Co., Ltd.) and supplemented with 10% (v/v) heat-inactivated fetal bovine serum (FBS; Shanghai ExCell Biology, Inc.), 1% (v/v) penicillin-streptomycin and 1% (v/v) L-glutamine. All cell lines were maintained in a humidified incubator at 37˚C with 5% CO₂. All cell lines tested negative for mycoplasma contamination.

To construct shRNA-expressing vectors, annealed complementary shRNA oligonucleotides were ligated into the PiggyBac-H1-2O2-shRNA-CopGFP-puro vector, which was constructed from the PiggyBac-CMV-CNR (human)-EF1a-CopGFP-T2A-Puro plasmid (cat. no. P31653; MiaoLingBio) by replacing the CMV promoter with the H1-2O2 promoter. The shRNA sequence was 5'-GCAAGAAGTACAAAGTGGAGT-3' and the sequence for the negative control group was 5'-TTCTCCGAACGTGTCACGT-3'. NCM460 cells were seeded into 6-well plates and transfected with a total of 4 µg of plasmid DNA using Lipofectamine 3000 reagent (cat. no. L3000008; Invitrogen; Thermo Fisher Scientific, Inc.). The plasmid mixture consisted of the shRNA-expressing PiggyBac vector and Super PiggyBac Transposase (cat. no. P0179; MiaoLingBio) at a mass ratio of 2.5:1. The DNA-lipid complexes were incubated at room temperature for 15 min before being added to the cells. After 24 h of incubation, the cells were selected with 1 µg/ml puromycin for 3 days, after which almost all surviving cells exhibited green fluorescence. Single clones were then obtained by serial dilution and maintained under a lower dose of puromycin (0.1 µg/ml) for 1 to 2 weeks to establish stable cell lines. Protein knockdown was finally confirmed by western blotting.

RNA was extracted from clinical samples using RNAiso Plus (cat. no. 9109; Takara Bio, Inc.). The integrity of the RNA samples was assessed by agarose gel electrophoresis. RNA concentrations were quantified using a Qubit fluorometer (Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. Subsequent RNA-seq library preparation and deep sequencing were performed by BGI Genomics (Shenzhen, China). The samples were subjected to paired-end (2x150 bp) sequencing on an Illumina NovaSeq 6000 platform (cat. no. 20012850; Illumina, Inc.). Sequencing data were trimmed using Trimmomatic (v.0.39) (7) and aligned to the reference genome hg38 using STAR (v.2.7.7a) (8) software. Subsequently, read counting and differential gene expression analysis were performed using Cufflinks (v.2.2.1) (9). The data generated from the present study have been deposited in the NCBI Gene Expression Omnibus database under the accession number GSE316039.

Cell viability assays were performed as previously described (10). For cell viability assays, cells were seeded into 96-well plates (1x10³/well) and then incubated at 37˚C under 5% CO₂. The viability was measured at four different time points: 24, 48, 72 and 96 h. At each designated time point, 10 µl of Cell Counting Kit-8 (CCK-8) reagent (cat. no. HY-K0301; MedChemExpress) was added to each well, followed by incubation for an additional 2 h (the incubation duration was determined by a pre-experimental optimization to ensure the OD values were within the linear range). The optical density was then measured using a microplate spectrophotometer (BioTek Synergy H1; Agilent Technologies, Inc.).

The colony formation assay was performed as previously described (11). A total of 1x10³ shMGST3 and shNC transfected cells were seeded into 6-well plates. After 14 days, colonies were fixed with 4% paraformaldehyde (room temperature, 30 min), stained with 0.1% crystal violet (room temperature, 20 min) and imaged. Colonies containing ≥50 cells were quantified.

The scratch healing assay was performed to assess the migration capacity of cells in vitro (12). Cells were seeded in 6-well plates and cultured in medium containing 3% FBS. A scratch was created across the monolayer using a 200 µl pipette tip. Cell migration was observed at 0 and 12 h post-scratch using a fluorescence microscope (magnification, x4) to monitor wound closure. Images were captured to evaluate the migration levels of the transfected cell groups. The wound healing rate was calculated using the formula: Percentage change=([width at 0 h - width at 12 h]/width at 0 h) x100.

The Transwell assay was performed as previously described (13). The Transwell chamber contains a porous membrane with a pore size of 8 µm. To assess invasion capability, the bottom membrane of the chamber was coated with Matrigel at a ratio of 1:40 (Matrigel: serum-free medium; 100 µl). The lower chamber was filled with 500 µl of serum-free medium. The plates were then incubated at 37˚C for 2 h to allow the Matrigel to solidify. After 2 h, 5x10⁴ cells transfected with shMGST3 or shNC were seeded into the upper chamber, while the lower chamber was supplemented with 500 µl of DMEM containing 10% FBS. After a 24 h incubation at 37˚C in a 5% CO₂ incubator, cells on the membrane were fixed with 4% paraformaldehyde for 30 min at 25˚C, followed by staining with 0.1% crystal violet at 25˚C for 10 min. Following three PBS washes, non-migrated cells on the upper side of the membrane were removed with a cotton swab. In total, five random fields (magnification, x20) were selected under a fluorescence microscope using the brightfield mode and the number of stained cells was counted to quantify invasion.

Protein extraction and western blotting assays were performed as previously described (14). Cells were lysed on ice using RIPA buffer (cat. no. 9806; Cell Signaling Technology, Inc.) supplemented with PMSF and protease inhibitors. The supernatant was collected and boiled for 10 min after adding 5x SDS protein loading Buffer (cat. no. 30215ES; Shanghai Yeasen Biotechnology Co., Ltd.). A total of 10 µg of protein per lane was resolved on 4-20% SDS-PAGE gels (cat. no. ET12420Gel; ACE Biotechnology Co., Ltd.) and subsequently transferred to a PVDF membrane (cat. no. ISEQ00010; MilliporeSigma). The membrane was blocked with 5% skim milk at room temperature for 1 h, followed by overnight incubation at 4˚C with primary antibodies against MGST3 (cat. no. ab192254; Abcam) and GAPDH (cat. no. GB15002-100; Wuhan Servicebio Technology Co., Ltd.) diluted to 1:1,000. Following incubation with the primary antibody, the membrane was washed three times with TBST (containing 0.1% Triton X-100) and then incubated for 1 h at room temperature with horseradish peroxidase-conjugated goat anti-rabbit IgG (H+L) (cat. no. A0208; Beyotime Biotechnology) and goat anti-mouse IgG (H+L) (cat. no. A0216; Beyotime Biotechnology), each at a dilution of 1:5,000. After another three TBST washes, protein bands were visualized using an ultrasensitive ECL chemiluminescent substrate (cat. no. P0018FM; Beyotime Biotechnology) and imaged with a Tanon system equipped with Tanon Bio Imaging Analysis Software (version 1.0.0000; Tanon Science & Technology Co., Ltd.).

Clinical variables including tumor node metastasis staging and the corresponding gene expression matrices profiles of patients with colon adenocarcinoma (COAD) were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov).

Enrichment analysis were conducted on Metascape (https://metascape.org/) with the following screening criteria: A minimum overlap of 3 genes, P≤0.01 and a minimum enrichment factor of 1.5. The signals of the 14 pathways in the high and low expression samples of the glutathione gene set were calculated using the PROGENy package (15) in R software (version 4.2.3; https://www.R-project.org/), with samples stratified into high- and low-expression groups based on the median expression of the glutathione gene set derived from TCGA colon cancer data.

Search Tool for the Retrieval of Interacting Genes (https://cn.string-db.org) (16) was employed for constructing PPI networks using the gene set from Cluster 3. The resulting network file was imported into Cytoscape software (version 3.9.1) (17), which generated multiple subnetworks. The MCODE plugin (version 2.0.3) was applied to identify the most significant subnetwork based on the node score and number of genes, defined as a module with a node score >5 and containing >10 genes. Subsequently, the hub genes were selected from this subnetwork using the cytoHubba plugin (version 0.1).

TMB was defined as the total number of somatic coding base substitutions, insertions and deletions per megabase of examined genome (18). Using the TCGA-COAD cohort, stratified by the expression level of the glutathione gene set, TMB analysis and distribution of tumor mutation types were calculated by R package maftools (version 2.14.0) (19).

The EMT gene sets in both the study by Tan et al (20) and the EMTome database (http://www.emtome.org) were searched and the enrichment scores of the samples with high and low expression of the glutathione gene set were calculated using the single-sample Gene Set Enrichment Analysis (GSEA) algorithm (21).

The tumor immune infiltration analysis was performed using the CIBERSORT (version 0.1.0) (22) and ESTIMATE packages (version 1.0.13) (23) within the R (version 4.2.3) programming environment and based on gene expression matrices obtain from TCGA database.

Statistical analyses were performed as follows: Fisher's exact test was used for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEGs from the RNA-seq data; gene expression across different tumor stages in TCGA database were analyzed by Kruskal-Wallis test followed by Dunn's multiple comparisons test; a two-tailed unpaired t-test was applied to analyze the following metrics: TMB, immune score and the results of the Transwell and scratch healing assays, comparison of ferroptosis-related gene expression and pathway signals (14 pathways) between samples stratified by high and low expression levels of the glutathione gene set; two-way ANOVA test followed by Bonferroni's multiple comparisons test was employed for proliferation analysis of shMGST3 and shNC cells; and the Kolmogorov-Smirnov test was used for GSEA of the glutathione gene set. P<0.05 was considered to indicate a statistically significant difference.

Colorectal tumor samples from 5 male and 5 female patients, all aged 40-60 years were classified into early-stage tumors (pathological T1/T2), advanced tumors (pathological T3/T4), villous adenomas and tubular adenomas. This analytical approach included RNA-seq of patient samples, functional phenotypic validation in cell models and bioinformatics analyses of tumor-related indicators using data from TCGA database (TCGA-COAD cohort) (Fig. 1). Specifically, differential gene expression analysis [P<0.01, |log₂(fold change)| ≥1] was performed comparing adjacent normal tissues to the defined tumor stages (Fig. 2A). Through clustering of significant DEGs, four distinct expression patterns were identified: i) Cluster 1 (601 genes) showed specific upregulation in advanced tumors; ii) Cluster 2 (482 genes) exhibited high expression in normal tissues but progressive downregulation across precancerous lesions (villous/tubular adenomas) and tumors, with maximal suppression in advanced cases; iii) Cluster 3 (574 genes) demonstrated elevated expression in normal tissues that sharply declined in tubular adenomas and showed further reduction in villous adenomas and tumors; and iv) Cluster 4 (447 genes) displayed consistently lower expression in normal tissues vs. all pathological groups.

Among these clusters, Cluster 3 exhibited the most uniform and distinct expression patterns in CRC, differing markedly from both adjacent normal tissues and adenomas. Given this clear differentiation, Cluster 3 was selected for further investigation. KEGG pathway enrichment analysis revealed that the top enriched pathways in this cluster were predominantly metabolic-related (Fig. 2B). PPI network analysis using Cytoscape identified the ‘Glutathione metabolism’ pathway as containing the highest proportion of hub genes (Fig. 2C), suggesting its potential central role in CRC pathogenesis. Glutathione metabolism has been functionally implicated in CRC pathogenesis (24). Based on this network analysis, a glutathione metabolism gene set comprising eight key enzymes: MGST3, PRDX6, HPGDS, GGT1, GPX3, GSTM4, GSTA1 and GSR was identified. These genes were selected based on their hub status in the PPI network and their enrichment in glutathione-related pathways.

To validate these findings, transcriptomic data from TCGA-COAD were analyzed, stratifying the samples by clinical stage (early and advanced tumors). The glutathione gene set showed significantly higher expression in adjacent normal tissues compared with both tumor stages (Fig. 2D). This consistent downregulation across tumor progression suggests that glutathione pathway dysregulation may contribute to colon carcinogenesis. To further investigate the significance of the glutathione gene set, its expression patterns were analyzed across multiple malignancies. Consistently, reduced expression in tumor vs. normal adjacent tissue was observed in several tumor types, including breast invasive carcinoma, lung squamous cell carcinoma and stomach adenocarcinoma (Fig. S1). These findings were consistent with existing literature reports (25-29) and collectively demonstrated that glutathione pathway downregulation represents an early and conserved event in tumorigenesis across diverse cancer types, including CRC.

Given the established interplay between CRC progression and the tumor microenvironment (TME) (30), the potential role of the glutathione gene set was investigated in this context. Using TCGA-COAD data, the expression profiles of these genes were extracted and their mean expression levels in each sample were calculated. Samples were then stratified into high and low expression groups based on median gene set expression, enabling comparative analysis of some characteristics between these subgroups. Next, differences in TMB, EMT markers (31) and immune cell infiltration (32) between the high and low expression groups were evaluated. Notably, it was found that tumors with low expression of the glutathione gene set exhibited significantly elevated TMB compared with tumors with high expression (Fig. 3A). Conversely, tumors with high expression showed significantly enhanced EMT activity (Fig. 3B and C). These results indicated that elevated expression of glutathione pathway genes may be associated with tumor aggressiveness, as previously reported (33).

To characterize the tumor immune microenvironment, immune cell infiltration profiling was performed using CIBERSORT. Infiltrating immune cells were characterized into six functional groups, including lymphocytes, macrophages, dendritic cells, mast cells, eosinophils and neutrophils, using the method described by Wellenstein and de Visser (34). Comparative analysis revealed significantly elevated infiltration of lymphocytes, macrophages and neutrophils in tumors with high expression of the glutathione gene set (Fig. 3D). These findings were validated using the ESTIMATE package, which similarly demonstrated enhanced overall immune activity in tumors with high expression (Fig. 3E). Collectively, these results demonstrated that elevated glutathione pathway activity may be associated with a more immunologically active TME, suggesting these genes may serve as potential biomarkers for immunotherapy response (35).

Within the glutathione-related genes, decreased MGST3 expression may contribute to compromised carcinogen detoxification processes and further facilitate cellular transformation into a malignant state (36). To elucidate the functional consequences of MGST3 dysregulation in CRC, with a specific focus on colon carcinogenesis, stable MGST3 knockdown was established in the NCM460 normal human colon mucosal epithelial cell line (37). Western blot analysis confirmed efficient MGST3 depletion (Fig. 4A). Subsequent functional characterization revealed that shMGST3 transfection significantly enhanced cellular proliferation compared with shNC, as demonstrated by both CCK-8 viability and colony formation assays (Fig. 4B). These results indicated that MGST3 deficiency promoted proliferative advantages in colonic epithelial cells, potentially representing an early event in malignant transformation.

To further evaluate the malignant potential conferred by MGST3 depletion, the invasive and migratory capacities were assessed. MGST3-knockdown cells demonstrated a significantly increased invasion and migration rate in compared with the shNC cells (Fig. 4D and E). These functional assays collectively established that knocking down MGST3 expression enhanced both the invasive and migratory properties of colon epithelial cells.

Given the established connection between glutathione metabolism and ferroptosis regulation (38), key ferroptosis-associated genes in CRC (glutathione peroxidase 4, heat shock protein B1; solute carrier family 40 member 1 and lipocalin 2) were investigated (39-42). Comparative analysis revealed significantly reduced expression of these ferroptosis regulators in tumors with low expression of the glutathione gene set compared with tumors with high expression of the glutathione gene set (Fig. 5A). This consistent downregulation of key ferroptosis inhibitors suggests enhanced susceptibility to ferroptosis in tumors with lower expression of the glutathione gene set. Given that ferroptosis is closely related to the levels of reactive oxygen species (ROS) and fatty acid metabolism, GSEA was performed to systematically evaluate these pathways in the expression groups. Analysis of the glutathione gene set revealed a significantly stronger enrichment of both ROS-related and fatty acid metabolic gene signatures in the low expression group compared with the high expression group.(Fig. 5B and C). These results demonstrated that reduced glutathione pathway activity may be associated with molecular profiles characteristic of ferroptosis susceptibility, marked by enhanced ROS production and dysregulated lipid metabolism. Furthermore, analysis revealed a significant enrichment of the p53 pathway in samples where the glutathione gene set was highly expressed compared with samples with low expression (Fig. S2).

These integrated analyses demonstrated that elevated glutathione gene set expression confers resistance to ferroptosis in colon cancer cells, which potentially implies a more unfavorable prognosis for patients. This effect, potentially mediated through fatty acid metabolism and ROS regulation, suggests that glutathione gene set expression may serve as a potential therapeutic target.