A total of 30 subjects met the Rome IV criteria for STC and 30 healthy controls were recruited from the outpatient department of Tianjin People’s Hospital from April 2020 through December 2020 (Drossman and Hasler, 2016; Rao et al., 2016). All subjects were evaluated by a physician. Inclusion criteria: STC patients confirmed by colonic transit stud with colonic transit time >48 h (Diamant et al., 1999); disease duration ≥ 3 months; age ≥ 18 years; body mass index (BMI): 18.5–25 kg/m². Exclusion criteria: congenital megacolon; constipation caused by secondary intervention (e.g., drugs, metabolic disorder, endocrine disorders, or neurological disorders); previous abdominal surgery and perianal surgery; history or current of gastrointestinal diseases (e.g., malignancy and inflammatory bowel disease); infected with enteric pathogens; used prebiotics, probiotics, or proton pump inhibitors within the past month; pregnant or breast-feeding women; long-term smoking and/or alcohol addiction; confirmed to have hepatic, renal, cardiovascular, respiratory or psychiatric disease; suffered from other disease that could affect intestinal transit and the gut microbiota.

The researchers obtained metadata about each participant, including sex, age, gastrointestinal symptoms, dietary supplements, medications, and allergy history. All participants did not take antibiotics, probiotics, and prebiotics in the 3 months before the fecal sample collection. The quality of life of all subjects was evaluated according to the Wexner constipation score standard and the Gastrointestinal Quality of Life Index (GIQLI) (Eypasch et al., 1995; Agachan et al., 1996). All subjects or their families signed informed consents.

Standardized instructions and kits for collecting stool were provided to subjects, and the feces were collected and frozen in liquid nitrogen for cryopreservation immediately. The feces of the subjects were scored according to the Bristol stool form scale (BSFS) (Chumpitazi et al., 2016).

