This study employed a comprehensive multi-omics approach to investigate the causal pathways linking circulating metabolites with CRC risk through immune cell mediation. The analysis began with metabolome-wide association studies using large-scale GWAS data to identify 233 metabolite candidates, followed by rigorous MR analyses with sensitivity testing to establish causal relationships with CRC susceptibility. Colocalization analysis was further conducted to verify shared causal variants between metabolites and CRC risk. To elucidate potential immune mechanisms, we implemented a novel two-step analytical framework (26, 30): first identifying CRC-associated immune cell traits from 731 immunophenotypes, then examining whether causal metabolites influence these specific immune characteristics. This approach revealed critical immunological mediators in the metabolite-CRC axis. To further characterize epigenetic mechanisms, we retrieved epigenome-wide association study (EWAS) data for the identified metabolites and mapped their associated CpG sites. Methylation QTLs (mQTLs) corresponding to these CpG sites were used as instruments in MR analyses to identify metabolite-driven CpG methylation sites associated with CRC risk. To investigate potential transcriptomic links between metabolites and CRC, we obtained expression quantitative trait loci (eQTL) data for 19,960 genes from the eQTLGen consortium. We applied Summary-data-based MR (SMR) combined with the HEIDI test to prioritize candidate genes associated with CRC. These genes were intersected with eQTL genes driven by metabolite-related mQTLs, identified using the FUMA GWAS platform. Genes present in both datasets were considered as potential overlapping targets mediating CRC risk. Finally, to functionally validate the candidate target genes, we analyzed their expression patterns, prognostic relevance, and immune-related features using the TCGA CRC dataset. In vitro and in vivo experiments were conducted to confirm their roles in CRC-relevant phenotypes. This systematic strategy bridges metabolic dysregulation, immune modulation, and CRC pathogenesis through vertically integrated omics layers, providing mechanistic insights into cancer metabolic immunology ( Graphical Abstract ).

Metabolite instruments came from Karjalainen et al.’s NMR-based GWAS of 136,016 Europeans (171 lipoproteins, 23 lipids, 19 fatty acids, 20 other metabolites; GCST90301941 - 90302173) (29). Files were harmonized in R for MR. Immune-trait GWAS (731 phenotypes, n = 3,757) were drawn from a Sardinian cohort (IEU IDs ebi-a-GCST90001391–90002121) (31). CRC summary statistics (8,801 cases; 345,118 controls) were obtained from FinnGen Release 11 (“finngen_R11_C3_COLORECTAL_EXALLC”) with non-cancer controls. Omega-3 PUFA EWAS data (670 samples; trait EFO_0005110) and ARIES mQTLs provided methylation instruments (± 1 Mb). eQTLGen (n = 31,684) supplied cis-eQTLs. This dataset provides high-resolution cis-eQTL mappings derived from blood RNA, enabling reliable identification of genetic variants regulating gene expression (32). TCGA (COAD, READ, COADREAD) furnished tumor expression, immune-gene profiles, and survival data; samples with <30-day follow-up were excluded (33). All datasets comply with original ethical approvals; no additional consent was required for these secondary analyses.

We assessed causal relationships between circulating metabolites, immune cell traits, epigenetic CpG loci (exposures), and CRC (outcome) using genetic instruments as proxies. The TwoSampleMR package (v0.6.14) was used for data preprocessing and MR analysis. CpG-related mQTLs served directly as IVs. For metabolites, SNPs were filtered using a Bonferroni-corrected threshold (P < 1.8 × 10⁻⁹), derived by dividing P < 5 × 10⁻⁸ by 28 principal components explaining >95% of trait variance. This corrected for over-conservative adjustments in multi-phenotype GWAS (27, 29). Independence was ensured by retaining SNPs with r² < 0.001 within 10 Mb windows. Weak instruments (F-statistics < 10) were excluded, and traits with fewer than three instruments were removed to maintain power, due to high intercorrelation in NMR data. To reduce confounding, SNPs associated with CRC risk factors (e.g., alcohol, smoking, BMI, diabetes) were removed using LDlink (34).

Outcome associations were retrieved with extract_outcome_data(), and data harmonized using harmonise_data() to align effect alleles. Ambiguous SNPs (e.g., A/T, C/G) were handled cautiously, with palindromic variants near allele frequency 0.5 excluded (action = 2) to avoid strand ambiguity. MR estimates were then computed, with multiple sensitivity analyses conducted to assess robustness.

