Altered bacteria and bacterial methyl donor-related functions in gut microbiome of CRC

Metagenomic sequencing was performed for fecal samples from 18 CRC patients and 18 healthy controls (HC) (Additional file 1: Fig. S1A). Alpha diversity indexes indicating community richness, diversity, and evenness were assessed via Richness, Shannon, Simpson, Pielou, Chao, and ACE indexes. No statistically significant changes in α diversity were observed in the comparison between CRC and HCs (Additional file 1: Fig. S2). Partial Least Squares Discriminant Analysis (PLS-DA) analyses were performed to evaluate β diversity. Overall fecal microbiome community of CRC patients was different from that of HCs as indicated by PLS-DA analysis (Fig. 1A). The comparison of relative species abundance between HCs and CRC patients revealed the enrichment of several previously reported species associated with CRC in CRC fecal samples, such as Fusobacterium nucleatum (LDA = 2.5, adjusted p = 0.044) and Bacteroides thetaiotaomicron patients (LDA = 4.2, adjusted p = 0.012) (Fig. 1B, Additional file 2: Table S4). Furthermore, the investigation of microbial functions associated with CRC identified 35 bacterial pathways that met our criteria (LDA > 2, α < 0.05) (Fig. 1C, Additional file 2: Table S5). Notably, multiple pathways related to methionine metabolism were prominently observed (Fig. 1C, highlighted in red). For instance, the COBALSYN − PWY (adenosylcobalamin salvage from cobinamide I, LDA = 3.8, adjusted p = 0.004) exhibits enrichment in the HCs. Adenosylcobalamin, also recognized as coenzyme B12, acts as a coenzyme, enabling the transfer of methyl groups during the synthesis of methionine. Nevertheless, the pathways associated with tetrahydrofolate biosynthesis, including both the super pathway of tetrahydrofolate biosynthesis (LDA = 3.1, adjusted p = 0.007) and the super pathway of tetrahydrofolate biosynthesis and salvage (LDA = 3.2, adjusted p = 0.008), exhibited notable over-representation in colorectal cancer (CRC) samples. Tetrahydrofolate, along with its derivatives collectively termed folates, serves as indispensable cofactors in one-carbon metabolism. These molecular entities play a pivotal role in transporting and donating C1-units, which are crucial for the synthesis of methionine and various other metabolites.

Fig. 1

The fecal microbiome in CRC patients and healthy controls. A β diversity (based on Bray–Curtis distances) evaluated by PLS-DA analysis. B Differential bacterial species between CRC and healthy controls (HC). C Differential bacterial pathways between CRC and HCs (pathways with LDA score > 3 are presented). Pathways related to methionine metabolism are highlighted in red. D The top 20 differential pathways identified by a meta-analysis based on the fecal microbiome of CRC and HCs derived from 10 public cohorts. E Forest plot of the S-adenosyl-L-methionine cycle I pathway

To further validate the findings of the dysregulation of methionine metabolism in CRC, we conducted a meta-analysis using fecal metagenomic datasets from 10 public CRC cohorts [30] (Additional file 1: Fig. S1B). The top 20 differential represented pathways between CRC and normal biopsies are presented in Fig. 1D. Consistent with our cohort analysis, we observed a significant alteration in methionine metabolism-related pathways. Particularly, the S-adenosyl-L-methionine cycle I pathway, where the S-adenosyl-L-homocysteine is recycled back to SAM, was reduced in CRC fecal samples across the cohorts (Fig. 1E). These findings further strengthen the evidence for dysregulated methyl donor-related pathways in CRC compared with HCs.

Altered methionine metabolism in tumor tissue microbiome

