Gut microbiota of obese and diabetic Thai subjects and interplay with dietary habits and blood profiles

Ethics statement

This study was approved by the Ethics Committee of Mae Fah Luang University (Ethics license: REH60075). The subjects were informed about the scope of the research project a day before participation using Thai-version information sheets. Written informed consents were obtained from all participants before sample and questionnaire collections.

Study subjects

The study included 30 subjects from Chiang Rai province located in Northern Thailand, which were considered to be representative of Thai population. Subject recruitment was conducted by voluntary participation through community clinics in the province and was carried out in July 2017. Subjects were divided into two major groups including diabetics (T2DM) and non-diabetics. Of these, seven subjects were placed in the T2DM group based on their morning fasting blood sugar level (cutoff level of > 126 mg/dL) (Reinauer et al., 2002) and irrespective of BMI. Twenty-three (non-diabetic) subjects were classified according to BMI using the criteria set by the World Health Organization Western Pacific Region (WHO, 2004), as follows: underweight (BMI < 18.5), normal or lean (18.5 ≤ BMI < 24.9), overweight (25.0 ≤ BMI < 29.9), and obese (BMI ≥ 30). Voluntary samples were taken from each group to meet a quota as follows: seven for T2DM (female (n = 6), male (n = 1)); eight for lean (female (n = 6), male (n = 1)), for overweight 8 (female (n = 6), male (n = 1)) for overweight, and seven for obese (female (n = 3), male (n = 4) and 4 men). Subjects that reported use of antibiotics (a duration of 6 months) and/or experienced diarrhea (a duration of 1 month) were excluded from the study. Average characteristics of the subjects that participated in this study are shown in Table 1. Statistical significance of each characteristic between groups (except gender) were assessed by One-way ANOVA test, followed by a post-hoc test for unequal sample size (Tukey–Kramer at a confidence interval of 0.95). The Fisher’s exact test was applied for a gender variable. The statistical analysis was performed using an R software package (stat) version 3.6.1 which Benjamini–Hochberg procedure was applied for multiple-test correction using multcomp package (version 1.4-10). Blood profiles, high-density lipoprotein (HDL) cholesterol and fasting glucose levels were also selected for further multivariate analysis.

Food frequency questionnaire

Dietary intake variables were collected using food frequency questionnaire (FFQ). The questionnaire contained 25 items of different food types including rice vermicelli, a traditional food of northern Thailand (La-ongkham et al., 2015). Records of frequency of consumption of yogurt/cheese/fermented milk, and fruits were missing for one subject (lean group). Frequencies were categorized into the following six levels: every day, 5–6 days a week, 3–4 days a week, 1–2 days a week, less than once a week, and never. The statistical significance of differences in the mean ranks among groups was determined using Kruskal–Wallis rank sum test with post-hoc analysis (Dunn’s test of multiple comparisons, p-value adjusted with the Benjamini–Hochberg method (hereafter referred to as q-value)). The frequencies of dietary consumption of each group are summarized in Table S1.

Fecal sample collection and DNA extraction

Fecal samples of all volunteers were collected in a sterilized container and immediately stored at −20 °C until further use. Total genomic DNA from fecal samples was extracted using the innuPREP Stool DNA Kit (Analytik Jena Biometra, Jena, Germany) following the manufacturer’s guidelines. Concentration and purity of DNA were evaluated on 1% agarose gels. Spectrophotometry was applied to determine the DNA concentration (ng/µl) by the Take 3 Micro-Volume Plate (Biotek, Winooski, VT, USA). Total DNA per gram of fecal wet weight was calculated and recorded.

Amplicon generation, library preparation and sequencing

