A total of 43 patients newly diagnosed with CRC at the 108 Military Central Hospital and 44 healthy controls from the community, who were in generally good health and showed no symptoms of CRC, were recruited for this study during the years 2022–2023. In addition, the following inclusion criteria were used for the selection of CRC and healthy subjects: aged 50–79 years, male or female, and body mass index (BMI; range: 18.5–30). If any of the following conditions applied, a participant was considered ineligible to participate in this study: complicating diseases, mental illness, evidence of a biochemical or structural abnormality of the digestive tract, irritable bowel syndrome, inflammatory bowel disease, taking any antibiotics within 4 weeks prior to recruitment, or adherence to a restricted diet.

Participants' demographic characteristics such as age, sex, BMI, health status, and medication history were collected and controlled for the data analysis. A family history of cancer and gastrointestinal and allergic disorders was also recorded. The results of blood chemistry tests and clinicopathological characteristics of the CRC patients, including tumor size, TNM (tumor, node, and metastasis) stage, tumor location, metastasis status, tumor recurrence, and survival status, were extracted from patients' medical records.

All participants were asked to self-collect fecal samples at home or at the hospital with a plastic specimen container. Fresh feces were transported to the laboratory on dry ice to prevent bacterial growth within 4 h or were kept in cold storage (~4°C) for 24–48 h before arriving at the laboratory. Upon arrival at the laboratory, the fecal samples were subsampled into aliquots to prevent further unnecessary freeze-thaw cycles and were frozen for long-term storage at −80°C.

Total genomic DNA was extracted from the fecal samples using the Qiagen DNeasy^® PowerSoil^® Pro Kit (QIAGEN, Germany) with a bead-beating step. Briefly, 0.25 g of fecal samples were homogenized mechanically in two intervals (10 min at a speed of 25 Hz) in a PowerBead Pro Tube containing CD1 solution and beads using a TissueLyser II (QIAGEN, Germany). The subsequent steps of genomic DNA extraction were then performed according to the manufacturer's extraction instructions. Genomic DNA was measured using the NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, USA). The genomic DNA was kept at −20°C for further use.

The extracted genomic DNA was used for PCR amplification with a bacterial primer pair (5′- CCTACGGRRBGCASCAGKVRVGAAT-3′ and 5′- GGACTACNVGGGTWTCTAATCC-3′) for the 16S rRNA gene V3–V4 hypervariable region. The PCR master mix contained 1 μM of each forward and reverse primer, 2 μl of dNTPs, 2.5 μl of TransStart buffer, 0.5 μL of TransStart Taq DNA polymerase, 20 ng of the extracted DNA, and 25 μl of nuclease-free water. The PCR program was as follows: 1. 94°C for 3 min, 2. 95°C for 5 s, 3. 57°C for 90 s, 4. 72°C for 10 s, and 5. 72°C for 5 min, with steps 2–4 repeated for 24 cycles. The amplicons from each sample were then incorporated with indexed adapters using limited cycle PCR and were purified using magnetic beads. The purified 16S amplicon libraries were quantified using the Tecan Infinite 200 Pro microplate reader and were pooled together at equimolar concentration. The 250 bp paired-end sequencing was processed on the pooled libraries using the Illumina MiSeq Platform (Illumina, USA) at Genewiz, Inc. (South Plainfield, NJ, USA). Raw 16S rRNA gene sequencing data is available at the NCBI Sequence Read Archive (SRA), under BioProject PRJNA1077687.

Analysis of 16S amplicon sequencing data was carried out as in our previous study (Tran et al., 2019). Briefly, after merging the 2x250 bp paired-end reads with FLASH v1.2.8 (Magoc and Salzberg, 2011), the merged reads were trimmed forward primers and reverse primers with cutadapt v1.8.3 (Martin, 2011) and truncate_reverse_primer.py script from Quantitative Insights Into Microbial Ecology (QIIME) v1.9.1 (Caporaso et al., 2010). The QIIME's split_libraries_fastq.py script was used to filter out low-quality reads with the quality threshold at Q19 and concatenate reads from different samples into one file. Next, the clustering of reads into operational taxonomic units (OTUs) was processed using USEARCH v8.1 software (Edgar, 2010). The reads were pre-processed using usearch -derep_fulllength, sortbylength, and sortbysize. A default threshold of 97% sequence identity was set for usearch -cluster_otus to form OTU representative sequences. Chimeras were eliminated from the OTU representative sequences using the usearch -uchime_ref with the ChimeraSlayer reference database. The usearch -usearch_global was applied to map all reads to the OTU representative sequences with an identity threshold of 97% to create the OTU table. Taxonomy assignments were obtained using mothur v1.36.1 (Schloss et al., 2009) with the Ribosomal Database Project (RDP) trainset 18 (Cole et al., 2009) and SPINGO v1.3 (Allard et al., 2015) for species classification. Alpha-diversity and beta-diversity were computed using QIIME's alpha_rarefaction.py and beta_diversity.py scripts on a rarefied OTU table at the minimum OTU reads across all samples as 70,034 reads per sample.

