We selected patients who were surgically treated for colorectal cancer between January 2020 and December 2021. According to the UICC/AJCC TNM Staging System for CRC (8th edition, 2017), 92 patients with CRC were divided into a metastatic group (n = 48) and a non-metastatic group (n = 44) by specialist physicians. Exclusion criteria were anal canal tumors, appendiceal tumors, neuroendocrine tumors, familial adenomatous polyposis and cases of additional surgery for perforation or bleeding complicated by endoscopic treatment, colorectal cancer with antibiotics, glucocorticoids or immunosuppressive drugs used within 1 month at the time of sampling, and other colorectal cancer patients who could not be entered into this cohort.

Clinical samples were collected including tumor and tumor-adjacent tissues from CRC patients (the area surrounding the tumor <3 cm was considered adjacent tissue).

The collection of relevant clinical parameters of CRC patients included general clinical information, routine blood tests, biochemistry, coagulation function, tumor markers, histopathology, and other indicators.

Tissue-associated bacterial DNA was extracted from samples by using the QIAamp DNA Microbiome Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol. In brief, approximately 100 mg of intestinal tissue was homogenized, host cells were lysed, and host DNA was digested by benzonase (human DNase) while leaving the bacterial cells intact. Then, bacterial cells were concentrated by centrifugation, and bacterial DNA was extracted.

The sequencing procedure was performed as previously described (Emery et al., 2020). Briefly, the V3–V4 hypervariable region of bacterial 16S rDNA was amplified using universal sequencing primers 341F 5′-CCTAYGGGRBGCASCAG-3′ and 806R 5′-GGACTACNNGGGTATCTAAT-3′ (Yuan et al., 2018). The amplicon was sequenced by the Illumina MiSeq PE300 platform at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China). Raw FASTQ files were de-multiplexed and quality-filtered by QIIME1 (V1.9.1).¹ The optimized sequences were clustered into operational taxonomic units (OTUs) using UPARSE 7.1² with 97% sequence similarity level. The most abundant sequence for each OTU was selected as a representative sequence. To minimize the effects of sequencing depth on alpha- and beta-diversity measure, the number of 16S rRNA gene sequences from each sample were rarefied to 20,000. Bioinformatic analysis of the microbiota was based on the OTUs information. Alpha diversity indices including Chao1 richness and Shannon index were calculated with Mothur v1.30.1 (Emery et al., 2020). The similarity among the microbial communities in different samples was determined by non-metric multidimensional scaling analysis (NMDS) based on Bray-curtis dissimilarity using Vegan v2.5-3 package. The linear discriminant analysis (LDA) effect size (LEfSe) (Segata et al., 2011)³ was performed to identify the significantly abundant taxa (phylum to genera) of bacteria among the different groups (LDA score > 3, P < 0.05).

The main differential bacteria and Bacteroides fragilis toxin gene were confirmed by qPCR analysis. Briefly, experiments were performed with a QuantStudio 3 Real-time PCR System (Thermo Fisher Scientific, USA). The qPCR reaction system was: 2 × SYBR Green premix [Takara Bio technology (Beijing) Co., Ltd. Beijing, China] 5 μL, 1 μM forward and reverse primer sets (Table 1) 2 μL, 20 ng/μL DNA template 1 μL, ddH₂O 2 μL. The conditions of qPCR reaction were as follows: initial denaturation was done at 95°C for 60 s; amplification by using 45 cycles including denaturation at 95°C for 5 s, annealing and extension at 60°C for 30 s; melting curve was done at 95°C for 15 s, 60°C for 60 s; 95°C for 30 s.

SPSS 20.0 statistical software was used for statistical analyses. Patient characteristics were compared using unpaired Student’s t-test, Wilcoxon rank-sum test, or χ2 test as appropriate. Student’s t-test was used to analyse the differential bacteria between metastatic and non-metastatic CRC tissue. ROC curve analysis was used to determine the diagnostic value of serum biomarkers or selected bacteria in patients with CRC. Other diagnostic parameters were also evaluated, including sensitivity, specificity, cut-off value, and area under the ROC curve (AUC) with 95% confidence interval (CI), to assess the discrimination power of individual or combined biomarkers. A p-value less than 0.05 was considered statistically significant.