To explore the microenvironment regarding microbiome in tumor tissues directly, we conducted microbiome profiling using 16S rRNA gene amplicon sequencing in 24 pairs of tumor and adjacent normal (AN) biospecies (Additional file 1: Fig. S1C). There are 1008 ASVs, which are annotated to 171 bacterial genera and 146 species, passed our criteria (> 0.0001 in more than 10% of samples). The paired analysis revealed differences in microbial structure and composition between tumor and AN tissues of CRC patients (Fig. 2A–C). Though the alpha diversity remained comparable between tumor and AN tissues (Additional file 1: Fig. S3), significant dissimilarities were observed in β-diversity (Fig. 2A) as well as microbial abundance at multiple taxonomic levels (Fig. 2B, C). Bacteria that were reported to be associated with CRC, such as Parabacteroides (LDA = 3.7, adjusted p = 0.030), Alistipes (LDA = 3.0, adjusted p = 0.002), and Ruminococcus (LDA = 3.1, adjusted p = 0.039), were significantly enriched genera in tumor tissues (Additional file 2: Table S6). Furthermore, Clostridium hathewayi (LDA = 2.7, adjusted p = 0.049), a newly recognized CRC pathogen [31], showed a significant over-representation in tumor tissues.

Fig. 2

The overview of the microbiome in tumor and adjacent normal tissues. A β diversity (based on BrayCurtis distances) evaluated by PCoA analysis. B, C The differential taxa between CRC and AN tissues at genus (B) and species (C) level. D, E The differential microbial enzymatic functions (D) and MetaCyc pathways (pathways with LDA score > 2.5 are presented) (E) between tumor and AN tissues. F, G Forest plot of the associations between the relative abundance of L-methionine biosynthesis III pathway (F) and Methionine synthase function (G)

To further investigate the functional implications of the altered microbiota in CRC, we utilized PICRUSt2 to predict MetaCyc pathways and enzymatic functions of the microbiome. Specifically, we discerned notable significant enrichment of functions and pathways associated with methionine metabolism (Additional file 2: Table S7-8). For instance, “Adenosylmethionine decarboxylase” (LDA = 2.8, adjusted p = 0.010) and pathways related to L-methionine biosynthesis (L-methionine biosynthesis I: LDA = 2.3, adjusted p = 0.047; superpathway of L-methionine biosynthesis: LDA = 2.4, adjusted p = 0.040; superpathway of L-lysine, L-threonine, and L-methionine biosynthesis I: LDA = 2.3, adjusted p = 0.049) were prominently enriched in the AN tissues, as depicted in Fig. 2D, E. Adenosylmethionine decarboxylase catalyzes the conversion of SAM to S-adenosyl methioninamine, playing a pivotal role in the methionine salvage cycle. In agreement with the analysis of gut microbiota, these findings further highlight the potential involvement of methionine metabolism in maintaining normal DNA methylation patterns in colonic tissues.

Next, to strengthen our findings, we performed a meta-analysis by incorporating our data with four additional 16S rRNA sequencing datasets from CRC and matched AN tissues (Additional file 1: Fig. S1D). The analysis was based on the predicted MetaCyc pathways and Enzyme Commission (EC) classifications. We noted an enrichment of the L-methionine biosynthesis III pathway in the tumor samples, whereas the methionine synthase displayed enrichment in the adjacent normal samples, as illustrated in Fig. 2F, G. These findings underscore a substantial alteration in methionine metabolism between the tumor and adjacent normal tissues.

The association between fecal and tissue microbiome in CRC

A pivotal debate within microbiome research concerns the correlation between gut microbiota composition in feces and that in tissues [32, 33]. We investigated this relationship by examining the association between the fecal microbiome of colorectal cancer (CRC) patients and the tissue microbiome derived from both tumor and AN tissues (Additional file 1: Fig. S1E). At the genus level, the richness of the fecal microbiome did not exhibit any discernible association with that of the tissue microbiome for either tumor or AN tissues (Fig. 3A). However, the Shannon index, which measures both the richness and evenness of the microbiome was significantly higher in tissue samples (Fig. 3B). Beta diversity analysis revealed that the microbiome of fecal samples is distinct from that of tissue samples (Fig. 3C), and the distance from fecal microbiome to tumor microbiome is marginally significantly higher than that to AN microbiome (Fig. 3D). Therefore, these data seem to indicate that microbiota analysis in feces can be considered only a partial representation of the colorectal tissue microbiota.

