Animal experimentation was approved by the Animal Committee of the Chinese Academy of Sciences Institutional Laboratory [WIVA042020003]. In total, 70 female C57BL/6 mice (6-weeks-old, 20–24 g) were used in this study. AOM was purchased from Sigma-Aldrich (No. A5486, USA), and DSS was purchased from MP Biomedicals (No. 160110, CA). The animal experiments were conducted in two stages.

In the first stage, 35 mice were randomly divided into seven groups (n = 5) and treated with AOM/DSS at different concentrations as follows: Group A, 10 mg/kg AOM, 2% DSS; Group B, 10 mg/kg AOM, 1% DSS; Group C, 10 mg/kg AOM, 0.5% DSS; Group D, 10 mg/kg AOM, 0.25% DSS; Group E, 7.69 (10/1.3) mg/kg AOM, 2% DSS; Group F, 5.92 (10/1.3²) mg/kg AOM, 2% DSS; Group G, 4.55 (10/1.3³) mg/kg AOM, 2% DSS. The mice were intraperitoneally injected with AOM on the first day. Then, 1 week later, the mice were treated with DSS solution for 1 week, followed by 2 weeks of normal drinking water, for three cycles. All mice were euthanized until 14 weeks. The animal modeling process is shown in Figure 1A.

In the second stage, 35 mice were randomly divided into an experimental group (n = 30) and control group (n = 5). Here, we hypothesized that DSS at a certain dose might result in tumors in 50% of mice, without tumors in the other 50% of mice. The expected dose of DSS was calculated as 0.5359% according to the method in a previous study (Sanchez et al., 2018). Mice in the experimental group were treated with 10 mg/kg AOM and 0.5359% DSS in drinking water, based on the aforementioned procedure. Meanwhile, mice in the control group were maintained under standard conditions for 14 weeks. Finally, the colon tissues and fecal samples of all mice were collected for further investigation.

The fecal genomic DNA was extracted using a Stool DNA Kit (Qiagen, Germany) according to the manufacturer’s experimental steps. After determining the DNA integrity and concentration, the qualifying DNA samples were used for amplification. Specific primers were designed for the 16S rRNA V3–V4 region (F: 5′-CCTACGGGAGGCAGCAG-3′; R: 5′-GGACTACHVGGGTATCTAAT-3′). PCR amplification was performed using High-Fidelity PCR Master Mix with GC Buffer (New England Biolabs, USA). Then, the PCR products were separated using 2% agarose gel electrophoresis and magnetic beads and purified with a Gel Extraction Kit (Qiagen, Germany). The amplicon libraries were constructed using a TruSeq^® DNA PCR-Free Sample Preparation Kit (Illumina, USA), then quantified with a Qubit and qPCR, and finally sequenced on a NovaSeq6000 (Illumina, USA) platform. The effective tags were analyzed and obtained with the assistance of Novogene Biotechnology (Guangzhou, China). Operational taxonomic units (OTUs) were clustered based on tags of more than 97% identity using the Uparse v7.0.1001 method. Further analyses, including alpha and beta diversity, were subsequently performed.

Untargeted metabolomics were investigated via LC-MS/MS analyses. Fecal samples (100 mg) were placed in Eppendorf tubes and quickly treated with liquid nitrogen. Then, the samples were resuspended well with 80% methanol and 0.1% formic acid. After incubation for 5 min in an ice bath, the mixture was centrifuged at 15,000 × g and 4°C for 20 min. The supernatants were transferred and diluted with LC-MS grade water with methanol at a final concentration of 53%. Following another centrifugation step at 15,000 × g and 4°C for 15 min, the resulting supernatants were collected for subsequent experiments.

LC-MS/MS analyses were performed using the Vanquish UHPLC system and Orbitrap Q Exactive™HF-X mass spectrometer (Thermo Fisher Scientific, Germany) provided by Novogene (Beijing, China). The samples were injected into a Hypesil Gold column (2.1 mm × 100 mm, 1.9 μm), and the flow rate was 0.2 ml/min. Eluent A was 0.1% formic acid in water, and eluent B was methanol for the positive polarity mode, whereas for the negative polarity mode, eluent A was 5 mmol/L ammonium acetate in water, pH 9.0, and eluent B was methanol. The solvent gradient was set as follows: 1.5 min, 2% B; 12 min, 2–100% B; 14 min, 100% B; 14.1 min, 98% B; 17 min, 2% B. A QExactiveTMHF-X mass spectrometer was used with the source conditions as follows: sheath gas flow rate of 40 Arb, aux gas flow rate of 10 Arb, spray voltage of 3.2 kV, capillary temperature of 320°C. The raw data files were generated based on UHPLC-MS/MS and processed using Compound Discoverer 3.1 (Thermo Fisher, USA).