The hypervariable region V3–V4 of the 16S rRNA gene was amplified using specific primers (16S V3–V4: 341F: 5′-CCTAYGGGRBGCASCAG-3′, 806R: 5′-GGACTACNNGGGTATCTAAT-3′) (Klindworth et al., 2013) with the barcode. All PCR reactions were carried out using Phusion® High-Fidelity PCR Master Mix (New England Biolabs, Ipswich, MA, USA). PCR products were run using electrophoresis on a 2% agarose gel for detection. Samples that showed a band between 400 and 450 bp were chosen for further experiments. PCR products were mixed in equidensity ratios. Then, mixture PCR products were purified with Qiagen Gel Extraction Kit (Qiagen, Germany). Sequencing libraries were generated using NEBNext ® Ultra DNA Library Pre Kit for Illumina, following manufacturer’s recommendations and index codes were added. Library quality was assessed using Qubit@ 2.0 Fluorometer (Thermo Scientific, Waltham, MA, USA) and Agilent Bioanalyzer 2100 system. Libraries were sequenced on the Miseq platform (Illumina, San Diego, CA, USA) at Novogene (Beijing, China) during September 2017 and 250 bp paired-end reads were generated. Further downstream steps included data analysis using Qiime (version 1.7.0), OTU clustering and taxa annotation, alpha and beta diversity analysis.

Data analysis

