Liraglutide-induced structural modulation of the gut microbiota in patients with type 2 diabetes mellitus

Study population

The 40 (including both male and female) T2DM patients were aged 25–83 years, with mean % HemoglobinA1c (HbA1c) concentration of 9.17% (76.67 mmol/mol) (5.3%) (34.43 mmol/mol) − 12.2% (109.84 mmol/mol), and body-mass index (BMI) scores >25 kg/m². HbA1c is the gold standard of glycemic control index and the standardization of HbA1c was considered to be important, which was proceeded by the National Glycohemoglobin Standardization Program (NGSP) and the International Federation of Clinical Chemistry and Laboratory Medicine (IFCC) (Weykamp, 2013). All participants had been treated with metformin as a monotherapy (stable doses for 2 months). Relevant exclusion criteria included: taking other glucose-lowering medications and those who had been taking antibiotics in the last 3 months. The recruited participants were randomly assigned to receive 1.2 mg of liraglutide for 4 months. BaselineHbA1c was assessed at before treated-with liraglutide. HbA1c change (liraglutide treatment response) was calculated based on treatment HbA1c value (after 4 months treated-with liraglutide) minus baseline HbA1c value. All 40 T2DM patients were switched from oral metformin to subcutaneous injections of liraglutide for 4 months once daily. The metformin-liraglutide substitution design was applied because a pre-study feasibility assessment reported that it was impossible to recruit sufficient subjects with liraglutide as their primary medication for T2DM (Hodge et al., 2016). All procedures were performed in accordance with the ethical standards of the Clinical Research Ethics Committee of the Nanyang Second General Hospital, and written informed consents were obtained from all participants included in the study.

Fecal DNA extraction and 16s rRNA gene sequencing

Total 80 fresh feces samples were collected in sterile collection tubes (Fisher Scientific, Waltham, MA, USA) before and after the 4-month liraglutide treatment and then stored at −20 °C until being transported to the research center in chilled styrofoam containers. They were subsequently stored on-site at −80 °C for further analysis. Microbial DNA was extracted from 200 mg fecal sample using the QIAamp PowerFecal Pro DNA Kit (QIAGEN) which contains a bead-beating step. Briefly, adding 200 mg fecal sample and 800 µl lysis buffer in a beads-containing tube and vortex at maximum speed for 10 min. Take 350 µL supernatant after 1 min centrifugation at 15,000 g and use it in the subsequent steps according to the kit instructions. DNA is finally eluted in 100 µl elution buffer for downstream applications. DNA was amplified using primers targeting the V4 region of the 16S rRNA gene (515F 5′-GTGYCAGCMGCCGCGGTA-3′, 806R 5′-GGACTACNVGGGTWTCTAAT-3′). PCR was run in a Veriti™ 96-Well Thermal Cycler PCR system (Thermo Fisher Scientific) using the following program: 95 °C for 3min, followed by 21 cycles of 95 °C for 30 s, 56 °C for 30 s, 72 °C for 30 s, with a final extension at 72 °C for 5 min. Mixed amplicons were pooled and the sequencing was conducted at Shanghai Biotecan Pharmaceuticals Co., Ltd. (Shanghai, China) using an Illumina Novaseq 6000 Sequencing system (Illumina, USA) according to the manufacturer’s instructions.

Sequence data analysis