Fig. 3

The relationship between microbiome derived from feces and tissues in CRC. A, B Alpha diversity comparison between feces and tumor/AN tissues from CRC patients. C, D Beta diversity calculated by PCoA of Bray–Curtis distance among microbiome of feces, tumor, and AN tissues. F, fecal; A, AN; T, tumor. E Positively correlated bacterial pathways between fecal microbiome and tumor or AN tissue microbiome. F Example pathways of the tumor-specific correlation

A taxa-wise correlation analysis between fecal and tissue microbiome revealed few significant associations (Additional file 2: Table S9). For example, the abundance of Eikenella is significantly correlated between fecal and tissue. Eikenella is a common inhabitant of the oral cavity and the intestinal and genital tracts, and a part of mucosal microbiota. It is considered to be an opportunistic pathogen leading to various infections. Next, we further assessed the congruence in microbial functions between the fecal and tissue microbiomes. Remarkably, we observed positive correlations in the abundance of 23 pathways between the fecal microbiome and tumor tissue microbiome. Likewise, 12 pathways displayed positive correlations between the fecal microbiome and AN tissue microbiome. These associations encompassed 13 pathways specific to tumor tissues, 12 shared pathways between both tumor and AN tissues, and 2 pathways exclusive to tumor tissues (Fig. 3E). This intriguing finding suggests that the functional attributes of fecal bacteria might partially reflect those of the tumor tissue. For instance, the L-glutamate degradation VIII pathway in fecal samples exhibited a robust correlation with the same pathway in tumor tissues, while this correlation was not observed in AN tissues (Fig. 3F). Notably, the meta-analysis conducted using publicly available fecal microbiome data indicated a significant increase in the L-glutamate degradation VIII pathway among tumor samples (Fig. 1D). This line of evidence suggests that specific functions within the fecal microbiome may indeed mirror those within CRC tumor tissues. This finding aligns with the perspective that insights gained from fecal samples could inform conclusions regarding the metabolic and functional profiles of the intestinal microbiota within the tumor microenvironment [34, 35].

Distinct DNA methylation profiles in CRC tumor and adjacent normal tissues

Genome-wide analysis of DNA methylation was performed using whole-genome bisulfite sequencing (WGBS) on 13 CRC tissue biopsies and their corresponding adjacent normal colonic tissue biopsies (AN) (Additional file 1: Fig. S1C). The CRC tissues exhibited a significantly lower overall methylation level compared to the matched AN tissues, as illustrated in Fig. 4A. Principal component analysis (PCA) analysis showed a sharp difference in the methylation profile between the tumor and AN tissues (Fig. 4B). These findings were further validated through additional analysis using 45 paired tumor and AN tissues from the TCGA-CRC cohort, as shown in Fig. 4C, D.

Fig. 4

The overview of DNA methylation profile in tumor and adjacent normal (AN) tissues. A The pair-wise comparison of the overall methylation level of the CRC tissues and its matched adjacent normal tissues. B PCA plot discriminated the tumor and AN tissues. C The average methylation level of TCGA CRC tumor and normal biospecies. D PCA plot discriminated the tumor and AN tissues in TCGA-CRC cohort. E Distribution of significantly hypo-methylated and hyper-methylated regions in a genomic region. F The heatmap of differentially methylated loci within promoter regions. G Functional enrichment of differentially methylated genes

To comprehensively investigate the genome-wide differential methylation patterns, we identified differentially methylated loci (DML) and regions (DMR). These regions were classified as either hypermethylated or hypomethylated based on an absolute average methylation difference (delta) of > 0.2 between CRC and AN tissues, representing a 20% change in methylation levels. A total of 210,171 hypomethylated regions and 1210 hypermethylated regions in CRC were identified. Notably, the promoter regions were highly hypomethylated, as demonstrated in Fig. 4E and F. Moreover, functional enrichment analysis revealed that the differentially methylated genes were enriched in functions of G-protein coupled receptor signaling pathway and ion transport (Fig. 4G), consistent with the analysis based on TCGA-CRC cohort data (Additional file 1: Fig. S4A).

