16S rRNA Sequencing and Metagenomics Study of Gut Microbiota: Implications of BDB on Type 2 Diabetes Mellitus

To characterize the hypoglycemic effect of BDB and whether it can improve the typical symptoms of type 2 diabetes, polydipsia, polyphagia and weight loss, the levels of blood glucose, food intake, water intake and body weight of each group were measured at the end of the experiment.

First, the hypoglycemic effects of BDB were verified after a 7-week oral administration of BDB in diabetic BKS db mice. As shown in Figure 1a, diabetic BKS db mice had hyperglycemia with higher blood glucose levels compared with non-diabetic BKS mice. Oral administration of BDB effectively reduced the blood glucose levels at the 7th week. The first-line T2DM drug metformin was used as the positive control. BDB had no effects on the body weight (Figure 1b) and food intake (Figure 1c), but significantly reduced the water intake of diabetic mice (Figure 1d).

16S rRNA pyrosequencing was used to analyze the overall composition of gut microbiota. We obtained a total of 1,108,806 clean tags from 1,176,921 raw tags. The rarefaction curves showed differences in microbial diversity and rationality of samples among the BKS, BKS db, metformin, and BDB groups (Figure 2a). These sequences had an average length of 400–440 base pairs, and clustered into 1549 operational taxonomic units (OTUs) at 97% 16S rRNA gene sequences identity, including 9 phyla, 92 genera, representing them as the same source of species. Six hundred and fifty-eight of the total OTUs were shared among all the samples, and some parts of them belong to unique OTUs. This showed that over 42% of the observed OTUs were not eliminated following drug treatment and that more than half of OTUs have changed, which illustrated that there were significant differences among model and BDB-treated groups (Figure 2b). According to Chao 1 and observed species richness index, the diversity of microbiota decreased significantly in the BDB-treatment group compared with other groups (p < 0.05), which was consistent with the previous study. Zhang et al., found berberine and metformin, significantly reduced both the richness and diversity of the gut microbiota in rats [21]. But no significant difference was observed in the metformin-treatment group compared with the diabetic model group (p > 0.05) (Figure 2c–d).

We found BDB- and metformin-treatment would affect the overall composition of the gut microbiota. According to our data, the most abundant phyla identified included Bacteroidetes and Firmicutes, which consistist with other resports [22]. The other abundant phyla included Verrucomicrobia, Proteobacteria, Actinobacteria, Tenericutes, Deferribacteres, Cyanobacteria and several unclassified bacteria (Figure 3a). Compared with the group of BKS db mice, metformin treatment significantly increased the relative abundance of Bacteroidetes (p < 0.05), while BDB treatment significantly increased the phyla of Verrucomicrobia (p < 0.001). The relative abundance of Firmicutes decreased significantly in the BDB-treated group (p < 0.05). Besides, the ratio of Firmicutes to Bacteroidetes (F/B) decreased in both metformin- and BDB-treatment groups (Figure 3b). The ratio of Firmicutes to Bacteroidetes could be used as an indicator of imbalance of gut microbiota in diabetic and obese mice [23,24]. With that, we analyzed the fecal samples at the family level, which showed dramatically increasing bacteria in the BDB-treated group, such as Verrucomicrobiaceae, compared with the BKS db group (p < 0.001). The relative abundance of Bacteroidales_S24-7_group showed a similar trend after the BDB treatment. Meanwhile, Eubacteriaceae, Lactobacillaceae and Rikenellaceae showed less abundance after 7-week BDB intervention (Figure 3c). At a deeper bacterial taxa level, a greater detailed difference was observed (Figure 3d). The heatmap showed the bacteria at the genus level in different samples, Akkermansia, Lactobacillus, Alistipes, and Lachnoclostridium were increased in the BDB-treated group. Meanwhile, the relative abundance of Lachnospiraceae_NK4A136_group, Lachnospiraceae_UCG-001, Prevotellaceae_UCG-001 (p < 0.05) decreased sharply. Thus, these results had given overall profiles of gut microbiota in different treated groups, which suggested that the intake of BDB and metformin may benefit the stability of certain gut microbiota and regulate microbial imbalance.

