Microbial and metabolic profiles unveil mutualistic microbe-microbe interaction in obesity-related colorectal cancer

Abstract

Obesity is a risk factor for colorectal cancer (CRC), and the involvement of gut microbiota in the pathogenesis of obesity and CRC is widely recognized. However, the landscape of fecal microbiome and metabolome distinguishing patients with obesity-related CRC from obesity remains unknown. Here, we utilize metagenomic sequencing and metabolomics from 522 patients with CRC and healthy controls to identify the characteristics of obese CRC. Our integrated analysis reveals that obesity-related CRC is characterized by elevated Peptostreptococcus stomatis, dysregulated fatty acids and phospholipids, and altered Kyoto Encyclopedia of Genes and Genomes pathways involving glycerophospholipid metabolism and lipopolysaccharide synthesis. Correlation analysis unveils microbial interactions in obesity, where the probiotic Faecalibacterium prausnitzii and the tumor-promoting species P. stomatis may engage in cross-feeding, thereby promoting tumorigenesis. In vitro experiments affirm enhanced growth under cross-feeding conditions. The mutualistic microbe-microbe interaction may contribute to the association between obesity and elevated CRC risk. Additionally, diagnostic models incorporating BMI-specific microbial biomarkers display promise for precise CRC screening.

Keywords: colorectal cancer; cross-feeding; gut microbiota; lipid metabolism; metabolome; metagenome; obesity; short-chain fatty acid.

Graphical abstract

Figure 1

Study design and global metagenomic and metabolomic characteristics of the fecal samples (A) A total of 539 fecal samples, including 228 healthy controls (CTRLs) and 311 patients with colorectal cancer (CRC), were collected and divided into four groups according to BMI: underweight group (gray), normal group (blue), overweight group (orange), and obesity group (red). The 17 samples in the underweight group were excluded for insufficient samples. 522 samples were randomly divided into the discovery cohort and the validation cohort. In the discovery cohort, we characterized the gut metagenome and metabolome, identified microbial and metabolic markers, and constructed a CRC classifier by the random forest models. The validation cohort was used to validate the diagnosis efficacy of the CRC classifier. Created with BioRender.com. (B) Alpha diversity measured by the Shannon index of patients with CRC and healthy controls in different groups. p values (two-sided Mann–Whitney U test) were calculated and shown in the figure. Data are delivered via the interquartile ranges (IQRs) with the median as a black horizontal line and the whiskers extending up to the most extreme points within 1.5× the IQR; outliers are represented as dots. (C and D) Principal-coordinate analysis (PCoA) of species (C) and metabolites (D). p values were calculated with adonis by 6,000 permutations. See also Tables S1 and S2.

Figure 2

Differential species in different BMI groups and the relationships between specific species and BMI in either patients with CRC or healthy controls (A) 296 differentially abundant bacteria species compared to the healthy controls in each of the three BMI groups are shown in the phylogenetic tree, grouped in the phyla Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria, and Fusobacteria. The innermost circle shows species log10 relative abundances averaged over all samples. In the outer rings, species are marked for significant (p adjusted by Benjamini-Hochberg < 0.05; multivariable regression model with age and gender as covariates) elevation (red) or depletion (blue). Bar plots show each group’s abundance log2 fold change of cancer-related species. Each plot shows bars in the four groups (normal, overweight, obesity, and the whole cohort) from left to right. Error bars are not applicable. (B) Co-occurrence relationships among species in normal weight, overweight, and obese groups performed by Spearman’s correlation. Only significant (adjusted p < 0.0005) correlations are shown. The red lines indicate positive species interactions, and the blue lines indicate negative interactions. See also Figure S1 and Table S3.

Figure 3

Metabolomics characteristics of patients with CRC and healthy controls in different BMI groups (A) Orthogonal partial least squares discriminant analysis (OPLS-DA) of metabolites in patients with CRC and healthy controls in each group. p values of the model are calculated by 1,000 permutations. (B) Volcano plot of differential metabolites in each group. p value adjusted by Benjamini-Hochberg < 0.1, absolute log2 fold change > 0.5, and variable importance for the projection (VIP) > 1 are considered differential metabolites. The x axis represents adjusted p value. The y axis represents log2 fold change. (C) Spearman’s correlation between differential species and metabolites in CRC and CTRL from normal weight, overweight, and obese groups. p values < 0.05 are marked with asterisks. ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001. See also Figure S2 and Tables S4 and S5.

Figure 4

CRC-associated functional alterations in KO genes (A) The heatmap shows the integrated meta-analysis that identified significantly changed KO gene expression in each KEGG pathway examined across three BMI groups. The cell color and intensity represent the log2 abundance fold change of KO genes. The cells identified and marked the significant differential KO gene (p adjusted by Benjamini-Hochberg < 0.05; multivariable regression model with age and gender as covariates). KO genes with absolute log2 fold change > 0.585 are shown in the figure. (B) Spearman’s correlation between differential species and KO genes. p value is adjusted by Benjamini-Hochberg. ∗Adjusted p < 0.05, ∗∗adjusted p < 0.01, and ∗∗∗adjusted p < 0.001. See also Figure S3 and Table S6.

Figure 5

The integrated network between species, KO genes, and metabolites and the cross-feeding diagram in the obese group (A) The integrated network between species, KO genes, and metabolites. Only significant (p < 0.05) links with absolute correlation coefficients above 0.2 are shown. The colors of nodes indicate the group of microbial features. The colors of the lines indicate the type of correlation (red, positive interactions; blue, negative interactions). (B) A schematic diagram illustrates the design of a cross-feeding experiment between F. prausnitzii and P. stomatis. (C) The bacteria concentration of F. prausnitzii and P. stomatis under cross-feeding conditions compared to monoculture controls. p values were calculated with Student's t test. ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001. See also Figures S4 and S5 and Table S7.

Figure 6

Microbial markers selected from each of the three BMI groups improved diagnostic efficacy (A) Heatmap of the importance of selected microbial features in each BMI group from the discovery cohort. Importance of each listed feature (belonging to the model constructed from each BMI group) for each group as estimated using the internal random forest “Gini importance” method. The microbial features not selected in one group are shown in gray. (B–D) Random forest models constructed from the selected features in (A). (B) The model constructed from N-CRC and N-CTRL of the discovery cohort was validated in the normal group, an overweight group, and an obese group of the validation cohort, and the AUROC values were calculated. (C) The model was constructed from Ov-CRC and Ov-CTRL of the discovery cohort and validated in each of the three groups from the validation cohort. (D) The model was constructed from the Ob-CRC and Ob-CTRL of the discovery cohort and validated in each of the three groups from the validation cohort. (E) The boxplots show the AUROC values of group-to-group validation for the models using selected features in (A). All boxplots represent the 25th–75th percentile of the distribution; the median is shown as a thick line in the middle of the box; the whiskers extend up to the most extreme points within a 1.5× IQR; and outliers are represented as dots. p values were calculated with Student's t test. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, and ∗∗∗∗p < 0.0001. See also Figure S6.