Microbial-promoter methylation association was disturbed in the tumors tissues

To gain insight into the interplay between the microbiota and DNA methylation in colorectal cancer (CRC) tissues, we examined the association between microbial abundance and host methylation levels. Firstly, we performed a comprehensive analysis that involved correlating the abundance of bacterial taxa at both the genus and species levels with the overall methylation levels in both the host tumor and adjacent normal (AN) tissues. The results of this analysis, as visually represented in Fig. 5A, clearly demonstrate that the abundance of bacteria genera and species associated with global CpG methylation levels exhibited notable distinctions between tumor and AN tissues. This intriguing finding suggests that the microbial composition and its influence on DNA methylation may vary significantly between tumor tissues and the AN tissues in the context of colorectal cancer.

Fig. 5

Interactions between intra-tissue bacteria and host DNA methylation. A The correlation between bacterial abundance at the genus (left panel) and species (right panel) level with the total DNA methylation level in AN and tumor tissues. B The enrichment of bacteria-associated methylation loci in different genomic regions. C, D The associations between bacteria abundance and promoter CpG loci in the AN (C) and tumor (D) tissues of our clinical samples

Subsequently, we explored associations between individual host CpG loci and bacterial taxa. We employed a lasso penalized regression model to identify specific microbial taxa whose abundance was correlated with the methylation levels of CpG loci across the entire genome (Additional file 1: Fig. S1F). These models were fitted in a loci-wise manner, with the methylation level of each host CpG loci as the dependent variable and the abundances of microbial taxa as predictor variables. To identify robust associations, a stability selection methodology described by Priya et al. was utilized [29]. Consequently, we identified 2212 significant and stability-selected host CpG methylation-bacteria associations in tumor tissues and 1200 such associations in AN tissues of the TCGA-CRC cohort (Additional file 1: Fig. S1G, Fig. S5B-C).

Previous studies have indicated that specific microbial taxa can regulate the promoter methylation of individual host genes, thereby influencing gene expression [19]. Hence, we annotated the microbial-associated CpG loci to genomic regions. Using the total probes as background, we observed an enrichment of microbial-associated CpGs in promoter-associated regions, including the “TSS200,” “TSS1500,” “1st exon,” and “5′ UTR” region, but a depletion in the gene body and 3′ UTR region in the AN tissues. Strikingly, this enrichment was absent in the tumor tissues (Fig. 5B). These findings indicated a disturbance of microbial-promoter methylation regulation in the tumor microenvironment.

Microbes regulate methylation level of promoter regions

Next, we focused on the association between microbes and promoter CpG loci, uncovering 27 and 41 associations between microbial taxa and CpG loci in the AN and tumor tissues in our cohort, respectively (Fig. 5C, D). Notably, in the AN tissues, we observed multiple associations between the health-beneficial bacteria and gene methylation (Fig. 5C). For instance, Faecalibacterium prausnitzii, which is known for its abundance in the healthy human microbiota and depletion in various intestinal disorders, including CRC [36,37,38], was associated with the methylation of PPFIA2, NR3C1, and EXO1. Particularly, the methylation of the EXO1 promoter showed a positive association with the abundance of both F. prausnitzii and genus Faecalibacterium. EXO1 is involved in checkpoint progression and several DNA-repair pathways [39], and its deregulation is commonly observed in tumors [40, 41]. Analysis of TCGA data also revealed significantly increased EXO1 expression in colorectal tumors (Additional file 1: Fig. S5). Furthermore, Blautia obeum, a potential probiotic known to affect the composition of intestinal microbiota and inhibit the colonization and proliferation of pathogenic bacteria [42, 43], exhibited a strong negative association with multiple methylated loci (chr7:55963861, chr7:55963871and chr7:55963874) within NIPSNAP2 promoter. NIPSNAP2 is a mitochondrial membrane protein acting as a mitophagy receptor, the defective mitophagy has been increasingly associated with various diseases, including CRC [44]. These examples suggested a possible role of certain beneficial bacteria in maintaining normal cell functions through epigenetic regulation of gene expression in normal tissues.