A total of 92 patients with CRC were included in the study, including 48 patients with metastatic CRC and 44 patients with non-metastatic CRC (Table 2). Statistical analysis of the basic clinical data indicated that age, gender, tumor size and location of tumor occurrence were not significantly associated with tumor metastasis. The differentiation degree was significantly related to CRC metastasis (p = 0.008), which is consistent with the understanding that if the tumor is less differentiated, it is more malignant and prone to metastasis (Derwinger et al., 2010). In addition, 38 of the 59 ulcerated CRC patients developed metastases, but only six of the 25 protuberant CRC patients were diagnosed with metastases (p = 0.001). This is because ulcerated CRC progresses deeper into the intestinal mucosa and is more likely to invade lymphatic and blood vessels, leading to CRC metastases (Bateman, 2022). Among the tumor markers alpha fetoprotein (AFP), carcinoembryonic antigen (CEA), and carbohydrate antigen 19-9 (CA19-9), the level of CEA was significantly associated with CRC metastasis.

To analyse the microbiota characteristics of tumor tissues in the metastatic and non-metastatic groups, we performed 16S rDNA amplicon high-throughput sequencing and subsequent bioinformatics analysis. The α-diversity of microbial communities is described by the Shannon and Chao indices (Ibrahim et al., 2019). The results showed that both indices of the two groups were not significantly different, indicating that the bacterial species diversity and richness were similar between metastatic and non-metastatic CRC tissues (Supplementary Figures 1A, B). Then, we applied non-metric multidimensional scaling analysis (NMDS) to analyse the β-diversity of the microbial communities. The results showed that the samples of non-metastatic group were not clustered together, and the β-diversity of two groups were not statistic significant (Supplementary Figure 1C).

We further analysed the composition of the microbial community from metastatic and non-metastatic CRC tissues. At the phylum level, the results showed that the intestinal bacteria in all tumor tissues were mainly from Firmicutes, Bacteroidota, Proteobacteria, Actinobacteriota, and Fusobacteriota (accounting for approximately 95%) (Figure 1A). The relative abundance of Bacteroidetes was significantly higher in the metastatic group than in the non-metastatic group (30.05 ± 21.20 vs. 18.35 ± 17.25%; P = 0.013), while the relative abundance of Proteobacteria was significantly lower in the metastatic group than in the non-metastatic group (9.87 ± 18.07 vs. 19.69 ± 29.13%; P = 0.009). In addition, Desulfobacterota, although the abundance was very low, was significantly increased in the metastatic group (0.82 ± 1.55 vs. 0.11 ± 0.35%; P < 0.001) (Figure 1B).

At the genus level, the bacterial composition in the tumor tissues was mainly composed of five genera: Bacteroides, Streptococcus, Escherichia-Shigella, Parvimonas, and Fusobacterium (Figure 1C). Among them, the relative abundance of Bacteroides was significantly higher in the metastatic group than in the non-metastatic group (21.67 ± 19.39 vs. 12.58 ± 12.93%; P = 0.049), while the relative abundances of Streptococcus and Escherichia-Shigella were significantly decreased in the metastatic group compared to the non-metastatic group (5.10 ± 11.9 vs. 23.12 ± 19.42%; P = 0.008 and 5.16 ± 14.65 vs. 11.66 ± 25.35%; P = 0.027, respectively) (Figure 1D).

The linear discriminant analysis (LDA) effect size (LEfSe; LDA score > 3.0) also found many differential bacterial species between the metastatic and non-metastatic CRC groups (Supplementary Figure 2). Interestingly, the metastatic CRC group had more relatively high abundance bacterial species than the non-metastatic CRC group, especially o_Bacteroidales, c_Bacteroidia, p_Bacteroidota, g_Bacteroides, f_Bacteroidaceae, g_Alistipes, f_Rikenellaceae, f_Oscillospiraceae, p_Desulfobacterota, and c_Desulfovibrionia, which had the highest scores, while c_Gammaproteobacteria, p_Proteobacteria, o_Enterobacterales, g_Escherichia-Shigella, f_Enterobacteriaceae, c_Bacilli, o_Lactobacillales, g_Streptococcus, f_Streptococcaceae, g_Curvibacter, o_Spirochaetales, and p_Spirochaetota were greatly enriched in the non-metastatic CRC group.

Bacterial composition analysis indicated that the abundance of g_Bacteroides was the highest among all the bacterial species and was increased in the metastatic CRC group. Then, we tried to identify the specific bacteria at the species level. In total, 37 OTUs were found to belong to g_Bacteroides, among which Bacteroides fragilis (OTU 2272), unclassified g_Bacteroides (OTU5969) and Bacteroides uniformis (OTU2249) were the three highest average abundance OTUs (Supplementary Table 1). Importantly, Bacteroides fragilis and Bacteroides uniformis showed a strong increasing trend in the metastatic CRC group compared to the non-metastatic group (p = 0.067 and 0.076, respectively) (Supplementary Table 1). As is reported that Bacteroides fragilis toxin (BFT) is the potential substance promoting tumorigenesis and metastasis (Zamani et al., 2019; Liu et al., 2020), we tested the bft gene frequency in non-metastatic and metastatic CRC tissue samples. The results showed that 28 out of 44 (63.6%) non-metastatic and 41 out of 48 (85.4%) metastatic CRC samples are bft gene positive (Supplementary Figure 3).

