Phocea, Pseudoflavonifractor and Lactobacillus intestinalis: Three Potential Biomarkers of Gut Microbiota That Affect Progression and Complications of Obesity-Induced Type 2 Diabetes Mellitus

Introduction

Type 2 diabetes mellitus (T2DM) is the most common metabolic disease around the world, which results in ineffective use of insulin and has accounted for 90% of diabetes.¹ The International Diabetes Federation reported that there would be 693 million people over the world with diabetes by 2045.² The progression of diabetes can lead to serious complications such as vascular disease and even death.³

With the worldwide prevalence of obesity, it is becoming more common to these obesity-associated diseases of insulin resistance (IR).⁴ Except for many well-known contributing factors of T2DM such as obesity, sedentary lifestyle, and aging, now the dysfunction of the gut microbiota is gaining more and more attention for its association with energy harvest.^5,6 Turnbaugh et al found the obesity can raise phylum-level changes in gut microbiota, which reflected in lower bacterial diversity, gene expression changes and influences of metabolic in twins.⁷ Studies so far have already proved that the difference between T2DM patients and the normal group of gut microbiota plays a vital role in IR and metabolic syndrome-developing progress.^8,9 And this may be due to the contribution of gut microbiota to the metabolites and glucose metabolism by regulating the composition of bile-acid pools and glucose metabolism modulation signaling.¹⁰ As is well-known, bile-acid plays an important role in regulating the metabolism of lipids, which serve as major structural components of biological membranes and participate in energy storage and signal transduction.¹¹ Lipid subclasses such as phospholipids, sphingomyelins, and triglycerides, have all been related to IR and T2DM in humans.^12,13 Since lipidomics related mass spectrometry-based techniques can conduct the qualitative and quantitative analysis of the lipidome profiling from cells or biological fluids, tissues and the whole organism, it naturally becomes practical tools to unravel pathological mechanisms and identify biomarkers for metabolic diseases.¹⁴ We use the extracting method of nonpolar substances, which can count in most lipids and a few other metabolites of plasma samples. Since compared with serum, metabolomics studies based on mass spectrometry using plasma samples can perform complete metabolic profiles and accurately detect metabolic dysregulations.¹⁵

At present, though many studies have been done in the IR field concerning about metabolism mechanism with T2DM or profiling the function of gut microbiota,¹⁶ concurrently combining gut microbial profile and plasma metabolomics is supposed to be more intuitive to demonstrate about the connection between gut microbial and the host’s metabolic changes.

Zucker diabetic fatty (ZDF) rat is an animal model characterized by leptin receptor deficiency that induces hyperphagia leading to inchoate fat.¹⁷ Naturally occurring glucose tolerance makes ZDF rats an ideal experimental subject to simulate the progress of human T2DM.¹⁸ ZDF model (fa/fa) group was induced to develop diabetes by a high-energy diet, and ZDF control (fa/+) group had no diabetes symptoms, but both groups were obese. In the present study, we used gene sequencing technology, combined with the ultra-performance liquid chromatography–mass spectrometry technology, to explore differences of the intestinal microbial metagenomics and plasma metabolites between groups. We found 12 potential biomarkers of gut microbiota and some of them were highly correlated with diverse perturbed metabolites, mostly lipids.

A detailed understanding of the pathophysiology of T2DM is a prerequisite for the development of preventive strategies. In particular, the identification of early metabolic alterations is promising in the study of etiological pathways and may further help to identify high-risk individuals. By further detecting the intestinal microbial condition of T2DM, and its impact on metabolites in blood, we have the chance to use gut microbiota as a new breakthrough in the treatment of diabetes and metabolic syndrome.

Materials and Methods

Animal

Six 8-week-old male ZDF model (fa/fa) rats fed by Purina #5008 Lab Diet (crude protein 23.5%, crude fat 6.5%) for 3 weeks and 6 age-matched ZDF control (fa/+) rats fed by normal rodent diet (Beijing Keao Xieli Feed Co., Ltd.) were purchased from Charles River Laboratories (Beijing, China). All the animals were housed at standard room temperature, with a 12-h light/12-h dark cycle with free access to food and water. Then, both groups received standard 1 week’s adaptive feeding with a basic rodent diet. The animal protocols were approved by the Institutional Animal Care and Use Committee of Shandong University.