Sequences were assigned to operational taxonomic units (OTUs) with 97% similarity (Greengenes database: http://greengenes.lbl.gov) into mothur (v.1.39.5), OTU taxonomy was assigned to the Greengenes database for the comparisons, using the Quantitative Insights into Microbial Ecology (QIIME) software package (Caporaso et al., 2010). Both Kyoto Encyclopedia of Genes and Genomes (KEGG) (Langille et al., 2013) and Clusters of Orthologous Groups of proteins (COG) (Liu et al., 2019) pathways were categorized using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PiCRUSt) and both were imported into STAMP (v.2.1.3) for visualization. To identify taxa with differing relative abundances between the two groups, linear discriminant analysis (LDA) effect size (LEfSe) analyses were performed on the website (http://huttenhower.sph.harvard.edu/galaxy). The cut off value was the absolute LDA score (log10) >3.0 with a p < 0.05. The alpha diversity “summary.single” script was used to calculate ACE, Chao1, Shannon and Simpson indexs in the mothur software package. “Beta_diversity.py” script was used to calculate the beta diversity in the QIIME software, and the Principal Coordinate Analysis (PCoA) was used to measure beta diversity in each group. Coabundance network analysis was based on Spearman correlation. Only genera present in at least 60% of samples were used for correlation analysis, and only connections with a r > 0.6 or < −0.6 and P < 0.05 were used for network building on the basis of Igraph R package and imported into Cytoscape (v.3.6.0) for visualization.

Statistical analysis

The data are shown as the mean ± standard error of the mean (SEM). Bacterial taxa that were differentially abundant were identified by using the Wilcoxon rank-sum test (for two groups) or Kruskal-Wallis test (for more than two groups) in R studio (v.3.6.1). Clinical data analyses between the two groups were conducted by Wilcoxon signed–rank test or t-test using Prism version 6.0 (GraphPad, San Diego, CA, USA). Linear regression analysis was performed using R studio (v.3.6.1). A p < 0.05 was considered statistically significant.

Baseline HbA1c is a key predictor of liraglutide treatment response

We recruited 40 individuals with T2DM to receive 1.2 mg liraglutide for 4 months. Clinical characteristics of the 40 individuals before and after the 4-month liraglutide treatment (termed L0 and L4, respectively) are presented in Table S1. As baseline HbA1c is a major predictor in blood glucose therapies, we assessed the relationship between baseline HbA1c and HbA1c change using linear regression analysis. The result indicated that baseline HbA1c explains 52.7% of liraglutide treatment response variation (R² = 0.527, β = −0.726, P < 0.0001, Fig. 1). In addition, BMI, HbA1c, homeostasis model assessment of insulin resistance (HOMA-IR), fasting blood glucose, 2-hour postprandial blood glucose, total cholesterol, triglycerides, HDL-C, and LDL-C were significantly lower in the post-liraglutide-treatment group than the pre-liraglutide-treatment groupusing Wilcoxon signed–rank test or t-test . However, there were no notable differences observed in fasting insulin, serum creatinine and the blood urea nitrogen between the two groups (Table S2).

Furthermore, due to many clinical characteristics may be associated with baseline HbA1c, adjusted for baseline HbA1c was suggested to avoid erroneous treatment response. So we used multiple linear regression to analyze the association between the baseline HbA1c and the baseline covariates (11 clinical characteristics), and the association between the baseline covariates and liraglutide treatment response after with or without adjustment for baseline HbA1c. However, the results indicated that none clinical characteristics were associated with both baseline HbA1c and liraglutide treatment response after without adjustment for baseline HbA1c. Only the blood urea nitrogen was associated with liraglutide treatment response after adjustment for baseline HbA1c (Table S3).

Characteristics of the pyrosequencing results

To characterize the effects of liraglutide on the gut microbiome, we performed high-throughput sequencing of the V4 regions of 16S rRNA genes from the 80 fecal samples collected from the two groups. A total of 8,822,088 high-quality sequences and 7488 OTUs (97% similarity) were obtained from the 80 samples, and the average number of reads was 100, 626 and OTUs was 296 for each individual (Table S4). The reads/OTUs were assigned to 34 different phyla, and the top 3 dominant bacterial phyla of the two groups were Firmicutes, Bacteroidetes and Proteobacteria.

Liraglutide modifies gut microbiota alpha and beta diversity in the T2DM

Since the participants have different duration of T2DM, we divided the 40 T2DM patients into three diabetes duration groups (short-duration group (<5 years); medium- duration group (5-10 years) and long-duration group (≥10 years)) and analyzed alpha and beta diversity in gut microbiota of the three groups. Alpha diversity is an index reflecting the variety of microbial species in stool samples and we also examined the patterns of beta diversity by calculating the dissimilarity in the composition of gut microbial (Rapacciuolo et al., 2019). The results indicated that the three duration groups had no significant difference in gut microbiota alpha and beta diversity (Figs. S1 and S2). Therefore, being treated with metformin for 2 months do not alter the fecal microbiome in T2DM patients. In order to evaluate the effect of liraglutide treatment on gut microbiota diversity, we analyzed the 80 fecal samples collected from pre- and post-liraglutide-treatment groups. The alpha diversity indices (the ACE, Chao1, Shannon and Simpson index) were used to describe alpha diversity. The community richness (ACE and Chao1) was significantly different between the two groups (Figs. 2A and 2B), showing a lower abundance in the post-liraglutide-treated group. However, there were no significant difference in the community diversity (Shannon and Simpson) of the two groups (Figs. S3A and S3B). As to the beta diversity, a principal-component analysis (PCoA) of a nonmetric multidimensional scaling plot (on a Bray-Curtis distance matrix) and unweighted UniFrac distances revealed significant qualitative differences in gut microbial community structure between the two groups (Figs. 2C and 2D), but there was no difference in weighted UniFrac distances (Fig. S3C).

Liraglutide alters gut microbiota composition

To assess the effects of liraglutide on gut microbiota composition, a taxonomy-based comparison was performed and a broad overview of our taxonomic data from the 80 samples was given in Fig. 3A. A heat map was constructed to visualize the top 30 genus in the two groups (Fig. 3A). We also conducted a coabundance network analysis to investigate how different gut bacteria could interact with each other. The results indicated that 4 months of liraglutide treatment promoted an increased number of positive correlations among microbial genera, especially those within Firmicutes and Bacteroidetes. Negative correlation was identified in Ruminococcus (Firmicutes) and Actinomyces (Actinobacteria) (Fig. 3B).

Phylogenetic and taxonomic profiles of gut microbiota

To analyze the statistical differences in microbial communities between the pre- and post-liraglutide-treatment of T2DM, the significant differences based on genomic characteristics were further confirmed by LEfSe analysis, which used LDA coupled with effect size measurement to identify bacterial taxa with sequences that differed in the abundance of the two groups.

Cladograms were obtained by the LEfSe analysis, which showed the most significantly difference at taxonomic levels between the two groups. The size of each circle represents the abundance of certain bacteria (Fig. 4A).

Upon analyzing the fecal samples, 36 phylotypes were identified as high-dimensional biomarkers (Table S5). Twenty-one species of bacteria were abundant in the pre-liraglutide-treatment group and 15 species were abundant in the post-liraglutide-treatment group. At the phylum level, the abundance of Fusobacteria was significantly higher in the pre-liraglutide-treatment group, while Verrucomicrobia and Actinobacteria were significantly higher in the post-liraglutide-treatment group. At the genus level, nine genera were found in significantly high abundances in the pre-liraglutide-treatment group. These genera included Acinetobacter (p = 0.016), Oscillospira (p = 0.013), Acidaminococcus (p = 0.021), Succinatimonas (p = 0.042), S24_7 (p = 0.008), Megamonas (p = 0.0005), Alistipes (p = 0.035), Fusobacterium (p = 0.017) and Megasphaera (p = 0.0002). The genera Collinsella (p = 0.011), Akkermansia (p = 0.002) and Clostridium (p = 0.002) were enriched in the post-liraglutide-treatment group (Fig. 4B).

Correlation between clinical features and gut microbiota

According to the LDA analysis (LDA value > 3.0), we discovered the relative abundance of 12 bacterial genera (Table S5) between the two groups. To investigate the inter associations between gut microbiota and clinical status of the host, we identified the statistical correlations between the 12 genera and the 11 clinical features, including age, sex, BMI, family history of diabetes, diabetes duration, family history of cardiovascular disease, family history of diabetic nephropathy (DN), family history of diabetic retinopathy (DR), family history of diabetic peripheral neuropath (DPN), smoking status and alcoholism status (Table S6).

First, we found that the abundance of Megamonas was significantly lower in the older age group (≥50) (p = 0.014), medium- and long-duration groups (p = 0.006) (Figs. 5A and 5B). What’s more, the levels of Clostridum genus was lower in the family history of diabetes (p = 0.043) group (Fig. 5C). Oscillospira was significantly different between the Non-DR group and DR groups, in which severe non-proliferative DR (NPDR) sub-group had the highest abundance of Oscillospira (p = 0.038) (Fig. 5D). We also found that the abundance of Megamonas was lower in the DR group (p = 0.005) (Fig. 5E). Last but not least, the higher abundance of Oscillospira was found in the DPN group, but not in the non-DPN group (p = 0.040) (Fig. 5F).

Effect of liraglutide-treatment administration on metabolic pathways

KEGG and COG pathway comparisons were performed to explore potential differences in the functional composition of the microbiome between the pre- and post-liraglutide-treatment group. Figure 6 highlights the significant differences in the distribution of metabolic pathways, and compares the microbiota functions between the two groups. As listed in Fig. 6A, seven differences were observed in KEGG pathways. Three KEGG pathways (including G protein-coupled receptors, Glycolysis/Gluconeogenesis and Selenocompound metabolism) were enriched in the post-liraglutide-treatment group, and 4 KEGG pathways (including tropane, piperidine and pyridine alkaloid biosynthesis, phenylalanine metabolism, porphyrin and chlorophyll metabolism, and styrene degradation) were enriched in the pre-liraglutide-treatment group.

As shown in Fig. 6B, 17 different functional COG pathways between the two groups were found. Six of them were enriched in the post-liraglutide-treatment group, including predicted ABC-type transport system involved in lysophospholipase L1 biosynthesis, acetate kinase, phosphatidylglycerophosphate synthase, predicted choline kinase involved in lipopolysaccharides (LPS) biosynthesis, predicted glutamine amidotransferase and putative cell wall-binding domain. In contrast, the pre-liraglutide-treatment group were enriched in ABC-type Co2+ transport system, periplasmic component, ABC-type Fe3+ transport system, periplasmic component, transcriptional regulator containing PAS, AAA-type ATPase, and DNA-binding domains, acetyl-CoA hydrolase, N-acyl-D-aspartate/D-glutamate deacylase, Na+/proline symporter, permeases of the drug/metabolite transporter (DMT) superfamily, precorrin-2 methylase, predicted metal-binding protein, pyruvate carboxylase and s-adenosylhomocysteine hydrolase.