Statistical analysis was performed using R v4.2.1 software (R Core Team, 2016). The significant differences in clinical measures and alpha-diversity indices were determined using the Mann–Whitney U-test for comparison of two independent groups and the Kruskal–Wallis test for comparison of more than two groups. Comparisons of categorical variables were calculated using Fisher's exact test for two groups and chi-squared test for >2 groups. Principal coordinate analysis (PCoA) was used to visualize the similarities in beta diversity between the samples. Permutational multivariate analysis of variance (PERMANOVA) tests were conducted to assess the significant differences in beta diversity between the groups. The differential taxonomic abundances between CRC patients and healthy individuals were calculated using DESeq2 (Love et al., 2014), controlling for age, sex, BMI, and hypertension. A p-value adjusted using the Benjamini–Hochberg correction, < 0.05, was considered statistically significant. The Spearman correlation between the relative abundance of OTUs and the clinical traits was conducted using the cor.test function in R.

The clinical descriptors of the CRC patients and healthy controls enrolled in this study are presented in Table 1. We included 43 CRC patients and 44 healthy controls, with an average age of 64 years (age range, 53–76 years) for CRC patients and an average of 61 years (age range, 51–75 years) for healthy controls. The ages were divided into three groups; the majority of healthy controls (56.8%) were 60–70 years old, while the patients with CRC were 50–60 and 60–70 years old, accounting for 65.7%. There was no significant difference between CRC patients and healthy controls with respect to BMI, sex, or family history of CRC (Mann–Whitney U-test, p-value > 0.05). A high prevalence of T2D and hypertension, which are common comorbidities in people aged over 50, was recorded in CRC patients compared to healthy individuals.

All patients with CRC suffered from abdominal pain. The coexistence of other symptoms included rectal bleeding, diarrhea, and/or weight loss. Among the 43 patients, five patients (11.6%) had either liver metastasis or peritoneal metastasis, while two patients (4.7%) experienced multiple organ metastases. A total of 62.5% of tumors in CRC patients were reported as 2–5 cm in size and 39.5% were 5–10 cm; all tumors were resectable. The majority of the tumors (46.5%) occurred in the left colon, followed by the right colon, the rectum, and the transverse colon, accounting for 23.3, 20.9, and 9.3%, respectively. According to the TNM classification, the recruited patients ranged from early to advanced stages, with the majority (72.1%) being in stages II and III. Cancer recurred in 11.6% of the patients, and three patients (7%) died within 2 years of diagnosis (Table 2). In addition to clinical features, blood biochemical characteristics were investigated in CRC patients. The mean values of most parameters in both male and female CRC patients were within the normal range. This indicates that CRC itself was not the cause of the blood test changes. On the contrary, carcinoembryonic antigen (CEA) level, a prognostic marker for CRC, was elevated in the serum of 21 patients (48.8%), especially patients with metastases (Supplementary Table 1).

Total genomic DNA was extracted from fecal samples of individual participants and underwent 16S rRNA gene amplicon sequencing for the identification of bacterial diversity and composition. We observed a significant increase in alpha diversity measured using the Chao1, Shannon's, Simpson's, and phylogenetic diversity indices in Vietnamese patients with CRC compared to healthy controls (Mann–Whitney U-test, p-value < 0.05; Figure 1A). The alpha-diversity indices did not differ by the participant's age category (50–60, 60–70 vs. 70–80 years), BMI category (normal weight vs. overweight), sex (male vs. female), and blood pressure (hypertension vs. non-hypertension), although the Simpson's index was slightly lower in the participants with T2D (median = 0.93; IQR = 0.90–0.94) than the participants without T2D (median = 0.94; IQR = 0.93–0.96; p-value = 0.098). Furthermore, alpha diversity in the CRC group was not influenced by any clinical traits, except for the BMI category and T2D status. In contrast, the Chao1 index differed among age categories in the healthy control group (Supplementary Table 2).