Random Blood Glucose Test

We used Random blood glucose (RBG) to examine diabetes in the fa/fa group. Glucose levels in tail blood samples were measured from weeks 8 to 9 using a glucometer (Sinocare Inc, China). Body weight and the result of RBG showed that we have successfully built the model of ZDF rats (Supplement Figure 1).^19,20

Sample Collection

Collection of both the feces and blood samples was taken at 12 weeks.²¹ The feces were stored in the sterile tube and then immediately snap-frozen at −80°C for microbial community analysis.²² And the blood taken from the abdominal aorta was collected under inhalation anesthesia with isoflurane then stored in the tube containing EDTAK₂. The plasma separated from collecting blood after centrifuging at 3000rpm for 20 mins and then frozen at −80°C.^22–25 Rats were killed by inhalational isoflurane for 60 seconds.

16S rRNA Sequencing and Analysis of Fecal Microbiota

Using the Sodium dodecyl Sulfate method, total genome DNA from samples was extracted. 16S rRNA genes of distinct regions were amplified, using specific primer 515F (5ʹ-GTGYCAGCMGCCGCGGTAA-3ʹ) and 806R (5ʹ-GGACTACHVGGGTWTCTAAT-3ʹ) with barcode. All PCR reactions were carried out in 30 µL reactions with 15 µL of Phusion^® High-Fidelity PCR Master Mix (New England Biolabs); 0.2 µM of forward and reverse primers, and 10 ng template DNA. Mix the same volume of 1×loading buffer (contained SYB green) with PCR products and operate electrophoresis on 2% agarose gel for detection. PCR products were mixed in equidensity ratios. Then, mixture PCR products were purified with Gene JET^TM Gel Extraction Kit (Thermo Scientific). Sequencing libraries were generated using Ion Plus Fragment Library Kit 48 rxns (Thermo Scientific) following the manufacturer’s recommendations. The library quality was based on the Qubit@ 2.0 Fluorometer (Thermo Scientific) and finally sequenced on an Ion S5^TM XL platform, ultimately generated 400 bp/600 bp single-end reads.

According to Cutadapt, we performed Quality filtering.²⁶ Then, we used UCHIME algorithm to compare reads with the reference database (Silva),²⁷ and Clean Reads was finally obtained. Sequence analysis was performed by Uparse software (Uparse v7.0.1001).²⁸ We defined the sequences with ≥97% similarity as the same operational taxonomic units (OTUs). We screened the representative sequence for each OTU based on Mothur algorithm (Silva Database). OTUs abundance information was normalized using a standard of sequence number corresponding to the sample with the least sequences on which subsequent analysis of alpha and beta diversity was all performed. Alpha diversity and beta diversity are calculated with QIIME (Version 1. 7. 0) and displayed with R software (Version 2.15.3). The difference of beta diversity index was achieved by using parametric test and non-parametric test between two groups, including T-test and Wilcox test. Principal component analysis (PCA) was applied to reduce the dimension of the original variables using R software. Metastat analysis performed permutation test among groups at various classification levels to obtain the p-value and then corrected the p value by Benjamini and Hochberg False Discovery Rate to obtain the q-value. The filter value of Linear discriminant analysis Effect Size (LEfSe) was set to 4 by default.

Plasma Metabolomics