Qiime software (Version 1.9.1) was used to calculate Observed-OTU, Chao1, Shannon, and Simpson indices. Differences in alpha diversity indices among groups were analyzed based on the rarefaction curve and rank abundance curve. A Wilcox test was used for alpha diversity and beta diversity analyses. ANOSIM analysis was performed to test for differences in the microbial communities among groups. Raw data of LC-MS/MS were analyzed using Compound Discoverer 3.1. The differences in metabolic patterns among different groups were revealed based on partial least squares discrimination analysis (PLS-DA). To study phenotypic changes that might be caused by changes in the host microbial community structure, correlation analyses between the microbiome and metabolome were performed based on Pearson’s correlation analysis, correlation network diagram analysis, and correlation Sankey diagram analysis. P < 0.05 was considered statistically significant.

First, the mice were treated with AOM/DSS at different concentrations. After the initiation of tumorigenesis for 14 weeks, the mice were euthanized and investigated. We observed that mice in the group administered 1 and 2% DSS lost significantly more weight than those in the other two groups (Figure 1B). Moreover, the nodular tumors were macroscopically visible in the distal colon of mice in different groups (Figure 1C). In addition, the mice had obvious tumors in group A and B, whereas group D had no tumors (Table 1 and Figures 1C, D). However, the mice treated with 2% DSS and different concentrations of AOM had poor survival outcomes (Table 1). Therefore, we selected DSS as a variable factor to establish the target mouse models.

Here, we hypothesized that DSS at a specific dose would cause 50% of mice in a test population to develop tumors, and the theoretical dose was 0.5359% based on the results of groups A–D. Next, the mice in the experimental group (n = 30) were treated with 10 mg/kg AOM and 0.5359% DSS, and a negative control group of mice (n = 5) was also used in parallel. Fourteen weeks later, the mice were euthanized and investigated. As a result, 12 mice had tumors and 18 mice had no tumors in the experimental group. Then, we randomly selected the nine mice with tumors (tumor group), nine mice with no tumors (non-tumor group), and five control mice (control group) for further experiments and analysis.

To explore whether tumorigenesis is related to the gut microbiome, 16S rRNA sequencing was performed to identify gut microbiota profiles. In total, 1,189 OTUs were obtained among the three groups, comprising 828 in the control group, 797 in the tumor group, and 883 in the non-tumor group (Figure 2A). Moreover, a relative increase in bacterial richness was found in the non-tumor group, as revealed based on the rarefaction curve, compared with that in the other two groups (Figure 2B). In addition, the rank abundance curve yielded similar results (Figure 2C). To investigate bacterial diversity, we analyzed the alpha diversity indices and observed that there were statistically significant differences in the observed species, Shannon, Simpson, and Chao1 indices among different groups (Figures 2D–G). The principal component analysis showed that there were three separations of gut microbiota distributions among the three groups (Figure 2H). Non-metric multi-dimensional scaling analysis also revealed different distributions of microbial communities among the three groups (Figure 2I).

The gut microbial community structures at the phylum and genus levels in the three groups were analyzed, and the top 10 differences are presented in Figures 3A, B. The differences in the microbial distribution among the groups were determined through ANOSIM analysis (Supplementary Figures 1A–C). Next, the differential component proportions of microbes in each group were revealed based on the heatmaps (Figures 3C, D). In addition, at the phylum level, we observed that the Cyanobacteria, Proteobacteria, and Fusobacteriota were significantly enriched in the non-tumor group compared to abundances in the tumor groups, and Verrucomicrobiota was more abundant in controls (Supplementary Figures 1D–F). At the genus level, Prevotella, Alloprevotella, Neisseria, and Akkermansia exhibited marked differences among the groups (Supplementary Figures 1G–I).