Upon assessment of the overall gut microbiota composition, beta diversity represented by unweighted UniFrac distance revealed two fairly distinct clusters based on CRC condition (PERMANOVA: R² = 0.035; p = 0.001), whereas no significance was detected between CRC patients and healthy controls based on the weighted UniFrac distance (PERMANOVA: R² = 0.017; p = 0.155; Figure 1B). In addition, there were no significant differences by the participant's age category, BMI category, sex, T2D condition, or blood pressure in either UniFrac distance. However, unweighted and weighted UniFrac distances were discriminatory in terms of BMI categories in the CRC cohort, with a low coefficient of determination (PERMANOVA: R² = 0.037 and 0.058, respectively; p < 0.05 for both; Supplementary Figure 1A and Table 3). Within the healthy subjects' samples, there was a significant separation among age categories based on the weighted UniFrac distance (PERMANOVA: R² = 0.098; p = 0.01), and sex had a relevant impact on the weighted UniFrac distance (PERMANOVA: R² = 0.043; p = 0.015; Supplementary Figure 1B and Table 3).

We observed high variability in microbiome composition at the family and genus levels among individuals in both the CRC and control groups. The most dominant families in the participants' gut were Lachnospiraceae, Oscillospiraceae, Bacteroidaceae, Enterobacteriaceae, and Peptostreptococcaceae, with a mean relative abundance of 74.6% (range = 44.01–82.18%, SD = 11.24; Supplementary Figure 2A). At the genus level, the relative abundance means of the top five genera were 16.81% (SD = 10.63) for Blautia, 5.69% (SD = 4.83) for Faecalibacterium, 4.68% (SD = 6.22) for Phocaeicola, 4.21% (SD = 7.46) for Escherichia/Shigella, and 4% (SD = 5.05) for Bacteroides (Supplementary Figure 2B). Next, we identified 6 classified families and 18 classified genera that significantly correlated with the CRC condition using DESeq2, controlling for age, sex, BMI, and hypertension (Figure 2A). Compared to the healthy individuals, Porphyromonadaceae, Turicibacteraceae, Eubacteriales incertae sedis XIII, Akkermansiaceae, and Peptoniphilaceae were found to be increased in CRC samples, while Micrococcaceae was negatively associated with CRC condition. Similarly, the genera Agathobacter, Parasutterella, Citrobacter, Muribaculum, Butyricicoccus, and Lactobacillus were underrepresented in the CRC patients. On the contrary, 12 out of 18 genera were positively associated with CRC, with the top genera showing the highest log2 fold changes in abundance in CRC compared to healthy controls, namely Mogibacterium, Porphyromonas, Peptostreptococcus, Prevotella, and Parvimonas (log2 fold change > 2.5, adjusted p-value < 0.001; Figure 2A and Supplementary Table 3).

A comparison of bacterial species detected a lower abundance of three species and a higher abundance of 13 species in the CRC patients than in the healthy controls (Figure 2 and Supplementary Table 4). Among those, Parvimonas micra, Peptostreptococcus stomatis, and Prevotella intermedia showed the greatest variance between the CRC and healthy groups (log2 fold change > 4, adjusted p-value < 0.001). In addition, we found significant enrichment of three Bacteroides (B.) species (B. caccae, B. cellulosilyticus, and B. salyersiae) and three Parabacteroides (P.) species (P. distasonis, P. merdae, and P. goldsteinii) in patients with CRC. On the contrary, CRC patients had lower abundance levels of two health-associated species including Lactobacillus johnsonii and Butyricicoccus pullicaecorum than the healthy individuals (log2 fold change < −2, adjusted p-value < 0.05).

