Patients of Chinese descent were enrolled from the Wuhan Union Hospital of Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China. The present study was approved by the Ethics Committee of Tongji Medical College of Huazhong University of Science and Technology (Approval no. 2014-041). All subjects provided written informed consent prior to their participation in the study. Samples and data of the participants were collected from hospital electronic medical records.

In the present study, 126 patients aged between 20 and 80 years (median age, 60.5 years) diagnosed histologically with primary colorectal cancer and possessing available clinical data, were enrolled. A total of 77 cases (61.11%) were female and 49 cases (38.89%) were male. The sample collection period spanned from June 2016 to December 2017. The patients with LM were assessed using liver magnetic resonance imaging and with no previous system treatment history, including chemoradiotherapy or targeted immunotherapy. Patients with a previous history of familial adenomatous polyposis, hereditary non-polyposis CRC, inflammatory bowel disease, metabolic diseases, infectious diseases, or immunodeficiency diseases, and also patients that had been administered antibiotics, corticosteroids or probiotics within 3 months prior to specimen collection were also excluded from the study.

A discovery cohort (cohort 2) and a validation cohort (cohort 3) were established. Colorectal cancer tissue samples were also concurrently tested and compared (cohort 1). The discovery cohort consisted of 18 LM and 36 NLM cases, and the validation cohort comprised 13 LM and 41 NLM cases. Samples from both cohorts were fresh feces.

To collect stool samples, ~40 g of feces were collected using a stool collection kit comprising a specimen receptacle, disposable sterile spoon, sterile conical tubes and disposable exam gloves. Participants promptly stored their stool samples at 4°C for a ≤3 h. Trained study staff aliquoted the stool samples and preserved them at −80°C for subsequent gut microbiome analyses. The fecal specimens were obtained during the initial hospitalization prior to anti-tumor treatment, including surgical intervention. Additionally, for the collection of tissue samples, a total of 18 primary carcinoma tissue samples (8 LM and 10 NLM; cohort 1) were collected. The tissue samples were fixed in formalin and embedded in paraffin (FFPE). Subsequently, five serial cuts (5 µm) per sample were placed in sterile microtubes and stored at room temperature until use in subsequent 16S rRNA MiSeq sequencing (Illumina, San Diego, U.S.A.). The samples were then frozen at −80°C for subsequent analyses.

DNA was extracted from the collected samples using the Omega Mag-Bind Soil DNA kit (Omega Bio-Tek, Inc.), following the manufacturer's instructions. The extracted DNA was quantified using a Nanodrop (Thermo Fisher Scientific, Inc.). Illumina sequencing was used to amplify the V3-V4 variable region of 16S rRNA of bacterial genome using reverse transcription-quantitative PCR (RT-qPCR). The sense and anti-sense primer sequences were 5′-ACTCCTACGGGAGGCAGCA-3′ and 5′-GGACTACHVGGGTWCTAAT-3′ respectively. The RT-qPCR amplification details have been indicated in a previous study by the authors (12).