Metabolites were prepared for LC-MS extracted with two rounds of 3000 rpm, 4°C centrifugation by methanol and dichloromethane. Chromatographic analysis was implemented with Waters ACQUITY ultra-performance liquid chromatography (UPLC) I-Class system. After separation by UPLC, mass spectrometry was performed by using a Xevo G2-S Q-TOF (Gentech™) with an electrospray ionization (ESI) source (Waters, Manchester, UK). In positive ion-mode, MS parameters were as follows: 2.5 kV capillary voltage; 24 V cone voltage; 100°C source temperature; 400°C desolvation temperature; 800 L/h desolvation gas flowing and 50L/h cone gas flowing. In negative ion-mode, MS parameters were as follows: 2.5 kV capillary; voltage; 25 V cone voltage; 100°C source temperature; 500°C desolvation temperature; 800 L/h desolvation gas flowing and 50 L/h cone gas flowing. We inputted the raw data into the Progenesis QI (Waters) for peak alignment, peak list containing the retention time, m/z, and peak area of each sample. Lipidmaps (http://www.lipidmaps.org), HMDB (http://www.hmdb.ca), NIST (https://chemdata.nist.gov) and an in-house lipid database were used for metabolite identification. The mass error used was 5 ppm. We completed the statistical with R software (Version 2.15.3). The threshold was set as VIP (Variable Importance in the Projection) >1.0, FC (Fold Change) >2.0 or FC<0.5, with p-value<0.05, so as to screen out the difference metabolites that met the setting conditions.

Correlation Analysis of Fecal Microbiota and Plasma Metabolites

P-value and rho of the correlation coefficients between the different metabolites and significant gut flora were calculated by Pearson statistical method (P-value was taken to 6 decimal places). P-value <0.05 (*) and |rho| ≥0.8 was considered to be statistically significant.

Results

OUTs Clustering and Diversity Analysis

A total of 925,442 clean reads with a mean length of 253 nt were obtained from all ZDF rats by Ion S5TMXL sequencing analysis after removing chimera. OUTs (Operational Taxonomic Units) were clustered based on clean data with 97% identity. Finally, 801 OTUs were found out and annotated from Class level to Species level. Figure 1A1 and A2 showed the composition of gut microbiota (top 10) at different levels based on relative proportions by each group. Alpha Diversity was used to analyze the richness and diversity within the community. The individual rarefaction curve tended to approach the saturation plateau (Figure 1B), implying high sampling coverage (beyond 99%) was accomplished in all samples. For beta diversity, Principal Component Analysis (PCA) showed that the OTUs composition within fa/fa group and fa/+ group was similar, respectively (Figure 1C).

Figure 1 (A) Relative Abundance. Relative gut microbiota abundance at class (A1) and family (A2) level in fa/fa group and fa/+ group. Sample Name gN.1 to gN.6 belonged to fa/+ group, gP.1 to gP.6 belonged to fa/fa group. (B) Rarefaction curve of each sample with 97% identity. (C) PCA Plot. The abscissa represents the first principal component (PC1). The ordinate represents the second principal component (PC2). The percentage represents the contribution value of PC1 or PC2 to the sample difference.

Taxonomic Composition of Gut Microbiota

The dominant phyla in both the fa/+ group and fa/fa group are Firmicutes (58.8% and 69.9%, respectively, mean) and Bacteroides (39.1% and 28.2%, respectively, mean) (Supplementary Table 1). Candidate bacteria belonging to the class Clostridia (25.8–42.5%), Bacilli (6.3–51.0%), and Erysipelotrichia (1.0–10.1%) were detected as major bacteria from within the phylum Firmicutes (Figure 1A1). Candidate bacteria belonging to the family Muribaculaceae (5.1–22.8%), Bacteroidaceae (1.0–18.3%), and Prevotellaceae (0.2–14.2%) were detected as major bacteria from within the phylum Bacteroidetes (Figure 1A2). At phylum level, the average percentage of Firmicutes and Proteobacteria was lower in fa/+ group (58.8%,0.27%) compared to fa/fa group (69.9% and 1.07%) while Bacteroidetes, Verrucomicrobia and Actinobacteria ratios were higher in the fa/+ group (39.1%, 0.95%, 0.146%, compared with 28.2%, 0.12%, 0.0996%)(Supplementary Table 1).

Intestinal Flora Biomarker of the fa/fa Group and fa/+ Group

Metastat analysis suggests that three genera of Phocea, Pseudoflavonifractor, and Flavonifractor from within the family Ruminococcaceae have significant differences between groups (Figure 2). LEfSe (LDA Effect Size) analysis showed statistically different biomarkers, including class Erysipelotrichia, order Erysipelotrichales, family Erysipelotrichaceae, Lactobacillaceae, Peptostreptococcaceae, genus Blautia, Turicibacter and species Lactobacillus intestinalis, Romboutsia (Figure 3A). Biomarkers only belonging to the fa/+ group include Peptostreptococcaceae, Romboutsia, Turicibacter (Supplementary Figure 2A). And Lactobacillus intestinalis and Blautia are only present in fa/fa group (Supplementary Figure 2B). Cladogram (Figure 3B) suggested that characteristic biomarkers from each group are homologous, respectively.

Figure 2 Metastat analysis uses the hypothesis test of species abundance data between groups to obtain the p-value, which was then corrected to obtain the q-value, then filtered species had significant differences according to it. *q < 0.05. Abbreviations: f, family; g, genus; s, species.

Figure 3 LEfSe (LDA Effect Size) can discover and interpret high-dimensional biomarkers, and emphasize statistical significance and biological correlation between groups. (A) Distribution histogram of LDA values shows the species only with LDA Score larger than 4, the length of which represents the impact size of different species (LDA Score). (B) Cladogram represents the taxonomic levels from phylum to species. The diameter of each small circle is proportional to the relative abundance. Abbreviations: c, class; o, order; f, family; g, genus; s, species.

Plasma Metabolites Profiling of ZDF Rats

The detection methods were used to extract substances with less polarity or non-polarity. Except for a small number of other metabolites, most of the substances identified were lipids. According to the results of plasma metabolites analysis, fa/+ and fa/fa groups could be differentiated readily using Principal Component Analysis (PCA), the first two components of which (i.e., PC1 and PC2) were finest separated, as performed by positive or negative EPI spectra in Figure 4A. Figure 4B illustrates that T2DM unbalanced the lipidic profiles of the gut microbiome, with 292 increased and 65 decreased metabolites. The structures of these metabolites were diverse, including lipids mostly, bile acids, fatty acids, amino acid derivatives and so on. Among these metabolites, as illustrated in Figure 5, the variation of several specific phosphatidylcholine (PC) and phosphatidylethanolamine (PE) is consistent with previous reports.^29,30 Deoxycholic acid (DCA), which is a secondary bile acid, downturned in T2DM.

Figure 4 (A) Controls were separated from model rats in metabolite profiles by PCA. (B) The plasma metabolic profile of model group changed with 357 molecular features being significantly changed compared with control group (p < 0.05).

Figure 5 Relative Quantification of (A) PC(16:0/16:1(9Z)) (B) PC(18:0/20:3(5Z,8Z,11Z)) (C) PE(16:0/22:6(4Z,7Z,10Z,13Z,16Z,19Z) deoxycholic acid (D) that differed significantly between fa/fa group and fa/+ group. *p<0.05; ***p<0.001. Abbreviations: PC, phosphatidylcholine; PE, phosphatidylethanolamine.

Correlation Between Gut Microbiome and Plasma Metabolites

Clear correlations could be identified between the biomarkers of gut microbiota and metabolism unbalance in plasma by the correlation matrix generated by counting the Pearson’s correlation coefficient (Figure 6). The microflora that affects most metabolites was genus Phocea, Pseudoflavonifractor and species Lactobacillus intestinalis, especially Phocea (Table 1). But there are still some metabolites that have a correlation with other flora. For example, DCA, which decreased 1.1-fold in fa/fa group, positively correlates with class Erysipelotrichia, including order Erysipelotrichales, family Erysipelotrichaceae, and genus Turicibacte, family peptostreptococcaceae and genus Romboutsia. In summary, T2DM induces a significant taxonomic perturbation in the intestinal flora, which substantially alters the metabolomic profile in blood.

Table 1 Results of Association Analysis of Differential Metabolites and Biomarkers of Bacterial Flora. The Results in the Table Only Contain the Metabolite with p Value Less Than 0.05 and Rho Was Greater Than 0.8, as Well as Their Associated Flora

Figure 6 Correlation plot showing the correlation between perturbed gut bacteria and altered plasma metabolites. Abbreviations: p, phylum; c, class; o, order; f, family; g, genus; s, species.

Discussion

Obesity is one of the main causes of T2DM and associates with IR.^31,32 Intestinal microorganisms can affect the energy homeostasis and metabolic level of the host through the transformation of food.³³ In ZDF rats, T2DM was induced by excessive consumption of a high-energy diet. Here, we took advantage of this rodent model to investigate the influence of intestinal microorganisms on the host in the process of obesity-induced T2DM.

Our study demonstrates that intestinal microbiome composition and plasma metabolites of ZDF rats in the fa/fa group displayed different patterns compared to the control group. Specifically, some of these differential metabolites are functional correlative with the perturbation of intestinal flora. As indicated in the results, not only the occurrence of T2DM would interfere with the growth of the intestinal microbial community, but might in turn further aggravate the disturbance of host metabolite homeostasis. These findings may provide novel insights into exploring the etiology and development mechanism of T2DM.

Our results first demonstrated significant alterations in proportion and abundance in the gut microbial community. Consistent with our findings, previous studies also have suggested the variations in the gut microbiota composition existed in the development of metabolic syndrome and obesity. For example, Hildebrandt et al reported that rapid changes, reducing Bacteroidetes and increasing Firmicutes, was observed in the high-fat-fed mice.³⁴ These changes were supposed to be associated with accelerating calorie absorption and gain weight by helping the host to digest polysaccharides or activate AMP-activated protein kinase (AMPK).^35,36 At phylum level, the abundance of Firmicutes is the most, Bacteroidetes is the second, and other phyla are few in both groups. At other levels, 12 gut flora biomarkers are obtained by statistical analysis of great significance. A previous study suggested that low abundance Erysipelotrichia can be used to predict the occurrence of T1D, with 69% prediction accuracy in NOD (no obesity diabetes) mice.³⁷ In this study, we found that class Erysipelotrichia and order Erysipelotrichales, family Erysipelotrichaceae, which belong to the offshoots of the same evolutionary line, all decreased in fa/fa group. Thus, these results supported the association with the low abundance of this evolutionary line and T2DM. At family level, besides Erysipelotrichaceae, we also noticed that Peptostreptococcaceae also decreased in fa/fa group, while Lachnospiraceae showed the opposite trend. Plenty of members of family Lachnospiraceae have been linked to obesity in human³⁸ and mice studies. They could accelerate the course of diabetes³⁹ and positively correlates with the degree of atherosclerotic lesions in apoE (apolipoprotein E) KO (knockout) mice,⁴⁰ which also suggested that high abundance of Lachnospiracea might potentially predict serious complications like cardiovascular disease. A population research suggested Peptostreptococcacea enriched in a normal group compared to T2DM group,⁴¹ which was consistent with our study. At genus level, Phocea, Pseudoflavonifractor, and Flavonifractor that all belonged to family Ruminococcaceae had higher abundance in fa/fa group. Many previous studies indicated that the variation trends of Ruminococcaceae and Lachnospiraceae, both of which belong to energy metabolic advantage families, were often the same, and their high abundance was related to metabolic disorders.^40,42,43 On the other hand, the relative abundance of Flavonifractor has a positive correlation with inflammation cytokines, like interleukin, cholesterol and LDL levels in plasma.^44,45 Flavonifractor may cleave the flavonoid C-ring and degenerate quercetin to increase oxidative stress and trigger inflammatory markers such as C5a, IL-6, IL-7, IL-8, IL-1β and IL-21.⁴⁶ Notably, by obstructing the insulin signal transduction, IL-6 can induce insulin resistance.⁴⁷ Moreover, IL-1β, IL-6, and IL-8 can lead to β-cell deterioration in pancreatic islets.⁴⁸ So an increase in Flavonifractor may be associated with the development of T2DM.

Some studies suggested that Phocea and Pseudoflavonifractor are probiotic,⁴⁹ but an increase of them in fa/fa group may be due to a decrease of else profitable strain. Other differential bacteria at genus level included Blauita, Turicibacter, and Romboutsia. The former increased, but the latter two declined. According to the previous reports, rich abundance Blauita is highly related to glucose disorders, low levels of insulin and arterial hypertension.^50,51 So in this case abundant Blautia and Lachnospiraceae were possibly related to cardiovascular and cerebrovascular complications. The abundance of Turicibacter significantly increases in the early stage of diabetes.²¹ The genus Romboutsia was isolated from the digestive samples of ileum for the first time in 2014.⁵² For the past few years, studies have shown that Romboutsia is negatively correlated with fasting glucose, insulin, and HDL-C, and positively correlated with serum TC, TG, and LDL-C.^53–56 A recent study found Romboutsia increased in T2DM mice model.⁵⁷ Our experiment firstly illustrated a clear correlation between Romboutsia and T2DM in a diabetes rat model. But we need further evidence to see if this strain is also comparably sensitive in human T2DM prediction. The augment of species Lactobacillus intestinalis is higher in fecal samples of diabetic patients than in those of control subjects.⁵⁸ Therefore, our study validated these biomarkers’ fundamental roles in predicting T2DM progression and associated complications.

Meanwhile, the plasma metabolite profile also shows variation between groups. PC(32:1) and PC(38:3) were increased in diabetes in a prospective investigation in European populations, which is consistent with our results.³⁰ In our analysis, PC(32:1) is accurate to PC(16:0/16:1(9Z)) and PC(38:3) is accurate to PC(18:0/20:3(5Z,8Z,11Z)). And a study found that some lipids, including PE(16:0/22:6), were potential biomarkers for T2DM patients.⁵⁹

Changes in the bile acid pool can cause a shift in some parameters of the metabolic syndrome (such as insulin resistance, hyperglycemia, hepatic steatosis, low HDL cholesterol levels, and cardiovascular risk).⁶⁰ DCA, a secondary bile acid, is positively correlated with GLP-1, which has a hypoglycemic effect.⁶¹ Although no study yet has confirmed that bile acid can be used as an early predictor and diagnostic marker for T2DM, the results of our study supported that DCA should be further studied. DCA was also proved to be one of the endogenous ligands of Farnesoid X Receptor (FXR) that improved various dysfunctions of glucose, cholesterol, lipid, triglyceride and cholesterol metabolism.^62–64 We already know that FXR-deficient mice have lower serum leptin concentrations, increased plasma-free fatty acid levels, and insulin resistance.^65,66 Mechanistically, the decrease of DCA is related to the reduced absolute abundance of Bacteroidetes, the bacteria involved in DCA generation.⁶⁷ Our research indicates reduced DCA also positively correlates with some members in Firmicutes, inclusive of the evolutionary branching of class Erysipelotrichia, family Peptostreptococcaceae and genus Turicibacter, Romboutsia. So it is reasonable to take DCA as a potential plasma metabolites biomarker in T2DM according to the correlation between gut microbiota abundance.

Except for DCA, the top 20 metabolites selected by p-value are most influenced by genus Phocea, Pseudoflavonifractor and species Lactobacillus intestinalis. The relative abundance of these three strains all decreased in the fa/fa group. Furthermore, genus Phocea, Pseudoflavonifractor both belong to family Ruminococcaceae which is related to energy metabolism and insulin sensitivity.⁶⁸ It may reflect that some members of family Ruminococcaceae disrupt the metabolic balance then aggravate the metabolic disorders of T2DM patients.

However, we cannot draw a definite conclusion on whether T2DM changes the intestinal microflora, which in turn accentuates further disease, or pathogenic factors, such as unhealthy diet and obesity, changes gut microbiota first, thus speed up the development of sub-health leading to T2DM, or both of these situations exist.

Our study makes the influential pathway of intestinal microbe on T2DM clearer, but we also need to study further about the causal relationship between intestinal microbe and T2DM, so as to find the target of T2DM treatment or preventive measures.

Conclusion

The results showed that compared with the control group, the intestinal microbiome and plasma metabolomics of ZDF fa/fa rats were changed. There are 12 potential biomarkers of gut microbes associated with T2DM complications and hundreds of different metabolites, including DCA, were identified. In addition, correlation analysis showed that some intestinal bacteria were highly correlated with changes in plasma metabolites. In summary, these data indicated that perturbation in the intestinal microbial community in T2DM hosts is related to disorders of metabolic level. Thus, these florae may be crucial biomarkers for future studies on the pathogenesis of T2DM and its complications.