T2D is a common comorbidity in patients with CRC, and this metabolic disease has been reported to be associated with alteration in the gut microbiota community (Umirah et al., 2021). Thus, we further evaluated CRC-related gut microbiota signatures in the presence of T2D. A comparison between CRC patients without T2D (n = 35) and those with T2D (n = 8) revealed significant differences in 17 genera and 15 species (Figure 3A). Among the significantly associated species, the top six species that showed the most significant fold change between the two groups were Ruminococcus flavefaciens, Peptococcus niger, Akkermansia muciniphila, Butyricicoccus pullicaecorum, Phascolarctobacterium faecium, and Porphyromonas endodontalis (log2 fold change > 3.9, adjusted p-value < 0.05; Figure 3A and Supplementary Figure 3).

Next, in comparison to the healthy control group, the families Turicibacteraceae and Peptoniphilaceae and six genera including Turicibacter, Mogibacterium, Parvimonas, Prevotella, Peptostreptococcus, and Porphyromonas were consistently enriched in CRC patients with and without T2D (log2 fold change > 1, adjusted p-value < 0.05; Supplementary Figure 4). At the species level, only Parvimonas micra, Peptostreptococcus stomatis, and Prevotella intermedia were associated with CRC patients with either non-T2D or T2D. On the other hand, the potential health-associated species were not shared between CRC/non-T2D and CRC/T2D. For example, the two species Lactobacillus johnsonii and Bifidobacterium longum were significantly increased in healthy individuals when compared to CRC/non-T2D. Meanwhile, the combination of CRC and T2D caused the change in the addition of species, increasing Turicibacter sanguinis, Alistipes finegoldii, Bacteroides fragilis, and Ruminococcus obeum and reducing Phascolarctobacterium faecium, Ruminococcus bromii, Ruminococcus flavefaciens, Parasutterella excrementihominis, Blautia glucerasea, Alistipes putredinis, Bacteroides uniformis, and Bilophila wadsworthia (Figure 3B and Supplementary Table 4). However, a larger sample size for the CRC/T2D group would be needed to confirm this finding. Together, these results suggested that both CRC and T2D contributed significantly to the alterations of the gut microbiota composition.

We examined whether the composition of the human gut of CRC patients differed based on clinical traits such as CEA level and tumor size. Since the T2D condition was found to be correlated with changes in the gut composition in CRC patients, we adjusted for the presence of T2D along with age, sex, BMI, and hypertension when conducting the differential abundance analysis of the gut microbiota for the level of CEA and the tumor size using DESeq2. Overall, 26 OTUs showed significant variation between the CEA-normal and CEA-elevated groups, while 25 OTUs were significantly associated with the tumor size (Figure 4). The CEA level and the tumor size shared three common microbial OTUs: OTU_1340, belonging to the species Bacteroides fluxus; OTU_101 and OTU_386, belonging to the genus Streptococcus and the genus Flintibacter, respectively.

Next, the Spearman correlation analysis in the CRC subgroups (CRC/non-T2D and CRC/T2D) revealed the differences in OTUs correlated with the CEA level and the tumor size (Supplementary Figure 5 and Supplementary Table 5). We observed strong correlations between microbial OTUs and the CEA level and the tumor size in the CRC/T2D group (Spearman's correlation coefficient, 0.7 < |ρ| < 0.9; nominal p-values < 0.05); meanwhile, moderate relationships were shown in the CRC/non-T2D group (Spearman's correlation coefficient, 0.3 < |ρ| < 0.6; nominal p-values < 0.05). Specifically, the relative abundances of several OTUs belonging to the Lachnospiraceae were inversely correlated with tumor size in both CRC/non-T2D (OTU_1275, OTU_1278, and OTU_1541) and CRC/T2D (OTU_83, OTU_176, OTU_385, OTU_1104, OTU_1198, and OTU_1541). Furthermore, the OTUs that correlated with either the CEA level or the tumor size were separated into three distinct clusters by the hierarchical cluster analysis. Cluster 1 was composed of five out of a total of 9 OTUs that were positively associated with the tumor size in CRC/T2D samples, whereas the majority of negative associations with CEA level and tumor size were observed in cluster 2 and cluster 3 in both the CRC subgroups, except for the association between the CEA level and its correlated OTUs in the CRC/T2D group. However, the microbial abundance heatmap of the correlated OTUs showed no clear separation of CRC samples according to TNM stage, tumor location, and metastasis (Supplementary Figure 5).