To further determine the specific gut microflora associated with colorectal tumorigenesis, linear discriminant analysis effect size was performed among the three groups. The branching maps containing six levels, from phylum to species, reveled the signature microbiota. We found that the family Prevotellaceae and class Gammaproteobacteria may have a great effect in the non-tumor group, whereas the families Pseudomonadaceae and Akkermansiaceae, orders Pseudomonadales and Verrucomicrobiales, and class Verrucomicrobiae might play important roles in the control group (Figure 3E). Moreover, based on linear discriminant analysis scores (log₁₀) > 4, the histogram showed that four bacterial taxa, including Proteobacteria, Gammaproteobacteria, Prevotellaceae, and Prevotella, were enriched in the non-tumor group, two bacterial taxa were enriched in the tumor group, and 10 bacterial taxa were enriched in the control group (Figure 3F).

To identify the signature metabolites from fecal samples among the groups, we performed untargeted LC-MS/MS-based metabolomics. The PLS-DA showed significant differences between the fecal samples and metabolic phenotypes of different groups in both positive and negative ion modes (Figures 4A, B). In total, 1,112 and 554 metabolites were found to be changed in the tumor group, non-tumor group, and control group, in positive and negative ion modes, respectively (Supplementary Table 1). The differences in metabolites among the three groups are shown in Figure 4C. Further, we focused on the differences between the tumor group and non-tumor group. The heatmaps revealed the metabolite differences across each sample within the two groups (Figure 4D). The volcano plots also showed the significant upregulated or downregulated metabolites in the tumor group compared with levels in the non-tumor group (Figure 4E). Briefly, in positive ion mode, levels of 31 and 22 fecal metabolites were up- and downregulated, respectively, in the non-tumor group, with statistically significant differences compared to those in the tumor group. Meanwhile, in negative ion mode, levels of 10 and nine fecal metabolites were significantly up and downregulated, respectively, in the non-tumor group, with statistically significant differences compared to those in the tumor group. The structures of these metabolites were diverse, with many of the metabolites being either directly generated or modulated by the gut bacteria, including homogentisic acid, 3-methyladenine, and 2’-deoxyguanosine (downregulated in positive ion mode); nicotinic acid mononucleotide, N-acetyl-L-leucine, and linoleoyl ethanolamide (upregulated in positive ion mode); glycoursodeoxycholic acid, 2’-deoxyuridine, and pentadecanoic acid (downregulated in negative ion mode); and hydrocinnamic acid, oxoadipic acid, and 3-methyladipic acid (upregulated in negative ion mode). Moreover, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was used to identify the enriched pathways associated with differential metabolites in the two groups, and the top 20 most enriched pathways are listed in Figure 4F. Among them, styrene degradation and amino sugar-nucleotide sugar metabolism, were significantly altered in positive and negative ion modes, respectively.

To investigate the association between differential microbiota and metabolites in fecal samples, we conducted correlation analysis based on top 10 different bacteria at the genus level and top 20 different metabolites between the tumor group and non-tumor group. As shown in Figures 5A, B, the correlation heatmaps revealed the association between metabolites and microbiota, based on Pearson correlation coefficient analysis, in positive and negative ion methods. To further reveal the key bacterial flora and metabolites, we generated correlation network diagrams and observed that the connections were multiple and consanguineous (Figures 5C, D). Moreover, the correlation Sankey chart analysis also visually demonstrated the association between the gut microbiota and metabolites (Figures 5E, F). Notably, in positive ion mode, we observed that D-α-tocopherol was significantly negatively correlated with most microbiota, including Actinobacillus, Capnocytophaga, F0058, Lautropia, and Peptostreptococcus. Similarly, N1-(5-methylisoxazol-3-yl)-2-tetrahydro-1H-pyrrol-1-ylacetamide also exhibited negative correlations with most microbes. Meanwhile, 1,2-di(3,4-dimethoxyphenyl)diaz-1-ene, 3-methyl-5-oxo-5-(4-toluidino)pentanoic acid, oxymatrine, pantothenic acid, progesterone, and styrene showed positive correlations with the vast majority of microbe–metabolite pairs. In negative ion mode, the results showed highly negative associations for several microbe–metabolite pairs, such as Actinobacillus/L-methionine sulfone, Capnocytophaga/L-methionine sulfone, and Parasutterella/2’-deoxyuridine. However, Lautropia/glycoursodeoxycholic acid, Capnocytophaga/glycoursodeoxycholic acid, Lautropia/4-hydroxyisoleucine, and F0058/LPG 15:0, among others exhibited opposite relationships. Taken together, these results revealed significant correlations with respect to key microbe–metabolite pairs in the tumor and non-tumor groups, suggesting their potential roles in modulating tumorigenesis.