In contrast, in the tumor tissues, we observed a greater number of associations involving potential opportunistic pathogenic bacteria (Fig. 5D). For instance, Bacteroides fragilis (B. fragilis) showed a positive association with C1D methylation, while Clostridium neonatale was positively associated with UNC5D promoter methylation. Impressively, Clostridium hathewayi, a recently identified pathogenic bacteria in CRC [19, 31], exhibited an association with the promoter methylation of SLC24A5 and ADHFE1.

Noncoding RNAs (ncRNAs) play a crucial role in gene expression regulation and are associated with various biological processes, including tumorigenesis in mammals. Dynamic DNA methylation patterns of ncRNAs have been reported during human embryonic development and disease progression [45,46,47]. WGBS identified a couple of methylation loci within the promoter of noncoding RNAs. In our analysis, we examined the association between noncoding RNA promoter methylation and tissue microbes, identifying 47 and 51 associations in AN and tumor tissues, respectively (Additional file 1: Fig. S7A-B). Notably, the methylation level of CCAT1, a well-known lncRNA involved in multiple tumors, including CRC [48], exhibited a negative association with Bacteroides ovatus in tumor tissues and a positive association with Veillonella parvulain AN tissues.

Furthermore, we observed associations between the methylation of different sets of small nucleolar RNAs (snoRNAs) and tissue microbiota in AN and tumor tissues (Additional file 1: Fig. S6A-B). snoRNAs are a group of regulatory RNAs mainly located in the nucleolus and guide the acquisition of 2′-O-methylation and pseudo-uridylation modify on ribosomal RNA (rRNA) and small nuclear RNA (snRNA). One of the two main families of snoRNAs, the C/D snoRNA genes, have been shown to be controlled by genomic imprinting [49]. Emerging evidence has demonstrated the significant roles of snoRNAs in cancer [50, 51]. For example, methylation of C/D snoRNA genes SNORD113-7 and SNORD114-1 showed associations with the abundance of Prevotella stercorea and Morganella morganii in the AN tissues. In tumor tissues, SNORD114-1 methylation was associated with Thermomonas fusca and Psychrobacter marincola abundance. Additional associations were observed between Parabacteroides gordonii and SNORD113-7, as well as between Bacteroides ovatus and SNORD114-23.

Collectively, our results highlight the potential involvement of microbial influences in the regulation of DNA methylation at gene promoter regions as well as noncoding RNAs. Moreover, the differential associations in the AN and tumor tissues implicate a distinct impact of microbial composition on the epigenetic landscape in these two contexts.

Tumor bacteria recurrently associated with host gene methylation

We proceeded to examine the recurrent association between microbes and promoter CpG methylation levels both in the TCGA-CRC cohort and our own data, revealing five genera that exhibited consistent associations with host promoter CpG methylation levels in tumor tissues. These genera include Fusobacterium, Stenotrophomonas, Bacteroides, Ralstonia, and Clostridium (Fig. 6A). Conversely, the abundance of Sutterrella, Dorea, Pseudomonas, Corynebacterium, and Prevotella was related to CpG methylation levels in the AN tissues across both datasets (Fig. 6B). Combining the gene methylation associated with these taxa in both TCGA-CRC and our datasets, we observed enrichment of MAPK cascade, ion transport, apoptotic process, and neuroactive ligand-receptor interaction-related functions in microbial-associated genes in tumor tissues (Fig. 6C). Conversely, in the AN tissues, genes were primarily enriched in carbohydrate binding functions (Fig. 6D). Most of the involved species were annotated to be contributed to metabolism of folate and methionine (Additional file 2: Table S10).

Fig. 6

Taxa recurrently associated with host gene methylation. A, B Taxa commonly associated with host promoter CpG methylation level in the TCGA-CRC cohort and our data (NJMU) in tumor (A) and AN (B) tissues. C, D The function enrichment of microbial-associated genes in AN (C) and tumor (D) tissues

