We recruited patients with CC from the Port of Veracruz, Mexico, in this open-label study. Inclusion criteria were: >18 years of age, self‐reported stool frequency of 3 or fewer bowel movements per week, self‐reported stool consistency of type 1-4 on the Bristol Stool Form Scale, and fulfilment of the Rome IV diagnostic criteria for functional constipation12,13.

Patients that were receiving pharmacologic treatment were asked to withdraw treatment for at least one week, so they could be evaluated. Exclusion criteria were gastrointestinal and other disorders likely to affect gut motility; prior abdominal surgery; ongoing therapy with medications known to affect gut motility; secondary causes of constipation; dyssynergic defecation; colonic inertia; anorectal problems; lactose intolerance, cow’s milk and soy allergy; regular consumption of probiotics, prebiotics, fiber supplements, and laxatives; antibiotic use within 4 weeks of study onset; and ongoing alcohol, drug, or medication abuse. The prevalence of constipation is almost twice as common in females, therefore we decided to only include women, to avoid potential confounding factors.

Patients were instructed to ingest one capsule of Alflorex® containing 1×109 CFUs/g B. infantis strain 35624 daily, for a consecutive period of 8 weeks, while continuing their regular dietary habits. B. infantis 35624 is also known as B. longum subsp. infantis (https://www.precisionbiotics.com/es/), but in addition, has been (more precisely) catalogued as B. longum subsp. longum (the original misassignment resulted from the sole use of the partial 16S rRNA-based taxonomic classification)14. The probiotic was donated by the Menarini Mexico company (Mexico City, Mexico), which only provided the drug and did not intervene in the study design or the result analysis.

Patients were informed that they were being invited to participate in a study whose aim was to investigate the effect of a probiotic on their symptoms and the composition of gut microbiota, without explicitly promising a beneficial effect. Considering the amount of product in each capsule (247 mg) and the concentration of the probiotic (1×109 CFUs/g), each capsule provided 2.47×108 CFUs. Patients filled out a questionnaire regarding clinical parameters, such as stool consistency, self-perceived bloating intensity, and days without symptoms, at the baseline and the end (final) of probiotic administration. Stool consistency was assessed by the patients using the Bristol scale15, in which 1 represents hard lumps, and 7 represents watery stools. Self-perceived bloating intensity was expressed, using a visual analogue scale from 0 to 100. Global improvement was evaluated, using the overall symptom question: in the last 7 days, do you consider that you had adequate relief from your constipation symptoms? Patients were instructed to use 17 g daily of polyethylene glycol, on a regular basis if necessary, if rescue therapy was required.

Fecal samples were collected by the patients at home, at the baseline and end of probiotic administration. The samples were placed in coolers with frozen gel packs and brought to our laboratory within one hour, where they were placed at –20°C, until processing. Because the homogenization of feces prevents the accurate evaluation of microbial communities16, we thawed the samples and obtained the aliquots for analysis, mostly from the middle of the fecal material available from each sample17. Total genomic DNA was obtained and purified from one aliquot (100 mg) from each fecal sample, using bead-beating, coupled with a commercial DNA extraction kit (Wizard® Genomic DNA Purification Kit Protocol, Promega Corporation, Madison, Wisconsin, USA). Two separate aliquots from the baseline sample from the last patient (sample ID: S19) were processed to assess the replicability of the procedure. The DNA samples (n=27, 26 from the 13 patients plus one duplicate from S19) were adjusted to 150ng/μl and shipped to Molecular Research LP (Shallowater, Texas, USA). PCR was carried out to amplify a semi-conserved region of the 16S rRNA gene using primers 515F (GTGYCAGCMGCCGCGGTAA) and 806R (GGACTACNVGGGTWTCTAAT). High-throughput DNA sequencing of the PCR amplicons was performed in a MiSeq instrument (Illumina, San Diego, California, USA), as described elsewhere18. Unless otherwise stated, only the results for the 13 sample pairs (n=26), using the first replicate of sample S19, are shown and discussed.

We used different tools in QIIME219 v.2020.2 to analyze the sequencing data, employing default parameters, unless otherwise mentioned. The sequencing laboratory provided free taxonomic assignment, using a curated database, but the results are only mentioned in the context of our taxonomic assignment in QIIME2, using the Greengenes 13.8 database. The phylogenetic method, UniFrac20, was used for comparing microbial communities, with both weighted and unweighted distances21.

It is common practice to analyze the different bacterial taxa separately. However, that approach is not necessarily better for understanding microbial fluctuations in complex ecosystems, where different types and numbers of microbes coexist. Therefore, we used the SAS® University Edition (SAS Institute, Cary, North Carolina, USA) to perform the multivariate analysis for different taxa at the same time. When needed, the Wilcoxon signed-rank test was used to compare variables (e.g., relative abundance of bacteria) between both periods of time (baseline and final). We also performed correlation analyses using the baseline and final values for the individual taxa, as well as between the different combinations of the time points and taxa (e.g., baseline values for taxa A vs. final values for taxa B), to investigate any relationship that could impact microbe variability over time. Plots were produced using SAS and PAST v.3.25.

A written statement of informed consent was required for participation in the present study. All procedures were carried out in accordance with the WHO code of ethics (Declaration of Helsinki) and approved by the local ethics committee of the Institute for Medical Biological Research (Universidad Veracruzana, 2018-013R). The authors declare that this article contains no information that can identify the patients.

Nineteen female patients were included, but only 13 were able to complete the protocol (Table 1). Prior to the study, 3 patients were not taking any medication or supplements to treat their constipation, 6 were on fiber supplements, and 11 used laxatives (6 polyethylene glycol, 3 lactulose).

Six of the 13 enrolled patients (46%) reported an improvement on the Bristol stool scale of more than 50%, compared with the baseline, and 10 of the 13 patients (77%) reported overall symptom improvement. Bloating intensity was significantly lower, after probiotic administration (67% vs. 42%, p=0.02) (Table 2), and the percentage of time (in days) without symptoms was also lower after the probiotic (49% vs. 39%), but the difference was not significant (p=0.48). Only 2 patients reported taking the rescue medication, intermittently. However, those results are likely to have limited clinical relevance due to a possible placebo effect and lack of a control group.

After demultiplexing, we obtained 2,069,625 16S sequences from 27 fecal samples (mean: 76,652 sequences/sample, min: 25,964, max: 170,946). The use of DADA222 for quality filtering yielded a total of 1,868,215 sequences (∼90% of the original), using 280 nucleotides as the truncating length.

Using a rarefaction depth of 24,000 sequences/sample, DADA2 revealed the presence of 1,123 amplicon sequence variants (ASVs) that belonged to 120 genera from 16 phyla. There were no significant differences in the number of ASVs, Shannon diversity values, Pielou’s evenness values, and Faith’s phylogenetic diversity values (Fig. 1). With the exception of Faith’s phylogenetic diversity values, there was a noticeable dissimilarity in the spread of the samples, between the baseline and final periods (Fig. 1).

Figure 1.

Boxplots with jitter showing alpha diversity metrics. There was no statistically significant difference between the baseline and final periods, but the values from most metrics showed a higher variation after probiotic use. The plots were constructed in PAST v. 3.25. ASVs: amplicon sequence variants.

(0.34MB).

A total of 99.3% of the 16S reads could be classified into a given phylum, but only 63% were assigned to a given genus. That percentage is not uncommon23, but most research articles prefer not to declare the percentage of reads assigned to the different taxonomic levels. Seven phyla accounted for ∼99% of all reads (Firmicutes: 55.8%, Bacteroidetes: 15.0%, Cyanobacteria: 10.2%, Tenericutes: 5.4%, Proteobacteria: 5.3%, Actinobacteria: 4.9%, and Verrucomicrobia: 2.0%). Those results partially agreed with the results provided by the sequencing laboratory that uses a different database for taxonomic assignments (Firmicutes: 60.6%, Bacteroidetes: 15.5%, Cyanobacteria: 10.3%, Proteobacteria: 6.5%, Actinobacteria: 4.9%, and Verrucomicrobia: 2.0%), with the exception of Tenericutes, which only accounted for 0.003% of all reads. With few exceptions, the variation in relative abundance of taxa between the two replicates of baseline sample S19 was lower, compared with the variation between the baseline and final values for all sample pairs.

The relative abundance for each phylum revealed a high inter-individual variation, as well as a wide difference in the variation over time for each subject (Fig. 2), as shown elsewhere. For instance, 7 of the 13 patients (54%) experienced a decrease in the relative abundance of Firmicutes of 37%, on average, whereas 6 of the 13 patients (46%) showed an increase in that taxon of 16%, on average (Table 2). In other words, the change in Firmicutes over time occurred in a similar proportion of patients (54% vs. 46%), but the change was more pronounced (37% decrease vs. 16% increase) in the patients that experienced a decrease in said taxon. Regardless, the difference in Firmicutes between the two time periods did not reach statistical significance (p=0.28). In the case of Bacteroidetes, 5 of the 13 patients (38%) showed a decrease in that taxon of 12%, on average, whereas 8 of the 13 patients (62%) showed an increase of 18%, on average (Table 2). The difference in Bacteroidetes between the two time points did not reach statistical significance, either (p=0.38) (Table 2). Two patients showed a small decrease in Proteobacteria (from 7.4% to 2.6%, and from 1.97% to 1.94%, respectively) after probiotic administration, whereas the other 11 patients showed an average increase of 8.5% in the taxon (Fig. 2). The difference in Proteobacteria between the baseline (3.2%) and final (10.1%) periods was statistically significant (p=0.01), although it was mostly driven by two patients (Fig. 2) and was not clearly associated with any particular taxon (Table 2).

Figure 2.

Percentages of the16S rRNA gene reads at the phylum level for all 13 patients at the baseline and the end (final) of probiotic administration. The purpose of this figure is to help visualize the high inter-individual variation (separation of dots within each time point), as well as a wide difference in the variation over time (change in the 16S reads from the baseline to the final periods) for each subject.

The plots were constructed in SAS.

(0.39MB).

The relative abundance of the genus Bifidobacterium was relatively high (the tenth most abundant genus) in our own analysis in QIIME2 (average 3.4%, min: 0.07%, max: 31.5%) and in the results from the sequencing laboratory (average 3.3%, min: 0.05%, max: 30.8%). We observed 9 ASVs that were related to Bifidobacterium, but most of them were low in abundance and showed low prevalence (Table 3).

We used the relative abundance of the 16S reads from the 6 most abundant phyla, in different combinations, to find the lowest value for the Wilks’ lambda. Accordingly, the combination of 3 phyla (Firmicutes, Proteobacteria, and Verrucomicrobia) was associated with the lowest value (supplementary table S1). This is interesting because it implies that the combination of those 3 taxa may be more precise for analyzing the differences between the two time points.

The relative abundance of Actinobacteria (R2=0.5, p=0.007) showed the strongest correlation between the baseline and final periods, but that was mainly driven by the data from one patient (supplementary table S2). The only other phylum that showed a trend for significance was Proteobacteria (R2=0.3, p=0.046, supplementary table S2). We also examined correlations between the different phyla at the different time points (supplementary table S3). Interestingly, there was an inverse correlation between the baseline values for Firmicutes and Proteobacteria (R2=0.73, p=0.0002), as recently shown in a study on milk microbiota24. Another interesting finding was the positive relationship between baseline Cyanobacteria and Proteobacteria (R2=0.82, p=2.4×10-5). Cyanobacteria are often considered contaminants, but Ley et al.25 raised the intriguing possibility that gut cyanobacteria may be descendants of non-photosynthetic ancestral cyanobacteria that adapted to life inside the digestive tract.

The use of PCoA and UPGMA clustering on unweighted UniFrac distances suggested little variation in the two replicates from S19 and the presence of two clear clusters of samples that were partly explained by probiotic administration (Figs. 3A and C). Because of the possible existence of enterotypes (clusters of samples identifiable by the variation in the presence of specific taxa and their levels26), we removed the different taxa, one by one, to investigate whether the absence of a given taxon would explain that clustering. However, we found no evidence to suggest that the clustering was due to the presence (or absence) of any particular taxon. Interestingly, the two clusters of samples disappeared, when analyzing data from separate groups of selected taxa (e.g., Firmicutes, Proteobacteria, and Verrucomicrobia), but the significance of that finding is still unclear to us. Unlike the results regarding the unweighted UniFrac distances, the analysis of the weighted UniFrac distances showed a stronger separation of replicates from S19 but revealed no significant clustering of samples (Figs. 3B and D). Regardless of any visual clustering, the PERMANOVA test did not show any significant differences in the unweighted (p=0.12) or weighted (p=0.81) UniFrac distances.

Figure 3.

PCoA plots and UPGMA trees using unweighted and weighted UniFrac distances. PCoA plots using unweighted (A) and weighted (B) distances, UPGMA trees using unweighted (C) and weighted (D) distances. In A and B, the baseline samples are represented by filled blue squares, and the final samples by a red X. The sample that was processed in duplicate (S19, baseline) is shown, using a different symbol (empty blue squares) for better visualization. The trees were constructed using 10,000 bootstrapped samples. The Euclidean similarity index was used for all plots in PAST v. 3.25.

(0.34MB).

Body mass index and age were the only continuous variables associated with the patients and did not significantly explain the number of ASVs observed. Some of the taxa were related to the number of ASVs, particularly the butyrate producer, Oscillospira (R2=0.55, p<0.0001), an enigmatic member of Firmicutes (family Ruminococcaceae) from the human gut microbiota that has been shown to be positively associated with leanness and health27 (supplementary table S4).

The concept of core microbiome refers to the species (i.e., ASVs) that are present in the majority of samples in a given environment, for example, human feces. Few operational taxonomic units (OTUs) were shared by 50% of all samples (122 of 1,123, or 10.9% of all ASVs), baseline (153, 19.1%) or final (105, 12.3%) samples, implying that most ASVs were unique to each individual.