The beta diversity metric was used to measure the extent of the similarity between microbial communities. Principal component analysis (PCA) and principal coordinates analysis (PCoA) analyses were used in the present study to observe interindividual variation and the microbiota separation. The PCA two-dimensional plot showed that the samples of the BDB-treatment groups were far from the other three groups, and basically gathered on the right side of the vertical axis, while samples of the BKS group were gathered on the left side of the vertical axis (Figure 4a), indicating BDB could drive gut microbiota change along PC2. The BKS db and metformin group had better grouping level on the horizontal axis and separated to some extent. The most significant gut microbial variations were induced by BDB which formed a separate cluster compared with other groups. The result also showed that the total diversity captured by the top two principal coordinates was 36.6%. Analysis of similarity (ANOSIM) statistical analysis of OTUs number of mice in each group showed that there were more remarkable differences in community structure between groups (p = 0.001, R = 0.6761), which indicates that the difference between groups is greater than that within groups, and the experimental group is meaningful. The R-value was greater than 0 between groups. The group composition of different samples was further studied based on unweighted (Figure 4b) and weighted ((Figure 4c) UniFrac PCoA analysis. The results showed that segregation degrees between weighted UniFrac PCoA analysis and the unweighted one were different, which revealed that the gut microbiota structure changed significantly by BDB treatment.

The differentially abundant taxon was identified by the Linear Discriminant Analysis (LDA) Effect Size (LEfSe) method using non-parametric factorial Kruskal–Wallis and Wilcoxon rank-sum test (p < 0.05). The histogram of the distribution of LDA values (LDA > 2.5) shows the statistically different biomarkers between the groups (p < 0.05) (Figure 5a). The evolutionary branch diagram shows the differential biomarkers from the phylum to the genus (Figure 5b). The BDB-treated group showed selective enrichment in 35 biomarkers, including Verrucomicrobia, Verrucomicrobiaceae, Verrucomicrobiales, and Akkermansia. The metformin-treated group showed a particular effect on Prevotellaceae UCG-001 and Prevotellaceae, Parabacteroides. Our results indicated that the gut bacteria dysbiosis could be extensively alleviated by metformin and BDB intervention, and certain genera belonging to the 3 predominant phyla may be used as identification biomarkers.

To investigate the effect of BDB on gut microbial function and metabolic pathway in diabetic BKS db mice, we extracted DNA from the feces of mice for metagenomic analysis. We obtained a total of 133,948,850 clean reads (40.2 Gb) and 754,022 open reading frames (ORFs). The Cluster of Orthologous Groups of proteins(COG) and Kyoto Encyclopedia of Genes and Genomes(KEGG) databases were used to annotate the processed data. LEfSe (LDA Effect Size) analysis was used to detect the functional COG categories with significantly different abundances between the BDB-treated mice and diabetic model mice. As shown in Figure 6a, 13 functional COG categories were observed with significantly overabundant reads in the BDB-treated group, which were assigned to the categories of RNA processing and modification (A), Nucleotide transport and metabolism (F), and Lipid transport and metabolism (I), etc. By contrast, the BKS db group had more reads involved in the Cytoskeleton (Z), Transcription (K), Cell motility (N), Defense mechanisms (V), Signal transduction mechanisms (T) and Replication recombination and repair (L) categories. Together, information storage and processing (cluster I) and cellular processes and signaling (cluster II) were summarized in diabetic model mice, and metabolism cluster (cluster III) was predominant in BDB-treated mice.

Furthermore, the KEGG pathway database was used to analyze changes in pathway function between BDB-treated mice and BKS db mice. The significant differences in the KEGG pathway were analyzed by LEfSe. Using the threshold values of LDA > 2.5 and p < 0.05, we found that, at level 2, the functional categories related to carbohydrate metabolism, metabolism of other amino acids, drug resistance, lipid metabolism, amino acid metabolism, and transport and catabolism were enriched in the fecal microbiota of BDB-treated mice (Figure 6b), these pathways are generally considered to be important determinants of gut microbiota composition and health [25]. At level 3, 23 KEGG pathways (including citrate cycle, carbon metabolism, valine-leucine-isoleucine degradation, and pyruvate metabolism, and others) and 15 KEGG pathways (including two-component system, glucagon signaling pathway, and bacterial chemotaxis, etc.) were significantly enriched between BDB-treated mice and BKS db mice, respectively (Figure 6c). Compared with BDB-treated group, we found the motility ability of bacteria increased in db/db group (based on the pathway of bacteria motility and flagellar assembly increased significantly). In addition, BDB may alleviate the metabolic disorder of diabetic mice by inhibiting the starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism. Interestingly, multiple representative functional pathways in BDB-treated mice were also involved in metabolism. According to the database of KEGG, the significantly different functional pathways were similar to that in COG categories. Briefly, genes related to the basic metabolism were enriched in the BDB-treated group compared with the BKS db group.

BDB can promote the abundance of short-chain fatty acid-producing bacteria in vivo, which was mentioned earlier in the article. In order to explore whether BDB can directly affect the intestinal flora and its impact on the growth of gut microbiota, we selected significantly increased key communities by BDB to test the interaction between drugs and bacteria in vitro, including A. muciniphila and L. johnsonii. E. coli, which is one of the most common opportunistic pathogens in the intestinal tract and can damage the intestinal mucosal barrier, was also examined. The result showed that BDB tended to promote the growth of A. muciniphila and L. johnsonii and inhibit the growth of E. coli, which indicates the role of promoting and inhibiting proliferation for some special bacterial species (Figure 7a–c). However, the effect of BDB on the growth of A. muciniphila and L. johnsonii is not significant, which may be due to in vitro culture cannot simulate the environmental conditions of the intestinal tract in vivo. The effects of BDB were probably taken by the inhibition of potential pathogenic/harmful bacteria. BDB-sensitive cells decrease, the proportion of BDB-resistant species, such as A. muciniphila and L. johnsonii increase.

Male wild type BKS mice (C57BLKS/J, the Jackson Laboratory stock number 000662) and BKS db mice (BKS.Cg^–Dock7^m+/+Leprdb/J, the Jackson Lab stock number 000642), aged six- to eight-weeks, were supplied by the Model Animal Research Center of Nanjing University (MARC). The mice were provided with individually ventilated cages in a specific pathogen free (SPF) animal rooms and were maintained temperature (24 ± 2 °C), humidity (60–80%), under 12 h light/dark cycle. All mice were fed sterile water and a normal chow diet ad libitum.

The mice for this study were administered according to the guidelines of protocol HAIFAJIZI-2013-3 (approval date: 9 December 2013) approved by the Institute of Oceanology committees for care and use of laboratory animals. A week after the mice adapted to the environment, all diabetic BKS db mice which accord with the experimental requirements were randomly divided into 3 groups: diabetic model group (db group), metformin-treated group (Met group, 100 mg.kg⁻¹.body⁻¹.day⁻¹), BDB-treated group (BDB group, 100 mg.kg⁻¹.body⁻¹.day⁻¹). Age-matched male wild type BKS mice were used as the normal control (BKS group). After 7 weeks, we measured the body weights, food intake, water intake and fasting glucose levels using an Accu-Chek Performa glucometer (Roche, Mannheim, Germany).

Fecal samples (n = 6) were collected under sterile conditions and snap-frozen in liquid nitrogen for half an hour and kept at −80 °C until used in 16S rRNA and metagenomics analysis.

Bacterial genomic DNA extraction using the DNA mini kit (Omega, Norcross, GA, USA) from 200 mg frozen fecal. Then the V3-V4 variable regions of bacterial 16S rRNA gene were PCR-amplified using the primers (338F 5′-ACTCCTACGGGAGGCAGCAG-3′, 806R 5′-GGACTACHVGGGTWTCTAAT-3′) and the gene of 16S rRNA was sequenced using the MiseqPE300 platform (service provided by Allwegene Technology Inc (Beijing, China). Raw Fastq files were quality filtering using Trimmomatic (v0.36), Pear (v0.9.6), Flash (v1.20) and Vsearch (v2.7.1). Vsearch was used to gain operational taxonomic units (OTUs) according to 97% sequence similarity and remove the chimeric sequences (2.7.1 https://github.com/torognes/vsearch↗). The taxonomy of OTUs was analyzed by RDP classifier algorithm (http://rdp.cme.msu.edu/↗) against the Silva (release 123) 16S rRNA database using confidence threshold of 70%. The relative proportion of each OTU was annotated at the Phylum, Class, Order, Family, Genus and Species levels.

The metagenomic sequencing method was used to investigate functional changes in the gut microbiota. A total of 6 fecal samples to extract the genomic DNA (3 samples from BDB-treated BKS db mice and others from placebo-treated diabetic mice).

The qualified DNA samples were randomly interrupted with the size of about 270 bp by covaris ultrasonic crusher. The metagenomic sequencing was carried out on an Illumina HiSeq sequencing platform (Illumina Inc., San Diego, CA, USA) by Allwegene Technology Inc (Beijing, China). The raw sequence reads with quality scores below 20 accounts for more than 50% of the total read length were trimmed. We filtered the low-quality sequences of raw reads obtained by sequencing, and the subsequent analysis was based on filtered clean reads. For the host contamination of the sample, compared with the host database to filter out the reads that may come from the host. The open reading frames (ORFs) were predicted by MetaGeneMark and removed redundant by using CD-HIT software to obtain the non-redundant-nucl (https://github.com/weizhongli/cdhit↗). The COG functional categories were annotated by the EggNOG database via Smith-Waterman comparison algorithm. The KEGG pathway annotation was performed using a BLAST search by the KEGG database.

Alpha, beta diversity, principal component analysis (PCA) and principal coordinates analysis (PCoA) phylogenetic distance matrix were conducted by the Quantitative Insights into Microbial Ecology (QIIME). The relative abundance of biomarkers among the mice of different groups and the difference of function analyzed with Linear Discriminant Analysis (LDA) Effect Size (LEfSe) method. Only when the p-value is less than 0.05, and an LDA value greater than 2.5 is considered to be significantly rich.

Lactobacillus johnsonii isolated from faeces of the mice, Akkermansia muciniphila (DSM22959) was purchased from the BeNa culture collection and Escherichia coli were provided by the first institute of oceanography, SOA. They were inoculated in MRS broth medium (de man rogosa Sharpe; Hopebiol, Qingdao, China), modified brain and heart infusion broth medium containing (in g/liter) yeast extract (5), cellobiose (1), maltose (1), cysteine (0.5), mucin (1)and hemin (0.01), and LB broth medium (Hopebiol, Qingdao, China) respectively at a proportion of 1%, and BDB was added to make the final concentration of 10μM in the anaerobic incubator for culture (37 °C, 5% hydrogen, 5% carbon dioxide and 90% nitrogen). OD₆₀₀ was measured continuously at a certain interval, and its growth curve was analyzed.

All the experimental data were presented as the mean ± standard error of the mean (SEM). One-way analysis of variance (ANOVA) followed by Tukey’s post hoc test was used to evaluate the difference between groups. A p-value < 0.05 was considered statistically significant. All statistical analyses were determined by SPSS 22.0.