Fusobacterium and Clostridium were well-known taxa associated with CRC progress. Xia et al. reported the role of Fusobacterium nucleatum and Clostridium hathewayi in driving promoter hypermethylation of tumor suppressor genes (TSG) in colorectal cancer [19]. In our analysis, Fusobacterium abundance was positively associated with the methylation level of C5AR1 promoter (cg10224107) in the TCGA-CRC cohort, and NT5E promoter (chr6:85449354) in our clinical samples. C5aR1 is a master regulator in colorectal tumorigenesis through immune modulation. C5a/C5aR1 signaling recruits myeloid-derived suppressor cells (MDSCs) into the inflamed colorectum, impairing CD8+T cells and modulating the production of critical cytokines and chemokines, thereby initiating CRC [52]. NT5E (CD73) functions as an inhibitory immune checkpoint molecule, suppressing cellular immune responses [53]. The positive association between Fusobacterium abundance and high methylation levels of the C5AR1 promoter is consistent with the potential role of Fusobacterium nucleatum in modulating immune responses in tumors, such as creating a proinflammatory microenvironment conducive to colorectal neoplasia progression [54] and enhancing the efficacy of checkpoint inhibitor blockade therapy [55, 56]. Similarly, Clostridium abundance was associated with hypermethylation of TSGs, including NFATC1 in TCGA-CRC cohort and UNC5D in our clinical samples. Stenotrophomonas, an opportunistic pathogen with increased colonization/infection in cancer patients [57], exhibited extensive positive associations with the hypermethylation of genes, including multiple TSGs in CRC, such as ESR1, RASSF5, DIRAS1, CADM1, ST5, GJB2, FAM172A, and HIVEP3. These findings further support the notion that bacteria in tumor tissue regulate gene expression through modulating methylation of promoter regions.

Microbial pathways related to host gene methylation

In pursuit of a deeper understanding, we investigated the intricate connection between microbial pathways and the promoter methylation status of host genes (Fig. 7A, B). Notably, a heightened complexity in correlations emerged within adjacent normal (AN) tissues, as compared to tumor tissues. For instance, the “superpathway of UDP-N-acetylglucosamine-derived O-antigen building blocks biosynthesis” exhibited associations with promoter methylation patterns across multiple genes, including C12orf45, GJD4, PNPLA8, MIR548X2, and LINC00380. Of significant interest, C12orf45 serves as a PAQosome cofactor, pivotal in facilitating the assembly of box C/D snoRNP [58]—an event observed to be elevated across various cancer types [59]. Equally noteworthy, the “coenzyme B biosynthesis” pathway demonstrated a positive correlation with NT5E promoter methylation. Importantly, several components of Coenzyme B, specifically B9 (Folate), B12, and B6, play a pivotal role in sustaining the one-carbon transfer cycles, which are integral for providing methyl donors essential for protein and DNA methylation [60].

Fig. 7

Interactions between microbes and microbial pathways and host DNA methylation. A, B Microbial pathways associated with host gene methylation level in AN (A) and tumor (B) tissues

As demonstrated above, an interrelation between bacterial functions within the fecal and tissue microbiome was established. Subsequently, we extended our analysis to explore associations between the abundance of fecal microbial functions and DNA methylation within tumor tissues. The outcomes unveiled three instances of negative associations between pathway abundance and gene promoter methylation levels (Additional file 2: Table S11). These instances encompassed the L-glutamate degeneration VIII pathway and OR8G5, the L-rhamnose degradation I pathway and C1D, and the superpathway of pyrimidine nucleobases salvage and NPFFR2. Notably, the abundance of the L-glutamate degeneration VIII pathway exhibited a distinctive positive correlation between fecal and tumor tissue microbiome, as depicted in Fig. 3B. This observation is particularly pertinent in light of downstream products resulting from glutamate metabolism, which have been posited to regulate chromatin modifications, including the intricate ten-eleven translocation (TET)-dependent DNA demethylation process [61,62,63]. Collectively, these findings underscore the multifaceted ways through which the microbiome potentially influences host DNA methylation.