The processing of raw sequencing data involved filtering, quality assessment, and the removal of query sequences. Initial filtering of Illumina MiSeq platform-generated data in FASTQ format utilized a sliding window method, employing a 10 bp window size and 1 bp step size. This process ensured an average sequencing accuracy of ≥99, a truncated sequence length of ≥150 bp, and the exclusion of ambiguous bases (N). Paired-end sequences from each library were overlapped using FLASH (version 1.2.7; ccb.jhu.edu/software/FLASH/), with criteria stipulating an overlapping base length of ≥10 bp and a mismatch base number less than 10% of the overlapping base length. Identification of mistaken sequences was conducted using QIIME software (version 1.8.0, http://qiime.org/), incorporating specifications such as a sequence length ≥160 bp, the absence of ambiguous bases (N), and constraints on 5′ primer mismatch base numbers and consecutive identical base numbers. Detection and removal of chimeric sequences were carried out through the combined use of USEARCH (version 5.2.236; drive5.com/usearch/) and QIIME software (version 1.9.1, http://qiime.org/).

Random forest algorithm was used to identify bacterial species that could distinguish the LM status. The top distinguishing phyla were included for further analysis. The ranking according to mean decrease in accuracy were obtained using default parameters in R (version 3.4.0, R core team, 2017) environment ‘Random Forest’. To further validate the predictive flora for LM, receiver operator characteristic (ROC) and area under curve (AUC) analysis were performed using MedCalc 16.1 (MedCalc Software bvba).

Quantitative Insights into Microbial Ecology 2.0 (QIIME2, version 1.9.1, http://qiime.org/) was used to analyze the obtained sequences and construct the taxa (13). Microbial alpha-diversity was analyzed using sampling-based operational taxonomic units OTUs and presented by Chao1, Faith's_pd and Observed_species indices. According to the total number of amplicon sequence variants (ASVs)/OTUs, R software was employed to generate rarefaction curve. Beta-diversity was assessed according to the microbial community structures using non-metric multidimensional scaling analysis 2 and principal coordinate analysis (PCoA) analysis. Linear discriminant analysis (LDA) effect size (LEfSe) was applied to analyze significant microbial characteristics. Additionally, the progression-free survival (PFS) and overall survival (OS) of the patients with CRC with high and low levels of Bacteroidetes and Proteobacteria in the primary cancerous tissues were compared.

Categorical data were analyzed using the Chi-squared test or Fisher's exact test, whereas the two independent samples t-test was used for continuous variables. Analysis was conducted using SPSS software (version 19.0, IBM SPSS). Kaplan-Meier survival analysis was employed, and the log-rank test was utilized to assess and compare the PFS and OS. P<0.05 was considered to indicate a statistically significant difference.

A flow chart of the experimental design is presented in Fig. 1. A total of 126 samples were analyzed in the present study. In particular, 54 fecal samples from the discovery cohort (cohort 2; LM=18 and NLM=36), 54 fecal samples from the validation cohort (cohort 3; LM=13 and NLM=41) and 18 primary tumor tissues from FFPE cohort as a supplementary discovery cohort (cohort 1; LM=8 and NLM=10) were included. The clinicopathological characteristics of the patients with LM and NLM from the three cohorts were analyzed (Table I). In the discovery cohort 1, there was no statistically significant differences between LM and NLM regarding sex, tumor location, and differentiation. Within discovery cohort 2, there were slight differences in sex and differentiation between LM and NLM. In validation cohort 3, there were slight differences in sex between LM and NLM. Our analysis of gut microbiota features influencing liver metastasis revealed that, given LM experienced liver metastasis while NLM did not, the clinical stages of most LM at the time of initial diagnosis consistently lagged behind (stage IV) in comparison with NLM (stage III) across all three cohorts.

DNA was extracted from the primary tumor tissues and fecal samples. Taxonomic profiling via 16S rRNA gene sequencing was performed. According to the correspondence between sequence similarity and bacterial taxonomic status, 97% similarity is generally considered for species division. The rarefaction curves displaying the sequencing depth reached a plateau, indicating a sufficient sequencing depth (Fig. 2A). The Venn diagrams demonstrated the number of shared and unique OTUs in each group. In total: i) 2,960 shared OTUs were observed among 42,613 OTUs in the LM and NLM groups of the discovery cohort; ii) 3,208 shared OUTs in a total number of 45,106 OTUs in the validation cohort and iii) 456 shared OUTs in a total number of 2,664 OTUs in the LM group and NLM group of the FFPE sample cohort (Fig. 2B). According to the Venn diagram, reduced numbers of total OTUs were detected in primary tumor tissues in comparison with the fecal samples.

Microbiota richness and diversity in the gut and primary tumors were subsequently compared between the LM and NLM groups. In the primary cancerous tissues of cohort 1, a discrepancy of alpha-diversity was observed between two groups, but only Chao1 index and Faith's_pd were not significant (P>0.05; Fig. 2C). In the fecal discovery (cohort 2) and the validation cohorts (cohort 3), the trend was consistent with four alpha-diversity index values (Chao1 index, Faith's_pd, Shannon and Observed_species between LM and NLM; Fig. 2C). Following beta-diversity analysis using PCoA, a clear separation of the community between the LM and NLM group in the discovery cohort (cohort 2) was observed (Fig. 3B). In the validation cohort (cohort 3), beta-diversity analysis also exhibited consistent results (Fig. 3C). These results suggested a higher microbiota richness and diversity in LM than in NLM group.