Total genome DNA from fecal samples was extracted using the CTAB/SDS method. DNA concentration and purity were monitored on 1% agarose gels. According to the concentration, DNA was diluted to 1 ng/μl using sterile water for further use. The V3–V4 variable region of 16S rRNA in fecal samples was amplified by PCR using a universal primer designed for the conserved region. According to the characteristics of the amplified 16S region, Illumina Miseq sequencing platform was used to construct a small fragment library and perform paired-end sequencing. Through the splicing and filtering of reads, cluster of operational taxonomic units (OTU) with Greengenes database (http://greengenes.lbl.gov/), species annotation, and abundance analysis, the species composition of the sample was revealed. Sequence analysis was performed by UPARSE software package using the UPARSE-OTU and UPARSE-OUTref algorithms. In-house Perl scripts were used to analyze α-diversity (within samples) and β-diversity (among samples). Sequences with ≥97% similarity were assigned to the same OTU. A representative sequence from each OTU was annotated taxonomic information using the RDP classifier (Wang et al., 2007).

We rarified the out table and drew four curves to evaluate the sequencing data including rarefaction curves, Shannon curves, rank-abundance curves, and species accumulation curves. The differences in α-diversity index between groups were analyzed by seven metrics: Observed Species, Shannon, Simpson, Chao 1, ACE, Coverage, and PD whole tree. Shannon index estimates flora diversity, and Chao1 estimates the species abundance. Principal component analysis (PCA) was applied to reduce the dimension of the original variables using QIIME (Bolyen et al., 2019). QIIME calculates both weighted and unweighted UniFrac distance, which are phylogenetic measures of β-diversity (Caporaso et al., 2010). The unweighted UniFrac distance was used for principal coordinate analysis (PCoA) to get principal coordinates and visualize them from complex, multidimensional data. At the same time, non-metric multi-dimensional scaling (NMDS) was also used to perform rank ordering and PCoA on the OTU data (Rivas et al., 2013). To confirm differences in the abundances of individual taxonomy between the two groups, STAMP software was utilized. Linear discriminant analysis effect size (LEfSe) was used for the quantitative analysis of differential microbes within STC and healthy subject groups (Segata et al., 2011). This method was designed to analyze data in which the number of species is much higher than the number of samples and to provide biological class explanations to establish statistical significance, biological consistency, and effect-size estimation of predicted biomarkers. Finally, PICRUSt software was used to predict the functional gene composition in samples, so as to analyze the functional differences between different samples or groups.

In order to understand the correlations between different genera or species, a co-occurrence network was constructed based on the 16S rRNA data (Wang et al., 2018; Dan et al., 2020). The Spearman’s correlation coefficient was used to analyze the bacterial correlations in the STC and C groups according to the relative abundance of each species or genus, respectively. The significant correlations were visualized by Cytoscape version 3.7.1 (http://www.cytoscape.org) (Shannon et al., 2003).

For fecal metabolomics, the fecal samples were slowly thawed at 4°C. An appropriate amount of samples was added to 800 μl of pre-cooled methanol/acetonitrile/water (2:2:1, v/v) solvent, vortexed to mix, sonicated at low temperature for 30 min, and stood at −20°C for 10 min. Subsequently, the samples were centrifuged at 14,000 g for 20 min at 4°C. The supernatant was dried in a vacuum centrifuge. For UPLC-Q-TOF-MS analysis, 100 μl of acetonitrile/water (1:1, v/v) solvent was used to dissolve the dry substances. The supernatant was taken for analysis, after being vortexed and centrifuged at 14,000 g for 15 min at 4°C. The samples were separated by UHPLC with Agilent 1290 Infinity LC (Agilent Technologies, Santa Clara, CA, USA). The 6600 Triple-TOF mass spectrometer (AB Sciex) was used to collect the first- and second-order spectra of the samples. The metabolite structure identification, data preprocessing, experimental data quality evaluation, and data analysis were performed in turn according to the peak area. Specific chromatography–mass spectrometry analysis conditions were carried out in accordance with the experimental procedures of Shanghai Applied Protein Technology (Shanghai, China).

Based on the fold change (FC) analysis and t-test/non-parametric test analysis, the difference analysis of all metabolites detected in the negative ion mode was performed.

After normalizing to total peak intensity, the processed data were analyzed by the R package, where it was subjected to orthogonal partial least-squares discriminant analysis (OPLS-DA) (Worley and Powers, 2013). The sevenfold cross-validation and response permutation testing were used to evaluate the robustness of the model. The variable importance in the projection (VIP) value of each variable in the OPLS-DA model was calculated to indicate its contribution to the classification. Metabolites with the VIP value >1 were further applied to Student’s t-test at univariate level to measure the significance of each metabolite; p-value less than 0.05 was considered as statistically significant. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed on the screened differential metabolites with VIP (VIP value ≥ 1 and p-value < 0.05).

The receiver operating characteristic (ROC) curve and area under the curve (AUC) index were applied to optimize and verify the identified biomarkers (including gut microbiota and/or metabolites) and determine whether these biomarkers have diagnostic significance for STC (Falony et al., 2016; Liu et al., 2020). Random forest was used to build the diagnostic prediction models using the MetaboAnalyst 5.0 (https://www.metaboanalyst.ca) (Pang et al., 2021). To determine the association between Wexner constipation score, GIQLI, BSFS score, potential diagnostic gut microbiota, and metabolites, a correlation analysis was constructed using Spearman’s correlations.

Tissue samples of the control group were obtained from normal para-cancer colon tissue (more than 3 cm from the tumor) from colon cancer patients undergoing total colon resection. The tissue samples were confirmed to be normal colon tissue by two pathologic examiners. Colon tissues of the STC group were derived from patients who had to undergo surgery. All subjects or their families signed informed consents. All sample collection was performed in the operating room of Tianjin People’s Hospital and the samples were immediately stored at −80°C.

Paraffin sections (5 μm thick) of colon tissues from STC patients or healthy subjects were prepared after immersion in 4% paraformaldehyde for 4 h and 70% ethanol. The brief steps of HE staining included dewaxing the sections, impregnating the cytoplasm with eosin, redyeing the nucleus with hematoxylin, and dehydrating and sealing the tablets. The first step of immunohistochemistry was dewaxing and hydration. Permeabilize tissue cells by incubating with 0.1% Triton X-100 in phosphate buffered saline for 15 min on ice. Antigens of Takeda G protein-coupled receptor (TGR5, Lot: ab72608, Abcam, England) and nuclear farnesoid X receptor (FXR, Lot: bs-12867R, Bioss, China) were unmasked by microwaving sections in 10 mmol/L citrate buffer, pH 6.0 for 15 min. After blocking nonspecific proteins and endogenous peroxidase, the primary antibodies (1:200) and fluorescein-conjugated secondary antibody (1:100) were respectively incubated in a dark humidity chamber. Observe the results using a microscope in a dark room. Six to eight different fields were randomly selected from each group, and the number of TGR5- or FXR-positive cells was analyzed and calculated by ImageJ software (National Institutes of Health, United States). The lysate was used to extract the RNA from the colon tissues of two groups [Lot: DP431, Tiangen Biochemical Technology (Beijing) Co., Ltd.], and the RNA purity and concentration were detected by a spectrophotometer. cDNA was synthesized from 1 mg of total RNA using a reverse transcriptase kit [Lot: KR106, Tiangen Biochemical Technology (Beijing) Co., Ltd.]. The relative gene expression levels of TGR5 and FXR were calculated by the 2^-△△CT method using the qRT-PCR method with GAPDH as the internal reference. Primer sequences are outlined in Figure 5B. The extraction of protein of the colon tissues was performed according to the kit instructions, and the protein concentration was determined using BCA assay. The quantified protein was pretreated and separated by SDS-PAGE, and transferred to polyvinylidene difluoride membranes. The membranes were blocked with bull serum albumin for 2 h at room temperature and incubated with rabbit antibodies against TGR5, FXR, or GAPDH (Lot: GB12002, Servicebio, China)(1:2,000) overnight at 4°C. The membranes were washed and further incubated with HRP-conjugated anti-rabbit or anti-mouse secondary antibody (1:10,000 dilution) for 1 h. The protein bands were developed using enhanced chemiluminescence reagent, imaged with a gel imaging system (Tanon, China), and quantified by NIH ImageJ software.

PCoA, Fisher’s exact tests, Spearman correlations, and Wilcoxon rank-sum test were used to identify differences and the relations within the clinical characteristics. For normally distributed data, Student’s t-tests were performed to evaluate the differences in taxonomic abundance, Wilcoxon rank-sum test was performed when data were not normally distributed. p-value < 0.05 was considered statistically significant. The Benjamini–Hochberg method was used to correct p-value as false discovery rate (FDR) (Green and Diggle, 2007).

A total of 30 subjects with a clinical diagnosis of STC (age, 53.2 ± 16.323; sex, male:female, 7:23) and 30 healthy individuals (age, 47.967 ± 13.803; sex, male:female 9:21) were recruited (Table 1). There were no statistical difference in age, gender, and BMI between the two groups.

The STC group and the healthy volunteer group displayed 43,697 and 24,730 unique OTUs, respectively. A total of 24,431 OTUs were shared by both groups (Figure 1A). The rarefaction curves, Shannon curves, rank-abundance curves, and species accumulation curves were generated from OTUs, with 97% identity achieved in all samples (Figure S1). This indicated that the testing samples were sufficient and the amount of data were reasonable for the investigation of fecal microbiota. Taking the genus level as an example, there was a difference in the relative abundance of species between the STC group and the healthy group. Bacteroides was significantly increased, while Agathobacter and Subdoligranulum were decreased in the STC group (p-value < 0.05 or 0.01) (Figure 1B). The results of 7 indexes in α-diversity analysis are shown in Table 2 and Figures 1C, D. It is indicated that the α-diversity of gut microbiota in STC patients was richer than healthy subjects, but there were no significant differences except for the Simpson index (p-value = 0.021). The community composition structure and aggregation similarity of two groups were different based on the PCoA and NMDS of their microbiota (Figures 1E, F). The analysis results of LEfSe included cladograms (phylogenetic distribution) (Figure 1G) and histogram of LDA value distribution (LDA>2, Figure 1H). The LEfSe analysis results based on the prediction of KEGG pathway are shown in Figure 1I. The red and green nodes in the branches indicate the functional items and important microbial groups that play vital roles in the healthy or STC group. The differential microflora in the STC group was closely related to carbohydrate metabolism, glycan biosynthesis and metabolism, immune system, energy metabolism, digestive system, metabolic diseases, excretory system, nervous system, transport and catabolism, metabolism of terpenoid and polyketides, and lipid metabolism. There were significant differences in membrane transport (p = 0.037), cell motility (p = 0.047), and digestive system (p = 0.014) between two groups.

To find and test the potential diagnostic gut microbiota biomarkers, a random forest model was applied based on the differential genus with relative abundance > 0 in at least 95% samples of two groups. With reference to the results of LEfSe analysis of different gut microbiota, we finally selected 14 different bacterial genera or families for performing the multivariate ROC curve based exploratory analysis under the automated important feature identification and performance evaluation (Figure 2A). After comprehensive consideration of the value of the AUC, five potential diagnostic gut microbiota biomarkers (g:Bacteroides, g:Parabacteroides, f:Desulfovibrionaceae, g:Ruminiclostridium 5, and g:Subdoligranulum) were screened out and were attempted to be constructed as a diagnostic model with AUC = 0.785 (Figure 2B). The relative abundances of these five bacterial genera in each group are shown in Figure 2C. In addition to g:Subdoligranulum, the other four diagnostic gut microbiota biomarkers were significantly upregulated in the STC group.

The human intestine tract has a complex microbial ecosystem, and the influence of individual groups in different microbial communities on the intestine tract is uneven (Yilmaz et al., 2019). In order to describe the potential relationship between the flora in the gut microbial community, we further constructed a co-occurrence network of five potential diagnostic gut microbiota biomarkers in the control and STC groups based on the significant Spearman correlations (Figure 2D). A different co-occurrence network was displayed between the five potential diagnostic markers. The variation trend relationship between them was visualized. There were 10 edges in total, whose density and colors represented the closeness of the relationship. The size of nodes represented its importance in the five bacterial communities.

Based on univariate analysis, differential analysis of all metabolites detected in negative ion mode was performed. The different metabolites with FC > 1.5 or FC < 0.67 and p-value < 0.05 were visually displayed in the form of a volcano plot (Figure 3A). There were more downregulated differential metabolites in the STC group than in the control group. The score plot of the OPLS-DA score showed that there was a large variation between the two groups (Figure 3B), and this evaluation model was stable with Q ² = 0.311 (Figure 3C). VIP obtained from the OPLS-DA model could be used to measure the influent intensity and explanatory ability of the expression pattern of each metabolite on the classification and discrimination of each group of samples. Generally, metabolites with VIP value ≥ 1 and p-value < 0.05 were considered as differential metabolites to have a significant contribution to model interpretation. The information of the 21 differential metabolites selected is displayed in Table 3.

In order to clarify the functional changes related to specific differential metabolites in the stool of STC patients, we performed the KEGG pathway on the 21 screened metabolites (Figure 3D). The results showed that the KEGG pathway was mainly concentrated in secondary bile acid (BA) biosynthesis (−log₁₀ p-value = 6.420), primary BA biosynthesis (−log₁₀ p-value = 5.825), bile secretion (−log₁₀ p-value = 4.266), cholesterol metabolism (−log₁₀ p-value = 3.085), linoleic acid metabolism (−log₁₀ p-value = 2.181), and biosynthesis of unsaturated fatty acids (−log₁₀ p-value = 1.381).

The gut microbiota can modulate intestinal motility through the release of short-chain fatty acids (SCFAs) (Barbara et al., 2005; Martin-Gallausiaux et al., 2021). SCFAs produced by gut bacteria, especially butyrate and propionic acid, have multiple beneficial effects on host health, including the maintenance of mucosal integrity, the prevention of symbiotic expansion of potentially pathogenic bacteria in the gut and the regulation of energy metabolism through the gut–brain axis (Wang et al., 2019; Zhuang et al., 2019). The quantitative comparison analysis (response strength, means ± SD) of propionic acid and isobutyric acid in the feces of healthy subjects was 1,221,545 ± 1,999,220 and 1,563,167± 969,880, while in STC patients, they were reduced to 255,230 ± 550,046 and 1,814,435± 1,329,376 (p-value = 0.013 and 0.406). This indicated that the levels of SCFA changed to some extent, especially the decrease of propionic acid during the occurrence of STC, which was consistent with the previous results (Li et al., 2020; Tian et al., 2021).

We also compared the differences in the content of other secondary BAs or its conjugated BAs. In STC patients, except for lithocholic acid, all of them were reduced to varying degrees compared with healthy subjects, but the statistical differences were not significant (Table 4). Combined with differential metabolites, it was indeed found that primary and secondary BAs were significantly reduced in the feces of STC patients. Based on the above results and the roles of BAs in lipid metabolism (Fiorucci et al., 2021), we selected nine different metabolites that were deeply involved in BA synthesis, metabolism, secretion, and lipid metabolism, and demonstrated their abundance in the two groups (Figure 3G).

To find the potential diagnostic metabolic biomarkers, the random-forest model was used to evaluate the accuracy of 21 differential metabolite abundances for the classification performance anmong STC patients and healthy subjects (Figure 3E). The AUC and 95% confidence intervals of the top 12 differential metabolites are shown in Figure 3F, which were regarded as potential diagnostic biomarker metabolites. Then, to create a biomarker model, the multivariate ROC curve-based exploratory analysis was performed. It is concluded that the prediction model containing nine differential metabolites [inculding taurocholate, taurochenodeoxycholate, taurine, deoxycholic acid, cyclohexylsulfamate, cholic acid, chenodeoxycholate, arachidonic acid (peroxide free), and 4-pyridoxic acid] deeply involved in BA synthesis, metabolism, secretion, and lipid metabolism showed a high discriminatory power to predict STC status with AUC = 0.81 (Figure 3H).

In healthy subjects, the epithelium of the colonic mucosa was intact. The glands were abundant and neatly arranged, and there was no obvious cytoplasmic expansion and inflammatory cell infiltration in the interstitium. The epithelium of the mucosal layer of STC patients was damaged or missing, with the number of glands reducing, whose cytoplasm was obviously enlarged. Compared with normal tissues, more inflammatory cells infiltrating were seen in the interstitium (Figure 4A).TGR5 is a cell membrane receptor located on human chromosome 2q35, and its homology is highly conserved in humans and mammals (>80%). TGR5 is highly expressed in immune cells, intestinal tract, and gallbladder (Duboc et al., 2014). FXR is a member of the nuclear receptor superfamily and a ligand-dependent transcription factor. It is mainly expressed in liver, intestine, kidney, and other tissues, which can regulate the metabolism and enterohepatic circulation of BAs (Fiorucci et al., 2010). The results of immunohistochemistry, qRT-PCR and WB experiments showed that compared with healthy controls, the gene and protein expression levels of TGR5 and FXR in the colon tissues of STC patients were significantly reduced (Figures 4B and 5).

Correlation analysis can help measure the closeness of significant diagnostic gut biomarkers (microbiota and metabolites) to clinical manifestations. The results are shown as heatmap in Figure 6. Red indicated positive correlation, green indicated negative correlation, and white indicated non-significant correlation. The color depth was related to the absolute value of the correlation coefficient. The significance of the correlation was related to the size of the point. The smaller the p-value was, the higher the significance was and the larger the point was. Compared with microbiota, Wexner constipation score, GIQLI, and BSFS score were more correlated with diagnostic gut metabolites. Wexner constipation score was negatively correlated with these diagnostic metabolites, while the other two were positively correlated. There was obvious positive correlation among metabolites. On the whole, there was a certain negative correlation between the differential microbiota and the differential metabolites, which provided enlightenment for us to study the regulatory relationship between the two. Comprehensive use of 14 diagnostic gut microbiota and metabolite biomarkers may accurately distinguish STC and healthy population with AUC = 0.877.