Causal relationship for circulating metabolites and CRC: To assess the causal impact of circulating metabolites on CRC, we applied MR methods. For exposures with multiple IVs, we used inverse variance weighted with multiplicative random-effects (IVW-MRE), assuming IV validity and independence. When only a single SNP was available, we used the Wald ratio (30). Given IVW’s sensitivity to directional pleiotropy, we also applied MR-Egger, weighted median, and radial MR to improve reliability. Additionally, RAPS (mr.raps v0.4.1) was used to adjust for weak instruments and pleiotropic bias (35). To address multiple testing across metabolites, Bonferroni correction was used; the significance threshold was P < 0.05 divided by the number of traits in each category. Associations passing this threshold in IVW or Wald ratio analyses were considered suggestive, and if confirmed in sensitivity tests, deemed robust. Significant metabolites were further tested for colocalization with CRC (coloc v5.2.3) (36). Strong evidence of shared genetic architecture was defined as PP.H4 > 0.90.

Two-step (Mediated) MR estimation (immune cells as mediators): To examine whether immune traits mediate metabolite effects on CRC, we used a two-step MR strategy (26, 30, 37):

Step 1: Identify immune traits causally linked to CRC (effect size = β1).

Step 2: Assess if these immune traits are influenced by CRC-associated metabolites (effect size = β2).

Assuming linear relationships, the indirect (mediated) effect of the metabolite on CRC through the immune trait was computed as: βmed=β1×β2. The direct (unmediated) effect was estimated as: βdirect=βtotal−βmed (Graphic abstract). This two-stage approach allowed us to identify metabolite-immune-CRC pathways.

Epigenetic MR: CpG methylation effects on CRC were analyzed via the Wald ratio for single-SNP sites or IVW-MRE for multiple IVs, supported by weighted median and MR-PRESSO. A Bonferroni threshold of P < 0.002 (0.05/23) determined significance.

SMR and FUMAGWAS integration: We integrated eQTL and CRC GWAS data using SMR (v1.3.1) to identify candidate genes. SNPs with MAF < 0.01 were excluded, and only cis-eQTLs (P < 5 × 10⁻⁸) within ±2 Mb were considered. Significant probes passed SMR testing and Bonferroni-adjusted HEIDI tests (P_HEIDI > 0.05). In parallel, FUMAGWAS SNP2GENE was used to map Omega-3 PUFA mQTL-linked genes. Mapping employed 1000 Genomes EUR data with positional, eQTL (P < 0.05), and 3D chromatin interaction strategies. The final gene list included those identified by both SMR and FUMA.

MR-based causal inference can be biased by directional pleiotropy, particularly in omics research. To evaluate robustness, we assessed instrument heterogeneity via mr_heterogeneity(), which tests for variability among instruments that could violate MR assumptions. Directional pleiotropy was tested using the MR-Egger intercept (mr_pleiotropy_test()), with P > 0.05 indicating low pleiotropic bias. To further detect and correct pleiotropic effects, we applied Radial MR (RadialMR v1.1), using both IVW and MR-Egger estimates. Radial plots helped identify and exclude outlier SNPs violating MR assumptions (38). The leave-one-out analysis was conducted to assess the significance of the driving effect of individual IVs.

Tumor vs. normal gene expression was analyzed across cancers using unpaired Wilcoxon rank-sum and signed-rank tests. Prognostic relevance was assessed via Cox proportional hazards models (coxph from survival v3.2-7), with log-rank tests for significance. The “pROC” package was also utilized to create receiver operating characteristic (ROC) curves and to calculate area under the curve (AUC). For combined diagnostic performance, a multivariable logistic regression model was utilized to create an integrated ROC curve. With the “decision_curve” function from the “rmda” package, decision curve analysis (DCA) was conducted, including for single-gene models and for multigene models integrating two genes. The analysis was under the assumption of a case–control study design, and probability thresholds from 0 to 1 were assessed. For immune profiling, gene expression was mapped to gene symbols, and immune cell infiltration scores (e.g., B cells, CD4⁺, CD8⁺ T cells) were computed using TIMER and deconvo_CIBERSORT (IOBR package). We also analyzed expression of SLC6A19 and 60 immune checkpoint markers (24 inhibitory, 36 stimulatory) (39), along with 150 immune pathway genes (41 chemokines, 18 receptors, 21 MHCs, 24 inhibitors, 46 stimulators). Pearson correlations between SLC6A19 and these immune genes were computed. All workflows (data harmonization, MR, and outlier filtering) were performed in R v4.3.1 with a reproducible pipeline.