Subsequently, it was investigated whether there were differences in microbiota composition between the LM and NLM group. Therefore, the general landscape of microbiome was first assessed. For the top 10 microbiota at phylum level, nine species exhibited consistent preponderant enrichment in both the discovery and validation cohorts with fecal samples, including Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria, Verrucomicrobia, Synergistetes, Fusobacterium, Cyanobacteria, Saccharibacteria and Chloroflexi (Fig. 3E and F). Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria, Verrucomicrobia, Cyanobacteria and Chloroflexi were also largely enriched in microbial populations in primary tumor tissues (Fig. 3D). In the validation cohort, the main microbiota composition analysis at phylum level were consistent with those of the discovery cohort (Fig. 3D-F).

To further investigate these findings, LEfSe analysis was conducted, in order to detect markedly changed species with LDA score=4.0. The results indicated that the LM group with fecal samples was significantly enriched in Veillonellaceae at the phylum level, whereas the NLM group was significantly enriched with Porphyromonas at the phylum level (Fig. 4D and F). Additionally, the LM group of the tumor tissues were rich in Burkholderiales and Betaproteobacteria at the phylum level (Fig. 4B).

A species composition heatmap was then used to discover the changed phyla between LM group and NLM. The results indicated that Bacteroidetes, Planctomycetes, Fusobacterium, Firmicutes and Bacteroidetes presented consistent ascending enrichment in LM group compared with the NLM in the two specimen methods (Fig. 5C-E). Proteobacteria presented consistent descending enrichment in LM than NLM group in both tissues and fecal specimen cohorts. Additionally, OD1 presented reversed trend in LM group in two specimens (Fig. 5C-E).

Prognostic analysis was performed in cohort 2 and 3 between the LM and NLM groups and the results did not reveal statistically significant differences (Fig. 5). For the detection of microbiota for prognosis, the patients with NLM we stratified into the high vs. low categories based on the median relative abundance of the phylum. Fusobacterium^high was negatively associated with a shorter PFS in cohort 2 [hazard ratio (HR, 4.000; 95% confidence interval (CI), 1.1210-203.4000] and cohort 3 (HR, 2.933, 95% CI, 0.6008-15.6400; Fig. 6). The findings further indicated that the tumor microbiota could serve as a predictor of survival outcomes, suggesting the potential relevance of the microbiome in mediating CRC progression.

Subsequently, it was investigated whether microbiome communities could be used as a predictive biomarker for LM based on OTU abundance at the phylum level. Using the stochastic forest RF classification model, it was demonstrated that the Actinobacteria, Bacteroidetes, Firmicutes, Fusobacterium and Proteobacteria were highly discriminative phyla with the highest mean decreased score in cohort 2 and cohort 3 (Fig. 6). Thereafter, to define reliable biomarkers of the gut microbial response to LM, the top altered bacterial cluster was used to perform AUC-ROC analysis. In the discovery cohort of fecal samples, a relative high accuracy of Actinobacteria (AUC=0.603), Bacteroidetes (AUC=0.611), Firmicutes (AUC=0.650), Fusobacterium (AUC=0.898) and Proteobacteria (AUC=0.505) in predicting LM was demonstrated (Fig. 6C and E). Similarly, in the validation cohort of fecal samples, a relative high accuracy of Actinobacteria (AUC=0.612), Bacteroidetes (AUC=0.557), Firmicutes (AUC=0.623), Fusobacterium (AUC=0.893) and Proteobacteria (AUC=0.538) in predicting LM was demonstrated (Fig. 6D and F). Notably, Fusobacteria phylum in both fecal sample cohorts displayed a relatively high predictive value.