Paired-end reads were assigned to samples based on their unique barcode and truncated by cutting off the barcode and primer sequence. Paired-end reads were merged using FLASH (V1.2.7, http://ccb.jhu.edu/software/FLASH/) (Magoč & Salzberg, 2011). Splicing sequences were called raw tags. Quality filtering on the raw tags was performed under specific filtering conditions to obtain high-quality clean tags (Bokulich et al., 2013) according to the QIIME (version 1.7.0, http://qiime.org/index.html) (Caporaso et al., 2012) quality-controlled process. Tags were compared with the reference database (Gold database, http://drive5.com/uchime/uchime_download.html) using UCHIME algorithm (UCHIME Algorithm, http://www.drive5.com/usearch/manual/uchime_algo.html) (Edgar et al., 2011) to detect chimera sequences, all of which were removed (Haas et al., 2011). The raw sequence data is available at NCBI SRA with BioProject accession number PRJNA610672 (BioSample accession numbers SAMN14309526–SAMN14309555).

OTU cluster and species annotation

Sequence analysis was performed using Uparse software (Uparse version 1.0.1001, http://drive5.com/uparse/) (Edgar, 2013). Sequences with ≥97% similarity were assigned to the same OTU. A representative sequence for each OTU was screened for further annotation. Sequences were queried against the Greengenes Database version gg.13.5 (http://greengenes.lbl.gov/) (DeSantis et al., 2006; Wang et al., 2007) to obtain taxonomic information. Newly generated OTUs were aligned using MUSCLE software (Version 3.8.31, http://www.drive5.com/muscle/) (Edgar, 2004) and phylogenetic trees were generated. The OTU annotation tree was visualized using a custom R package (developed by Novogene Co., Ltd., Beijing, China). OTU abundance information was obtained by normalizing the sequence number corresponding to the sample with the least sequences (OTU counts rarefied to 103,744 reads per sample). Subsequent analysis of alpha diversity and beta diversity were all performed basing on this output normalized data. The relative abundance of gut bacteria between sample groups was compared by the unpaired two-samples Wilcoxon test and multiple comparisons were adjusted with Benjamini–Hochberg method (q < 0.05) using R software package (stats) version 3.6.1.

Alpha diversity

Alpha diversity was applied to analyze complexity of species diversity for each sample using the following indices: ACE, Chao1, observed-species, Shannon, Simpson and Good’s coverage. The unpaired two-samples Wilcoxon test with Benjamini–Hochberg procedure (q < 0.05) was used to compare alpha diversity indices between groups for statistical differences. Comparisons were visualized as a box plot by ggplot2 (Wickham, 2009) in R software (version 3.6.1).

Beta diversity analysis

Beta diversity analysis was used to evaluate differences of fecal samples in bacterial community structure between BMI and T2DM groups. Beta diversity was calculated using both weighted and unweighted unifrac. Principal Coordinate Analysis (PCoA) was performed to get principal coordinates and visualize from complex, multidimensional data. PCoA analysis was displayed by WGCNA package, stat package and ggplot2 package (Wickham, 2009) in R software (version 2.15.3). Multi-response permutation procedure (MRPP) was used to determine dissimilarities of microbial community structure between groups implemented in the R package vegan (version 2.5–6) (Mielke & Berry, 2001). Unweighted Pair-group Method with Arithmetic Means (UPGMA) Clustering was performed as a type of hierarchical clustering method to interpret the distance matrix using average linkage. The relative abundance of OTUs that most likely explain the differences between groups was evaluated by LEfSe (linear discriminant analysis (LDA) Effect Size) analysis (Segata et al., 2011).

Firmicutes/Bacteroidetes ratios

The non-parametric Wilcoxon rank–sum test was performed to compare Firmicutes/Bacteroidetes ratios between groups (L, OV, OB and T2DM) (p < 0.05). The comparisons were visualized as a boxplot by ggplot2 (Wickham, 2009).

Multivariate statistical analysis

In order to gain a deeper insight into the dietary consumption profile of individuals in the different groups, we used multiple factor analysis (MFA) in FactorMineR version 1.42 (R software, version 3.6.1) (Lê, Josse & Husson, 2008). MFA is beneficial for simplifying and structuring variables into groups on the components (PCA-based eigenvalues). The MFA of dietary consumption in different BMI and T2DM groups included 25 variables (frequency of intake) that belonged to 4 variable groups (protein (8 variables), carbohydrate (6 variables), fiber (5 variables), and beverage (6 variables)). For evaluating the association between dietary habits, blood profiles and fecal gut bacteria in BMI and T2DM groups, 37 variables were included in the MFA, which the analysis of multidimensional distance between subjects was based on 25 variables pattern of food consumption (frequency of intake), 2 variables described blood profiles (HDL cholesterol and fasting glucose levels), and 10 variables belonged to relative abundance of fecal gut microbiota at genus level. An integration of confounding factors (gender and age) with other concerned variables (blood profiles, dietary habits and fecal gut microbiota) also implemented in the MFA analysis. The age variable was categorized into six groups (20–29, 30–39, 40–49, 50–59, 60–69 and over 70) according to Png et al. (2016). Therefore, the variations were revealed by the influence of each variable on the principle components. Results of multivariate data analyses were extracted and visualized by Factoextra version 1.0.5. We used the mixOmics R package version 6.10.2 (Rohart et al., 2017) to determine associations between microbial communities, dietary consumption and blood profiles in the different groups. We applied sparse Partial Least Square (sPLS) analysis to explore relationships between these variables in each study group (L, OV, OB and T2DM). The sPLS “canonical mode” was used to specify the microbial OTU that most correlated with diets and/or blood profiles (Lê Cao et al., 2008, 2009). The high-dimensional data sets were visualized with clustered image maps (González et al., 2013).

Dietary consumption in different BMI groups and T2DM subjects

Dietary habits of the different groups were determined to enquire whether the observed data could support the gut microbial profile. Based on the frequency of consumption records, none of the subjects in any of the groups differed in the intake of dietary protein, carbohydrate, fiber and beverage with some exceptions: significant differences were noted regarding consumption of chicken (OB-T2DM comparison, q < 0.01), rice vermicelli (L-T2DM comparison, q < 0.05) and fermented fruits or vegetables (L-OV and L-T2DM comparisons, p < 0.05 for both).

The results of MFA revealed that individuals with T2DM displayed a notable variation on frequency of food consumption from the rest of the groups. This was discernible on the factor map indicating the first two dimensions accounting for 31.6% of variance (Fig. S1). The ellipse of L, OV and OB groups had a strong overlap compared to the ellipse representing T2DM group with 95 % confidence. These results indicate a lower variability of dietary consumption among the different BMI groups when compared to the T2DM group (coordinate = −1.42, p < 0.001). The frequency of fermented fruit/vegetable consumption significantly described the first dimension (r = 0.60, p < 0.001), and that variation in food consumption among groups supported the reduction of this particular type of dietary fiber intake in subjects having type 2 diabetes.

Blood profiles in BMI and T2DM groups

Information regarding subject status of the study is shown in Table 1. There were statistical differences between the average age, BMI, weight, HDL cholesterol, and fasting glucose in four study groups assessed by Tukey–Kramer post-hoc test (q < 0.05). HDL cholesterol was significantly lower in the OB groups compared to L (q < 0.001), OV and T2DM groups (q < 0.05), and was also significantly lower in OV vs. L groups (q < 0.05). Unsurprisingly, fasting glucose level increased with increasing BMI with highest level in T2DM group (q < 0.05).

Composition of prokaryotic fecal microbiota in BMI and T2DM groups

A total of 3,408,383 reads were obtained from 16S rRNA amplicon sequencing with an average of 113,613 reads per sample for a total of 30 samples. Using 97% identity criterion for determining OTUs, we obtained 504.33 ± 33.15 OTUs per sample (range 454–579 OTUs). Gut microbiota of all samples was classified into 995 OTUs, 145 genera, 82 families, 53 orders, 31 classes and 18 phyla. Shared and unique OTUs among different groups are shown in Venn diagram (Fig. 1). The total number of OTUs presented in the L, OV, OB and T2DM groups were 832, 811, 753, 852, respectively. The number of shared OTUs in all groups was 588; 729 OTUs (22.44%) were shared between L and OV groups, 752 OTUs (23.15%) between L and T2DM groups, 671 (20.65%) between L and OB groups, 730 (22.48%) between OV and T2DM groups, 657 (20.22%) between OV and OB groups, and 629 (19.37%) between T2DM and OB groups (Fig. 1A). Regarding the number of non-shared OTUs, non-diabetic subjects (merging OTUs of L, OV and OB groups) had four times more than diabetic subjects (T2DM). Specifically, 139 non-shared OTUs (17%) were associated with the non-diabetic group, whereas 33 OTUs (4%) were uniquely present in T2DM subjects (Fig. 1B). The rarefaction curves of microbial diversity estimators for thirty samples reached plateau phase, indicating that most microbial species had been detected in all samples (Figs. S2A−S2C).

Fecal microbiome community diversity (richness and evenness) in the four groups was characterized using ACE, Chao1, observe-species, Shannon, Simpson and Good’s coverage. Sequencing data and alpha diversity indices in each sample are presented in Tables S2 and S3. Significant differences of overall bacterial community structure across the four groups were found in the ACE, Chao1, and observe-species indices (Fig. 2). Specifically, microbial communities of L and T2DM groups had significantly greater species richness as compared to those in OB. No significant differences in the diversity of communities (species richness and evenness) were found across the four groups by Shannon and Simpson indices, suggesting a similar pattern of the community composition in all groups. Nevertheless, microbial communities of all four groups had high species-level diversity as indicated by Simpson index, the value of which approached 1. This implies that as species richness and evenness increased, diversity also increased.

The top ten phyla of microbial communities across the four groups were Bacteroidetes, Firmicutes, Proteobacteria, Fusobacteria, Actinobacteria, Verrucomicrobia, Cyanobacteria, Tenericutes, Elusimicrobia, and TM7. No significant differences were detected among groups (p < 0.05) (Figs. S3A and S3B). Top ten of bacterial genera with high relative abundance were used to construct phylogenetic relationships (Fig. S3C). Based on the similarity threshold, some bacterial species related to Prevotella genus were clustered in (Prevotella) as their discrete lineages distinct from other known species within this genus. The representative OTUs assembled in (Prevotella) consist of Prevotella tannerae (1 OTU), Uncultured bacterium (7 OTUs) and Prevotellamassilia timonensis (2 OTUs). The latter is a newly identified bacterial species in the human gut (Ndongo et al., 2016).

We further compared the median differences of the Firmicutes:Bacteroidetes ratios between groups with 95% confidence interval at phylum level (Fig. S3D). No statistically significant difference was noted among any of the groups. Among the top ten dominant genera, Prevotella and Bacteroides accounted for the largest proportion in all sample groups (Figs. 3A and 3B). Several significant differences among the dominant taxa were also found without the Benjamini–Hochberg method. Faecalibacterium showed significant differences in OV-OB (p = 0.044) and OB-T2DM (p = 0.011) comparisons (Fig. 3C), while (Prevotella) differed in L and OV comparison (p = 0.034) (Fig. 3D). Gut microbial alterations at species levels in four groups of samples (Figs. S4A and S4B) present significant differences of butyrate-producing bacteria (Faecalibacterium prausnitzii) in OV-OB (p = 0.044) and OB-T2DM (p = 0.011) comparisons (Fig. S4C), whereas there was a statistically significant difference between OV and T2DM group (p = 0.026) with respect to the relative abundance of Prevotella copri (Fig. S4D). Bacteroides coprophilus was also significantly different in OV-OB (p = 0.026) and OV-T2DM comparisons (p = 0.015). Analysis of statistical differences of microbial species abundance among groups by LEfSe showed marked differences between subject groups (Fig. 4). In the histogram, the colors represent taxa that were found to be more abundant in OB group compared to the other groups (both for the positive (Figs. 4A, 4C, and 4D) and the negative score (Fig. 4E)). In this regard, there were six taxa (Veillonellaceae, Prevotella, Dialister, Megamonas, Faecalibacterium, and Faecalibacterium prausnitzii), two genera (Clostridium and Prevotella) and one species (Bacteroides coprophilus) enriched in OB, T2DM and OV groups, respectively, while no enrichment of species was noted in the L group.

Beta diversity analysis of fecal microbiota in BMI and T2DM groups

Qualitative (unweighted UniFrac) and quantitative beta (quantitative measure) diversity measures yielded substantially different perspectives on the factors (BMI groups and/or the disease) that may be involved in structuring bacterial diversity. Unweighted UniFrac showed less distance of samples than Weighted UniFrac (Fig. S5B). PCoA based on Unweighted UniFrac revealed clearer patterns of microbial variation (Fig. 5B). Gut microbial communities in OB group were significantly different from those of the T2DM group (p < 0.05), whereas no dissimilarity was observed in comparison with other groups. Furthermore, more similar community composition was observed in OB and OV groups with some overlaps also identified by UPGMA (Fig. S5B). Clustering analysis suggested association of BMI and/or a disease with the variations in bacterial community compositions among subjects. Conversely, PCoA and weighted UniFrac did not clearly discriminate microbial communities among groups implying that there was no strong association with BMI and/or T2DM in this population (Fig. 5A). When taking the relative abundance of each type of OTUs into account, the results displayed similarities in bacterial composition, suggesting that indistinguishable communities may result from the number of organisms collected in the dominant phylum (Firmicutes, Bacteroidetes and Proteobacteria) (Fig. S5A).

Associations between dietary habits, blood profiles, and fecal gut microbiota in BMI and T2DM groups

MFA revealed the variables that mostly contributed in explaining the variations regarding dietary habits, blood profiles, and the relative abundance of gut microbiota of subjects in different BMI groups as well as the T2DM group. The factor map of the MFA generated by data integration of all variables showed the significance of Dim 1 and 2 that explained 14.3% and 12.8 of observed variability, respectively (Fig. 6). The distinct or similar profile of individuals, indicated by the ellipses on both axes of the MFA, mainly resulted from the variation in the blood profiles: HDL cholesterol levels were lower in the OB group than those observed in the other groups (coordinate = −1.20, p < 0.05) and higher fasting glucose levels were firmly correlated with T2DM (coordinate = −1.54, p < 0.001) in comparison with the BMI groups. In addition to the blood profiles, the genera of bacteria that associated with Dim 1 included Fusobacterium (r = 0.55, p < 0.01), Bacteroides (r = 0.52, p < 0.01), Prevotella (r = −0.54, p < 0.01), and Faecalibacterium (r = −0.47, p < 0.01). Fusobacterium was associated with Dim 2 (r = 0.47, p < 0.01). For dietary consumption, dairy products and mixed rice variables were negatively correlated with Dim 1 (r = −0.63, p < 0.001) and Dim 2 (r = −0.64, p < 0.001), respectively. The second dimension was described by a beef variable with a correlation coefficient of 0.53 (p < 0.01). The relationships of the relative abundance of gut bacteria and the frequency of food intake were nevertheless interpreted as moderate correlation to both dimensions (moderate variance). Concerning the contribution of all variables to describe the differences between individuals, the current analysis suggested that blood profiles seemed to have most influence on the variability of 30 subjects with different BMIs or with T2DM. Furthermore, the association between dietary intake, blood profiles, and microbial OTUs at the genus level from sPLS analysis are additionally summarized in Table 2.

Gender and age highly contributed to the variations of subjects in different BMI groups and T2DM group (Fig. S6). Integration of five variables (gender, age, blood profiles, dietary consumption, and the relative abundance of gut microbiota) by the MFA revealed that gender and age were the top two variables that highly explained individual variation in terms of blood profiles, dietary consumption, and the relative abundance of gut microbiota. The opposed pattern of T2DM females (coordinate = −3.12, p < 0.001) at age levels of four (50–59) and five (60–69) to Dim 1 was explained by the consumption of fish (r = −0.82, p < 0.001), whereas the consumption of fermented fruits or vegetables (r = 0.52, p < 0.01) and dairy products (r = 0.51, p < 0.01) as well as the abundance of Faecalibacterium (r = 0.49, p < 0.01) tended to be prevalent in OB females (coordinate = 1.92, p < 0.05). HDL cholesterol levels were negatively correlated to Dim 2 (−0.78), particularly among L females (coordinate = −3.0, p < 0.001) at the age level of three (40–49) (coordinate = −2.58, p < 0.001) that displayed higher HDL profile than OB and T2DM. In Dim 3, Escherichia was moderately correlated to the dimension (r = 0.51, p < 0.01) and was predominant in OB male at the age of 71 (coordinate = 3.25, p < 0.01). The distinct pattern of T2DM females marked in Dim 4 (coordinate = 1.52, p < 0.05) resulted from high fasting glucose levels with the correlation dimension of (r = 0.51, p < 0.01).

When all blood tests were taken into account, several important variables were maintained (Fig. S7), however, their variations were described by different dimensions as compared to the MFA in Fig. S6. Different profiles between OB females and T2DM females were illustrated in Dim 1, where the consumption of chicken (r = 0.70, p < 0.001), dairy products (r = 0.59, p < 0.001), and the abundance of Faecalibacterium (r = 0.53, p < 0.01) were less prevalent among T2DM females (coordinate = −2.99, p < 0.001) in comparison with OB females (coordinate = 2.21, p < 0.01). A distinct cluster of L females observed in Dim 2 (coordinate = −2.65, p < 0.001) was mostly described by HDL cholesterol levels (r = −0.72, p < 0.001) as being analogous to Dim 2 of Fig. S6. A moderate correlation of Escherichia dimension (r = 0.42, p < 0.05) also displayed in Dim 2 in parallel with OB male at the age of 71 (coordinate = 3.22, p < 0.01). Faecalibacterium (r = 0.44, p < 0.05) and FBS (r = 0.37, p < 0.05), however consistently presented in Dim 3, was predominant in OB females (coordinate = 1.90, p < 0.05). Interestingly, most OV individuals were negatively correlated with Dim 2 as resulted from LDL cholesterol levels (r = −0.45, p < 0.05) that tended to be higher in both OV male (coordinate = −2.92, p < 0.05) and females (coordinate = −1.51, p < 0.001), especially at the age level of two (30–39) (coordinate = −2.10, p < 0.001).

Exclusion of gender and age variables showed a strong overlap of non-diabetic (BMIs) and diabetic (T2DM) subjects, though no specific group was defined in relation to the variations of blood profiles and gut microbiota, indicating dispersal of variables across BMI and T2DM groups (Fig. S8). In Dim 1, the distribution of subjects on the MFA map was mainly influenced by the abundance of three major genera including Bacteroides (r = 0.68, p < 0.001), Prevotella (r = −0.65, p < 0.001), and Faecalibacterium (r = −0.65, p < 0.001). Cholesterol (r = 0.49, p < 0.01), triglyceride (r = 0.43, p < 0.05), and diastolic blood pressure (r = 0.43, p < 0.05) levels moderately correlated with the dimension. Furthermore, systolic blood pressure (r = 0.68, p < 0.001) and HDL cholesterol (r = −0.54, p < 0.01) levels mainly described the individual variance in the second dimension.