The single-cell landscape of SLC6A19 in CRC was profiled with the GEO dataset GSE166555. All preprocessing and downstream analyses were carried out in Seurat v5.1.0. After SCTransform normalization, highly variable genes were identified and subjected to linear dimensionality reduction; the first 20 principal components were then embedded in two dimensions with UMAP. Cluster-specific differentially expressed genes were obtained with FindAllMarkers and manually cross-checked against the CellMarker 2.0 reference sets. Gene-level expression values were extracted with FetchData function. Human intestinal epithelial cells NCM460 and CRC cell lines HCT116, SW480, and CACO2 were obtained from Procell (China). NCM460, HCT116, and SW480 were cultured in DMEM with 10% FBS, while CACO2 cells were maintained in DMEM/F-12 with 15% FBS. All cells were incubated at 37°C in 5% CO₂ and authenticated by STR profiling within the last three years. For SLC6A19 knockdown, three siRNAs and a non-targeting control (si-NC) were used (GenePharma, China). Cells (5 × 10⁵/well) were seeded in 6-well plates and transfected at 70 – 80% confluence with 50 nM siRNA and 5 μL siRNA-mate plus in 250 μL Opti-MEM. After a 15-min incubation, the mix was added to cells; media was refreshed at 6 h, and incubation continued for 48 h. Knockdown efficiency was verified by RT-qPCR. For overexpression, 2 μg of pcDNA3.1-SLC6A19 plasmid was transfected with 5 μL GP-transfect-Mate in 250 μL Opti-MEM. After 6 h, media was changed, followed by 48 h of incubation. Controls included si-NC and empty vectors. Western blotting was performed on denatured lysates using SDS-PAGE and PVDF membranes (0.45 μm, Millipore). Membranes were blocked in 5% milk/TBST, probed with anti-SLC6A19 (Proteintech, Cat# 27575 - 1-AP, 1:1000) and HRP-conjugated secondary antibody (Cat# SA00001 - 2, 1:3000). Protein bands were visualized with ECL reagents.

For the Cell Counting Kit-8 (CCK - 8) assay, transfected CRC cells were seeded into 96-well plates at a density of 2,000 cells/well, with six replicates per group. At 0, 12, 24, and 48 hours, 10 μL of CCK - 8 reagent (Dojindo, Japan) was added to each well, and absorbance at 450 nm was measured after 2 hours of incubation to assess cell viability. Wound healing assays were performed by scratching monolayers of overexpressing cells at 90% confluence in 6-well plates using a 200 μL pipette tip. After PBS washing, serum-free medium was added, and wound closure was imaged at 0 and 48 hours under an inverted microscope. For Transwell migration and invasion assays, cells were serum-starved for 12 hours, trypsinized, and resuspended in serum-free medium at 4 × 10⁴ cells/mL. For migration, 300 μL of cell suspension was added to the upper chamber (8 μm pore size; Corning) with 600 μL of medium containing 20% FBS in the lower chamber. For invasion, 100 μL of Matrigel (Corning) was pre-coated on the upper chamber membrane and solidified at 37°C for 1 – 2 hours. After 24 hours, cells were fixed with 4% paraformaldehyde, stained with 0.1% crystal violet, imaged, and counted. All in-vivo protocols were approved by the Institutional Animal Care and Use Committee of Shantou University Medical College. Six-week-old male BALB/c nude mice were kept under specific-pathogen-free conditions. Two HCT116 derivatives, one bearing an empty control vector (HCT116-Vector) and the other engineered to over-express SLC6A19 (HCT116-OE-SLC6A119), were injected subcutaneously at 5×10⁸ cells in 150 µL (n = 4 per group). Animals were sacrificed 1 – 4 weeks after implantation. Tumors were excised, fixed in 10% neutral-buffered formalin for SLC6A19 immunohistochemistry, and portions were snap-frozen at -80°C for further analyses. The entire experiment was repeated independently on three occasions. Statistical analyses were performed using GraphPad Prism (version 10.4). Two-group comparisons were analyzed using independent sample t-tests, while comparisons among multiple groups were conducted using one-way ANOVA. A P-value < 0.05 was considered statistically significant.

We retrieved genetic IVs for 233 circulating metabolites. Two traits—conjugated linoleic acid to total fatty acids (CLAbyFA) and glycerol—were excluded due to having only two valid instruments, leaving 231 metabolites for MR analysis. These included 171 lipoproteins, 23 lipids, 18 fatty acids, and 19 amino acid/other metabolites, each represented by 3 – 81 IVs ( Supplementary Table 1 ). IV strength was supported by F-statistics: 36.18 – 2152.62 (lipoproteins), 36.19 – 1915.87 (lipids), 36.19 – 5088.20 (fatty acids), and 36.33 – 7610.08 (amino acid/others), all well above the threshold of 10.

Causal effects on CRC risk were estimated using MR IVW. No lipoproteins met the Bonferroni threshold (P < 2.92 × 10⁻⁴; 0.05/171) ( Figure 1A , Supplementary Table 2 ), nor did any of the 23 lipid traits (P < 2.17 × 10⁻³) ( Figure 1B , Supplementary Table 3 ). However, three metabolites—two fatty acid-related and one amino acid/other—showed significant associations. Higher DHAbyFA [OR = 1.26, 95% CI: 1.10 – 1.43, P = 5.58 × 10⁻⁴] and FAw3byFA [OR = 1.22, 95% CI: 1.13 – 1.31, P = 2.51 × 10⁻⁷] levels were associated with increased CRC risk (Bonferroni cutoff: P < 2.78 × 10⁻³) ( Figure 1C , Supplementary Table 4 ). Conversely, lower pyruvate (Pyr) levels acted as a protective factor against CRC [OR (95% CI) = 0.59 (0.49 – 0.70), P = 6.81 × 10⁻⁹], passing the corrected threshold for the amino acid/other group (P < 2.63 × 10⁻³) ( Figure 1D , Supplementary Table 5 ). Robustness checks (weighted median, Egger, IVW radial, RAPS) confirmed the direction and significance of effects ( Figures 2A, B ). RAPS estimates aligned with IVW: DHAbyFA (OR = 1.28, P = 5.35 × 10⁻⁵), FAw3byFA (OR = 1.22, P = 3.32 × 10⁻⁵), and Pyr (OR = 0.58, P = 4.52 × 10⁻³).

No significant pleiotropy or heterogeneity was detected for these metabolites (P > 0.05; Table 1 ), supporting valid instruments. One outlier SNP (rs4601123) was identified in DHAbyFA-CRC analysis ( Figure 2C ). After removing it, causal associations remained consistent across methods, with no residual pleiotropy ( Figure 2D ). Colocalization analysis showed strong evidence for shared loci with CRC for FAw3byFA (PP.H4 ≈ 0.97), but not for DHAbyFA (PP.H3 ≈ 0.81) or Pyr (PP.H1 ≈ 0.72), suggesting the latter two act via independent loci ( Figures 2E, F ).

Together, these findings suggest FAw3byFA as a likely causal biomarker for CRC, while DHAbyFA and Pyr may be biologically relevant but genetically distinct contributors.

Given the immune system’s central role in CRC pathogenesis, we next tested whether genetically predicted immune cell traits mediate the metabolite–CRC link. IVs were identified for 731 immune phenotypes, of which 614 met the inclusion criteria, with F-statistics ranging from 29.85 to 5062.70—indicating strong instrument validity. Thirty-nine immune traits were nominally associated with CRC (P < 0.05); 16 showed positive associations and 23 showed inverse associations ( Figure 3A , Supplementary Table 6 ). RAPS-based MR confirmed 30 associations: 18 protective and 12 risk-related phenotypes ( Figure 3B , Table 2 , Supplementary Table 7 ).

We then conducted two-step MR to explore mediation by immune cells in the FAw3byFA–CRC relationship. In step 1, 30 immune phenotypes with confirmed CRC links were retained. In step 2, IVW MR found FAw3byFA to be positively associated with Effector Memory CD4⁺ T cell Absolute Count (β = 0.16, P = 0.03), showing nominal significance ( Figure 3C , Supplementary Table 8 ). Sensitivity methods, including MR-Egger (β = 0.33, P = 0.01), supported this result ( Figures 3C, D , Supplementary Table 9 ). No evidence of heterogeneity or pleiotropy was observed ( Supplementary Tables 10 , 11 ), reinforcing the robustness of the findings.

These results suggest that increased FAw3byFA levels are genetically linked to elevated Effector Memory CD4⁺ T cell counts, which are independently associated with higher CRC risk. To quantify mediation, we decomposed the total effect: FAw3byFA’s total effect on CRC was β_total = 0.20, with indirect effect β_med = β₁ × β₂ = 0.02 (β₁ = 0.12 for immune cell to CRC; β₂ = 0.16 for FAw3byFA to immune cell), and a direct effect of β_direct = 0.18. This indicates a mediation proportion of 10% (0.02/0.20).

Overall, these data identify Effector Memory CD4⁺ T cell Absolute Count as a partial genetic mediator of FAw3byFA’s causal effect on CRC, accounting for a quantifiable share of its impact.

Recent studies have highlighted that DNA methylation, as an epigenetic modification, significantly influences metabolite levels and may consequently disrupt gene expression and contribute to disease pathogenesis. In this study, we investigated the EWAS findings for omega-3 PUFAs and performed epigenetic MR analysis. A total of 47 CpG sites associated with omega-3 PUFA levels were identified ( Supplementary Table 12 ). From the ARIES mQTL database, we further identified 996 SNPs representing mQTLs associated with 29 of these CpG sites, among which 23 were available for outcome analysis ( Supplementary Table 13 ). Using IVW and Wald ratio methods, we found that nine CpG sites exhibited causal associations with CRC, surpassing the multiple testing correction threshold (P < 0.002; 0.05/23) ( Figure 4A ). Additional analyses using RAPS, weighted median, and MR-PRESSO methods consistently highlighted strong causal links for cg05181941, cg06817802, and cg22456785. Specifically, cg22456785 (OR = 0.89, 95% CI: [0.89 – 0.89], P = 1.00E - 302) and cg06817802 (OR = 0.95, 95% CI: [0.94 – 0.96], P = 7.04E - 24) were found to be protective, whereas cg05181941 (OR = 1.02, 95% CI: [1.02 – 1.03], P = 2.05E - 14) was identified as a risk factor ( Figure 4B , Supplementary Tables 14 , 15 ). Sensitivity analyses revealed no significant heterogeneity in the MR estimates. In pleiotropy assessment, cg22456785 showed no evidence of horizontal pleiotropy based on the MR-Egger intercept test (P = 0.987), whereas cg05181941 and cg06817802 displayed potential pleiotropic effects ( Table 3 ). To determine whether any single mQTL was disproportionately influencing the causal effect estimate, we performed a leave-one-out analysis. Results indicated that the observed associations for cg05181941 and cg06817802 were not driven by any individual SNP ( Figures 4C, D ).

Overall, elevated methylation at cg05181941 was associated with increased CRC risk, while lower methylation at cg06817802 and cg22456785 was linked to reduced CRC susceptibility ( Figure 5A ). These findings underscore the potential of omega-3-related DNA methylation as biomarkers for CRC susceptibility.

To identify genes underlying omega-3 CpG methylation, we mapped CpG-associated mQTLs to eQTLs using transcriptomic data from 19,960 genes. SMR analysis (threshold P < 0.0025; 0.05/20) identified 98 genes ( Figure 5B ), 31 of which were excluded due to pleiotropy (HEIDI P < 0.05), leaving 67 CRC candidate genes ( Figure 5C , Supplementary Tables 16 , 17 ).

Using FUMA’s SNP2GENE, we assessed whether mQTLs for the three significant CpGs also served as eQTLs. Among 428 mQTLs, 224 unique eQTL target genes were found across whole blood and 118 tissues ( Figure 5D , Supplementary Table 18 ). Overlap with the 67 SMR-identified genes revealed seven shared candidates: C1QB, LPCAT1, NDUFS6, SLC6A19, SRSF10, SRSF6, and YWHAH ( Supplementary Table 19 ). Four—LPCAT1, NDUFS6, SLC6A19, and SRSF10—were expressed in T cells and colon tissue, suggesting possible immune-mediated mechanisms ( Figure 5D ; Table 4 ).

Differential expression analysis showed consistent upregulation of LPCAT1, NDUFS6, and SRSF10, while SLC6A19 was markedly downregulated in COAD, READ, and COADREAD subtypes. Notably, SLC6A19, which encodes an amino acid transporter, showed the most significant reduction across CRC types and was expressed in both colon and T cells ( Figures 5D, E ). Immunoblotting in CRC cell lines (HCT116, SW480, CACO2) confirmed reduced SLC6A19 protein levels compared to normal epithelial cells (NCM460) ( Figure 5F ), suggesting a potential functional role in CRC development.

SLC6A19 was identified as a causal CRC gene via SMR and mapped to the omega-3-linked CpG site cg06817802, sharing common mQTL/eQTL SNPs ( Figure 6A , Supplementary Table 20 ). Cox survival models showed that low SLC6A19 expression correlated with poorer outcomes in READ (HR = 0.81, P = 0.03) ( Figure 6B ). ROC analysis showed that both SLC6A19 and the established CRC biomarker CEA had moderate diagnostic accuracy for distinguishing CRC from normal tissues, with CEA performing slightly better (AUC = 0.829 vs. 0.789) in the READ cohort ( Figure 6C ). DCA indicated that the combined SLC6A19-CEA model provided greater net clinical benefit than either marker alone across most threshold probabilities ( Figure 6D ). Consistently, log-rank survival analysis confirmed that patients with low SLC6A19 expression had markedly worse survival outcomes (HR = 0.24, 95% CI: 0.08 – 0.752, P = 7.3 × 10⁻³) ( Figure 6E ). Building on prior mediation findings implicating CD4⁺ T cells in the omega-3–CRC pathway, we analyzed immune cell infiltration. Using TIMER and CIBERSORT, SLC6A19 expression was positively associated with activated CD4⁺ memory T cells in COAD and COADREAD ( Figure 6F ) and correlated with CD4⁺ T cell infiltration in all subtypes (r = 0.12 – 0.22) ( Figure 6G ). In addition, we investigated the association between SLC6A19 expression and immune regulation. Pearson correlation analyses were performed between SLC6A19 and 150 immune pathway marker genes (representing chemokine, receptor, MHC, immunoinhibitor, and immunostimulator classes), as well as 60 immune checkpoint genes (24 inhibitory and 36 stimulatory). The results indicated varying degrees of expression correlation across CRC subtypes, suggesting that SLC6A19 is broadly involved in immune modulation ( Figure 7 ). Together, these data support a link between SLC6A19 expression and CD4⁺ T cell activity, reinforcing its role in the omega-3–CRC axis.

To explore what SLC6A19 might be doing in colorectal cancer, we first looked at its expression at single-cell resolution. We found that the gene shows up mainly in normal epithelial cells, and also in specific subsets of T cells within the tumor microenvironment ( Figures 8A, B ). To validate its function, we overexpressed and silenced SLC6A19 in HCT116, SW480, and CACO2 cells. CCK - 8 assays showed that overexpression inhibited proliferation, while knockdown enhanced it ( Figure 8C ). Wound healing and Transwell assays revealed reduced migration and invasion upon overexpression. Specifically, wound closure was significantly delayed in all lines ( Figure 8D ). Similarly, in the Transwell invasion assays, SLC6A19 overexpression significantly impaired the invasive capacity of CRC cells. The number of cells traversing the Matrigel-coated membrane was notably reduced. Quantitative data revealed significantly lower invasion rates in HCT116 (P < 0.001), SW480 (P < 0.01), and CACO2 (P < 0.001) cells ( Figure 8E ). We then tested this in mice. HCT116 cells with stable SLC6A19 overexpression, or empty vector controls, were implanted into BALB/c nude mice. Tumors in the SLC6A19 group showed stronger SLC6A19 staining in immunohistochemistry ( Figure 8F ) and reduced tumor size ( Figure 8G ). Tumors in the overexpression group grew more slowly ( Figure 8H ), were lighter at the endpoint ( Figure 8I ). These findings confirm a tumor-suppressive role for SLC6A19 in CRC.