Characteristics of gut microbiota in anastomotic leakage patients in cold zones post-colorectal cancer surgery: A high-throughput sequencing and propensity-score matching study

2.1 Fecal sample collection and inclusion criteria

CRC patients admitted to the Colorectal Cancer Surgery Department of the Second Affiliated Hospital of Harbin Medical University from July 2018 to February 2019 were included as study subjects. Fecal samples were collected on the day of admission and prior to any medical treatment. The inclusion criteria were as follow: (1) The patient was diagnosed as CRC by pathology and imaging; (2) The patient underwent primary tumor resection in our department by the same surgery team and underwent one-stage anastomosis (end-to-end anastomosis for rectal cancer, side-to-side anastomosis for right-sided colon cancer, and end-to-side anastomosis for left-sided colon cancer).

The exclusion criteria were as follow: (1) The patient had a history of taking corticosteroids, antibiotics or probiotics in the last 3 months; (2) The patient had a history of abdominal surgery (including open and minimally invasive surgery) and other invasive operations or treatments in the last 3 months; (3) The patient underwent colonoscopy or mechanical intestinal preparation (taking a cathartic agent, cleaning enema, etc.) in the last week; (4) The patient had a family history of tumor, a personal history of tumor, or inflammatory bowel disease; (5) The patient followed a special diet (such as vegetarians); (6) The patient had received neoadjuvant radiotherapy and chemotherapy before the operation; (7) The patient had a history of fecal bacteria transplantation; (8) The patient’s information was incomplete or they did not agree to participate in this study^[11-12].

A total of 98 patients met the inclusion and exclusion criteria. These patients represent the CRC patients in the frigid zone of northeast China. The collected fecal samples (at least 300 mg per patient) were stored in a -80°C refrigerator. All patients included in this study provided informed consent and signed an informed consent form. This study was approved by the Ethics Review Committee of the Second Affiliated Hospital of Harbin Medical University (Ethical permission code: KY2019-2008). All methods were performed in accordance with the relevant guidelines and regulations.

2.2 Propensity-score matching analysis

Among the 98 patients who met the inclusion and exclusion criteria, 10 patients suffered postoperative AL, comprising the anastomotic leakage group (AG). The remaining 88 patients without AL were designated as the normal group (NG). The diagnosis of AL was based on clinical manifestations, imaging evaluation, and laboratory examination results^[13-14]. The PSM method was used to balance the basic information between groups to exclude the differences in gut microbiota analysis results caused by mismatches in basic information. The assignment of basic information variables of the 98 patients was completed through a logistic regression model, and the propensity score was the obtained P value. A 1:1 matching was performed according to the nearest matching method, with a caliper value set at 0.2. The patient’s age, sex, American Society of Anesthesiologists (ASA) score, body mass index (BMI), tumor site, and TMN stage were selected as variables. Among these, age, BMI, sex, tumor location, and TMN stage are the basic demographic information of the patients. The ASA score, widely used in clinical work and research, is a simple and subjective quantitative score used to assess patients’ physical condition and surgical risk^[15]. After PSM, fecal samples from two groups of patients with consistent basic information were included in this study for further processing and subsequent analysis. The study flow chart is shown in Fig. 1.

Fig. 1

The study flow chart. CRC, colorectal cancer; AL, anastomotic leakage

2.3 DNA Extraction and PCR Amplification

Microbial DNA was extracted from fecal samples using the E.Z.N.A.@ soil DNA Kit (Omega Bio-tek, Norcross, GA, U.S.) according to the manufacturer’s protocols. The final DNA concentration and purification were determined using a NanoDrop 2000 UV-Vis spectrophotometer (Thermo Scientific, Wilmington, USA), and DNA quality was checked by 1% agarose gel electrophoresis. The V3-V4 hypervariable regions of the bacterial 16S rRNA gene were amplified with the following primers:

Forward: 338 (5'-ACTCCTACGGGAGGCAGCAG-3')

Reverse: 806 (5'-GGACTACHVGGGTWTCTAAT-3')