To confirm the high-throughput sequencing results, we performed a qPCR experiment to quantify the main differential species in tumor tissues. The results showed that at the phylum level, the abundance of Bacteroidota was increased in metastatic CRC tissue, while that of Proteobacteria was decreased (Figures 2A, B), but the results of Desulfobacterota were lacking because its abundance was lower than the limit of detection by qPCR in this study; at the genus level, the Bacteroides abundance increased, but the Streptococcus and Escherichia-Shigella abundances significantly decreased in the metastatic CRC group (Figures 2C–E). In addition, we tested the abundance of Bacteroides fragilis and Bacteroides uniformis, and the results showed that they were greatly increased in metastatic CRC tissues (Figures 2F, G).

First, we examined the diagnostic efficiency of the serum markers AFP, CEA, and CA19-9 in CRC metastasis. As shown in Figures 3A–C, the area under the curve (AUC) of CEA (0.652, 95% CI: 0.5387–0.7652, p = 0.012) was the largest, with a sensitivity of 0.479 and specificity of 0.8409 at the optimal cut-off value of 8.875. Next, we examined the diagnostic efficiency of differential bacteria in CRC metastasis. The AUCs of the ROC curves of p_Bacteroidota, p_Proteobacteria, p_Desulfobacterota, g_Bacteroides g_Streptococcus, and g_Escherichia-Shigella were 0.6709 (95% CI: 0.5609–0.7810, P = 0.0048), 0.5978 (95% CI: 0.4801–0.7155 P = 0.1065), 0.6906 (95% CI: 0.5822–0.7990, P = 0.0017), 0.6402 (95% CI: 0.5269–0.7534, P = 0.0207), 0.6449 (95% CI: 0.5313–0.7584, P = 0.0168), and 0.5630 (95% CI: 0.4430–0.6829, P = 0.2985), respectively (Figures 3D–I). These results indicated that the differential bacterial levels of CRC tissue groups possessed a moderate diagnostic efficiency for CRC metastasis.

Then, we attempted to improve the diagnostic efficacy by combining CEA with selected bacteria. The combination ROC curve of CEA with p_Bacteroidota, p_Proteobacteria, p_Desulfobacterota, g_Bacteroides and, or g_Streptococcus was drawn (Figures 4A–E), and the AUC was 0.6974 (95% CI: 0.5903–0.8046, P = 0.0011), 0.6723 (95% CI: 0.5628–0.7819, P = 0.0044), 0.7027 (95% CI: 0.5942–0.8111, P < 0.001), 0.6785 (95% CI: 0.5690–0.7880, P = 0.0032), 0.7055 (95% CI: 0.5996–0.8114, P < 0.001), respectively. Therefore, combination analyses obtained a higher diagnostic efficiency for CRC metastasis.

To explore whether the differential bacteria only existed in tumor tissue or existed in other normal intestinal tissues, we analysed the microbial composition of tumor-adjacent tissues from these patients. The results showed that at the phylum level, the relative abundances of Bacteroidota and Desulfobacterota were significantly higher in the adjacent tissue of the metastatic CRC group than in the adjacent tissue of the non-metastatic CRC group (30.78 ± 20.13 vs. 21.46 ± 20.42% and 0.56 ± 1.19 vs. 0.17 ± 0.4%, respectively) (Figures 5A, B), while the relative abundance of Proteobacteria showed a decreased trend in the adjacent tissue of the metastatic group compared to the non-metastatic group (6.02 ± 6.88 vs. 16.03 ± 25.72%) (Figures 5A, B), although there was no statistical significance. At the genus level, the relative abundance of Bacteroides in the metastatic group was significantly higher than that in the non-metastatic group (23.89 ± 19.31 vs. 18.16 ± 19.66%) (Figures 5C, D), while Streptococcus was significantly lower than that in the non-metastatic group (3.44 ± 5.6 vs. 5.9 ± 7.43%) (Figures 5C, D). These results indicated that the differential bacteria in the adjacent and tumor tissues of metastatic and non-metastatic CRC were consistent, meaning that the CRC metastasis associated bacteria are not specifically enriched in tumor tissues alone but are present in a larger area of the intestinal tract of CRC patients.