The amplification was performed using a thermocycler PCR system (GeneAmp 9700, ABI, USA). The PCR reactions were conducted using the following program: 3 minutes of denaturation at 95°C, 27 cycles of 30 seconds at 95°C, 30 seconds for annealing at 55°C, 45 seconds for elongation at 72°C, and a final extension at 72°C for 10 minutes. PCR reactions were performed in triplicate in a 20 μmol/L mixture containing 4 μmol/L of 5 × FastPfu Buffer, 2 μmol/L of 2.5 mmol/L dNTPs, 0.8 μmol/L of each primer (5 μmol/L), 0.4 μmol/L of FastPfu Polymerase, and 10 ng of template DNA. The resulting PCR products were extracted from a 2% agarose gel, further purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA), and quantified using QuantiFluor™-ST (Promega, USA) according to the manufacturer’s protocol^[11-12].

2.4 Illumina MiSeq Sequencing

Purified amplicons were pooled in equimolar concentrations and paired-end sequenced (2 × 300 bp) on an Illumina MiSeq platform (Illumina, San Diego, USA) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China).

2.5 Processing of Sequencing Data

Raw fastq files were quality-filtered using Trimmomatic and merged using FLASH with the following criteria: (1) Reads were truncated at any site with an average quality score < 20 over a 50 bp sliding window; (2) Sequences with an overlap longer than 10 bp were merged according to their overlap, allowing for no more than 2 bp mismatches; and (3) Sequences from each sample were separated according to barcodes (exact matches) and primers (allowing 2 nucleotide mismatching), and reads containing ambiguous bases were removed. Operational taxonomic units (OTUs) were clustered with a 97% similarity cut-off using UPARSE (version 7.1 http://drive5.com/uparse/) with a novel ‘greedy’ algorithm that performed chimera filtering and OTU clustering simultaneously. The taxonomy of each 16S rRNA gene sequence was analyzed by the RDP classifier algorithm (http://rdp.cme.msu.edu/) against the Silva (SSU123) 16S rRNA database with a confidence threshold of 70%.

2.6 Beta Diversity Analysis

Beta diversity analysis compares microbial community compositions to assess differences. It yields a distance matrix representing the dissimilarity between every pair of samples in the community. (Non-metric multidimensional scaling (NMDS) analysis was selected for comparing sample similarities between two groups. NMDS simplifies, analyzes, and categorizes research objects (samples or quantities) from a multi-dimensional space into low-dimensional spaces while preserving the ability to analyze raw relational data among objects. It treats the similarity or dissimilarity data between objects as a monotonic function of point distance. While maintaining the original data order relationship, the original data are replaced with new identical data columns for metric multidimensional scaling analysis.

Principal Coordinates Analysis (PCoA) is a non-constrained data dimensionality reduction analysis method that uses various distance algorithms for operations, providing better evaluation than principal component analysis based solely on Euclidean distance algorithm. Analysis of similarities (ANOSIM) is a nonparametric test used to assess whether the difference between groups is significantly greater than the difference within groups, and its test value is used for PCoA analysis.

2.7 Statistical Software

The software Mothur (version_1.30.2) was used for alpha diversity analysis. for non-metric multidimensional scaling (NMDS) analysis, Qiime (version_1.9.1) was employe to calculate the distance matrix of beta diversity. Subsequently, the R package (version_3.4.3) was used for analysis and visualization. In addition, Linear discriminant analysis Effect Size (LEfSe; http://huttenhower.sph.harvard.edu/galaxy/root?tool_id = lefse_upload) was used for multi-level species difference discriminant analysis.

2.8 Statistical Analysis

Categorical variables were compared between groups using the Chi-square test. For continuous variables, Student’s t test or the Wilcoxon rank-sum test was employed depending on whether the variables followed a normal distribution. All statistical analyses were conducted using the statistical software package SPSS 25.0 (IBM Corp, Armonk, NY, USA). A P < 0.05 was considered statistically significant, and correction of the P value was performed for the false discovery rate (FDR).

3.1 Propensity-score matching analysis

Gender, age, body mass index (BMI), American Society of Anesthesiologists (ASA) score, tumor location, and TNM stage were selected as variables for PSM (Table 1). The results indicated no statistically significant differences between the two groups in baseline data after PSM (P > 0.05). Furthermore, there were no statistically significant differences in pathological characteristics between the two groups after PSM (P > 0.05) (Table 2). With the interference of different baseline information on the patients’ gut microbiota characteristics removed, the sequencing data of 10 patients in the AG and 10 patients in the NG were further analyzed using bioinformatics.

3.2 Basic information of sequencing

After quality filtration and removal of chimeras, a total of 1,160,840 high-quality sequences were obtained, with the average sequence length concentrated in the 406 bp to 449 bp region. For further analysis, operational taxonomic units (OTU) cluster analysis was conducted for all valid sequences, and the minimum number of sample sequences was used for leveling. Finally, 322 OTUs were obtained after running the program. Subsequently, the Shannon index was plotted against the sequencing data volume to generate Shannon curves (Fig. 2). Initially, the curve rose rapidly, indicating an increase in newly discovered strains with the sequencing depth. As the sequencing depth increased, the curve gradually approached a plateau, suggesting that the sequencing depth had met the requirements, and all strains had been covered.

Fig. 2

The Shannon curves reflect sequencing depth reach the requirement

3.3 Alpha diversity analysis and community composition analysis

Alpha diversity analysis is mainly used to evaluate the diversity and richness of species composition within a sample^[16]. At the OTU level, the samples from the two groups were compared for differences in alpha diversity. The Ace index, Chao index, and Sobs index were selected to evaluate differences in species richness between the two groups, while Shannon index and Simpson index were selected to evaluate differences in species diversity. The results indicated no statistically significant difference in species composition diversity and richness between the two groups (P > 0.05) (Table 3). Subsequently, the top 15 genera in the abundance ranking were selected to construct a community heatmap at the genus level. Visual inspection of the heatmap allowed for the examination of community composition differences between the two groups (Fig. 3). It was observed that the abundance levels of each genus varied between the two groups, suggesting potential differences in community composition. Further beta diversity analysis deemed necessary to identify specific differences.

Fig. 3

The community heatmap at the genus level

3.4 Beta diversity analysis

Beta diversity analysis compares microbial community compositions to evaluate the differences between microbial communities^[17]. In this study, we adopted principal coordinate analysis (PCoA) and non-metric multidimensional scaling (NMDS) analysis^[11-12]. Two algorithms were used for PCoA in this study. The weighted_unifric algorithm (R = 0.246, P = 0.012) and unweighted_unifric algorithm (R = 0.207, P = 0.008) based on the OTU level revealed significant differences in the principal components (Fig. 4). To assess the advantages and disadvantages of NMDS analysis results, stress was calculated. The stress values for analyses at the phylum (stress = 0.024), class (stress = 0.075), order (stress = 0.081), family (stress = 0.141), genus (stress = 0.143), and species levels (stress = 0.149) were all stress below 0.2, indicating significant differences in beta diversity between the two groups (Fig. 5).

Fig. 4

Principal coordinate analysis based on operational taxonomic units (OTU) level

(A) unweighted_unifric algorithm, R = 0.207, P = 0.008; (B) weighted_unifric algorithm, R = 0.246, P = 0.012

Fig. 5

Non-metric multidimensional scale analysis on different level

(A) Phylum level, stress = 0.024; (B) Class level, stress = 0.075; (C) Order level, stress = 0.081; (D) Family level, stress = 0.141; (E) Genus level, stress = 0.143; (F) Species level, stress = 0.149

3.5 Significant test of differential microbiota between two groups

At the phylum level, Proteobacteria showed a significant difference between the AG group and the NG group (P = 0.021) (Fig.6). At the genus level, there were statistically significant differences in the richness of Streptococcus (P = 0.045), Citrobacter (P = 0.008), and Klebsiella (P = 0.002) between the AG and NG groups, with higher abundance observed in the NG group (Fig.7). To further analyze the contribution of the abundance differences of these three genera to the inter-group differences, we used LEfSe multi-level species difference discriminant analysis to generate a linear regression analysis (LDA) discriminant column (Fig. 8). A higher LDA score indicates a greater influence of the abundance of representative species on the observed difference effect. The LDA threshold was set at 3.5. The results suggested that Streptococcus, Citrobacter, and Klebsiella made the highest contributions to the differences between the groups.

Fig. 6

Differential microbiota on phylum level

Fig. 7

Differential microbiota on genus level

Fig. 8

LDA discriminant column