Antibiotic-Induced Alterations in Gut Microbiota Are Associated with Changes in Glucose Metabolism in Healthy Mice.

The gut microbiome plays an important role in health and disease. Antibiotics are known to alter gut microbiota, yet their effects on glucose tolerance in lean, normoglycemic mice have not been widely investigated. In this study, we aimed to explore mechanisms by which treatment of lean mice with antibiotics (ampicillin, metronidazole, neomycin, vancomycin, or their cocktail) influences the microbiome and glucose metabolism. Specifically, we sought to: (i) study the effects on body weight, fasting glucose, glucose tolerance, and fasting insulin, (ii) examine the changes in expression of key genes of the bile acid and glucose metabolic pathways in the liver and ileum, (iii) identify the shifts in the cecal microbiota, and (iv) infer interactions between gene expression, microbiome, and the metabolic parameters. Treatment with individual or a cocktail of antibiotics reduced fasting glucose but did not affect body weight. Glucose tolerance changed upon treatment with cocktail, ampicillin, or vancomycin as indicated by reduced area under the curve of the glucose tolerance test. Antibiotic treatment changed gene expression in the ileum and liver, and shifted the alpha and beta diversities of gut microbiota. Network analyses revealed associations between Akkermansia muciniphila with fasting glucose and liver farsenoid X receptor (Fxr) in the top ranked host-microbial interactions, suggesting possible mechanisms by which this bacterium can mediate systemic changes in glucose metabolism. We observed Bacteroides uniformis to be positively and negatively correlated with hepatic Fxr and Glucose 6-phosphatase, respectively. Overall, our transkingdom network approach is a useful hypothesis generating strategy that offers insights into mechanisms by which antibiotics can regulate glucose tolerance in non-obese healthy animals. Experimental validation of our predicted microbe-phenotype interactions can help identify mechanisms by which antibiotics affect host phenotypes and gut microbiota.

Introduction

The human gastrointestinal tract contains a multitude of microbiota, including bacteria, viruses, and fungi (Utzschneider et al., 2016). Their genome, although variable between individuals (Human Microbiome Project Consortium, 2012), is capable of a diverse set of functions that may influence the host’s metabolic and immune systems (Tremaroli and Backhed, 2012; Greer et al., 2013; Sanz et al., 2015), including normal homeostasis (Utzschneider et al., 2016). Changes in the gut microbes have recently been associated with various diseases (Qin et al., 2012; Karlsson et al., 2013; Wu et al., 2015). For example, changes in Lactobacillus, Clostridium, Ruminococcus sp., E. coli, Bacteroides, Akkermansia muciniphila are observed in diabetic and obese patients (Qin et al., 2012; Karlsson et al., 2013; Murri et al., 2013; Chakraborti, 2015; Kasai et al., 2015; Sanz et al., 2015). These diverse results indicate a need for a better understanding of the mechanistic roles specific taxa play in the regulation of host metabolic functions.

Antibiotics add an interesting dynamic to the host-microbiome relationship. Although, antibiotics are well-known to cause short (Perez-Cobas et al., 2013; Pallav et al., 2014; Panda et al., 2014) and long-term (De La Cochetiere et al., 2005; Jernberg et al., 2007; Dethlefsen et al., 2008; Jakobsson et al., 2010; Dethlefsen and Relman, 2011; Raymond et al., 2016) alterations in the gut microbiome, there is a lack of consensus on their effects on glucose tolerance, body weight and other metabolic parameters (Francino, 2015; Mikkelsen et al., 2016). Moreover, effects of antibiotics in lean, normoglycemic mice as compared to mouse obesity models have not been widely investigated. An intervention study in healthy, glucose tolerant young human males treated with 4-days broad-spectrum antibiotics cocktail showed shifts in the cultivable gut microbiota but no changes in postprandial plasma glucose and serum insulin (Mikkelsen et al., 2016). Due to the use of a broad-acting antibiotic cocktail in a short course as well as the use of fecal samples for culture-based bacterial assessment, this study provides limited insight on a comprehensive picture of changes in intestinal microbes and on associations between individual antibiotics and specific intestinal microbes. Understanding antibiotic-microbiome interactions and their effects on glucose metabolism in healthy mammals is critical for identifying initial changes in microbiota that eventually may lead to diseases such as obesity and diabetes.

In this study, we aimed to understand the regulatory mechanisms by which individual antibiotics and their cocktail influence the cecal microbiome and host phenotypes in lean mice, namely, gene expression and metabolic parameters. By treating lean mice with different antibiotics we sought to: (i) study the effects on body weight, fasting glucose, glucose tolerance, and fasting insulin, (ii) examine the changes in expression of key genes of the bile acid and glucose metabolic pathways in the liver and ileum, (iii) identify the shifts in the cecal microbiota, and (iv) infer interactions between gene expression, microbiome, and the metabolic parameters. We repeated the entire experiment twice and performed meta-analyses to increase the confidence of our results.

Materials and Methods

Mice and Antibiotics Treatment

Eight weeks old adult male Swiss Webster mice were initially purchased from Taconic Biosciences (Germantown, MD, United States). Mice were housed at the Laboratory Animal Resource Center at Oregon State University for 3–5 days for acclimation under standard 12-h light cycle with free access to food (5001, Research Diets) and water. Experimental procedures were carried out in accordance with protocols approved by the Oregon State University Institutional Animal Care and Use Committee. Mice were given single, cocktail, or no antibiotics for 4 weeks to create a stable altered microbiome. Antibiotics were administered in autoclaved drinking water individually, or in a cocktail for 4 weeks in the following concentrations: ampicillin (1 gl-1), metronidazole (1 gl-1), neomycin trisulfate (1 gl-1), and vancomycin (0.5 gl-1). This time course is consistent with standard antibiotic administration used in multiple studies for altering microbiota (Rakoff-Nahoum et al., 2004; Morgun et al., 2015; Greer R.L. et al., 2016). Each group consisted of five mice per experiment, total 30 mice per experiment, except for four mice in the cocktail group from the second experiment. Water consumption was monitored over the 4 weeks treatment period and all groups showed consumption equivalent to control water.

Glucose Tolerance Testing

Mice were fasted for 6 h during the light phase with free access to water. A concentration of 2 mg kg-1 glucose (Sigma–Aldrich) was injected intraperitoneally. Blood glucose was measured at 0 (immediately before glucose injection), 15, 30, 60, and 120 min with a Freestyle Lite glucometer (Abbot Diabetes Care).

Serum Collection and Hormone Measures

Mice were fasted for 6 h during the light phase with free access to water. Serum was collected via submandibular bleed using BD microtainer serum separator tubes. Fasting insulin was measured by ultrasensitive ELISA (Crystal Chem) according to manufacturer’s protocols.

Bacterial DNA Extraction, 16S rRNA Gene Library Preparation and PCR

Unflushed cecal tissue and content was suspended in 1.4 ml ASL buffer (Qiagen) and homogenized with 2.8 mm ceramic beads followed by 0.5 mm glass beads using an OMNI Bead Ruptor (OMNI International). DNA was extracted from the entire resulting suspension using QIAamp DNA Stool Mini Kit (Qiagen) according to manufacturer’s protocol. DNA was quantified using Qubit broad range DNA assay (Life Technologies). The V4 region of 16s rRNA gene was amplified using universal primers (515f and 806r) (Caporaso et al., 2012). Individual samples were barcoded, pooled to construct the sequencing library, and then sequenced using an Illumina Miseq (Illumina, San Diego, CA, United States) to generate pair-ended 250 nt reads. Quantitative PCR was performed for A. muciniphila as described in Schneeberger et al. (2015) with DNA for standard curve isolated from the cultivated microbe.

RNA Preparation and Gene Expression Analysis

Liver and ileum (flushed out of content) were collected and snap frozen prior to RNA extraction. Liver was homogenized using OMNI Rotor-Stator Homogenizer in Trizol and RNA was extracted using Trizol/chloroform extraction followed by the RNeasy Mini kit (Qiagen). Ileum RNA was extracted using OMNI Bead Ruptor and 2.8 mm ceramic beads (OMNI International) in RLT buffer followed by Qiashredder and RNeasy kit using Qiacube (Qiagen) automated extraction according to manufacturer’s specifications. Total RNA was quantified using Nanodrop (Thermo Scientific). Complementary DNA was prepared using iScript reverse transcription kit (Bio-Rad) and qPCR was performed using QuantiFast SYBR mix (Qiagen) and StepOne Plus Real Time PCR system and software (Applied Biosystems). Primers used for qPCR are listed in Supplementary Table .

Statistical Analysis of Phenotypic Data

An outlier value per group per experiment was removed (if p-value < 5%) for each phenotype (metabolic parameters and genes) using the default Grubb’s test from R package outliers v0.14 (Komsta, 2011). The data was log2 transformed and differential phenotypes (antibiotics vs. control) were detected using limma (Ritchie et al., 2015) (Bioconductor 3.4, BiocInstaller 1.24.0, R 3.3.2) per experiment. A combined Fisher’s p-value was calculated for each phenotype from the p-values for the limma t-statistic from each experiment. A false discovery rate (FDR) was calculated on the combined p-values. Change in phenotype was considered statistically significant if the phenotype had same direction of (abx/control) fold change in both experiments, individual p-value < 20% in each experiment, Fisher’s combined p-value (Fisher, 1932) <5% and FDR < 10%. The dot plots for the phenotypes were generated using R package ggplot2 (Wickham, 2009) and the GTT curves were generated using GraphPad Prism software v7.03.

Analyses of 16S rRNA Gene Sequencing Data

The samples were demultiplexed and forward-end fastq files were analyzed using QIIME v. 1.9.1 (Caporaso et al., 2010). The default quality filter parameters from QIIME’s split_libraries_fastq.py were applied to retain high quality reads (Phred quality score > = 20 and minimum read length = 75% of 250 nucleotides). A closed reference OTU picking with 97% sequence similarity was performed using UCLUST (Edgar, 2010) and Greengenes reference database v13.8 (DeSantis et al., 2006; McDonald et al., 2012) to cluster 16S rRNA gene sequence reads into OTUs and assign taxonomy. The reference sequence of an OTU from the Greengenes database was used to obtain species level taxonomic assignment using Megablast (Altschul et al., 1997; Morgulis et al., 2008) (top hit using default parameters). A threshold of 99% cumulative abundance across all samples in an experiment was used to retain abundant microbes, thus removing OTUs with approximately <0.01% abundance across all samples in that experiment. The read counts were normalized using cumulative sum scaling (Paulson et al., 2013) followed by quantile normalization.

The normalized OTU tables were used for diversity and statistical analysis. Briefly, a sampling depth of 200,000 sequences per sample was used for rarefaction. The alpha diversity metrics were calculated on unrarefied and rarefied OTU tables (Supplementary Table ). The Shannon diversity index (from rarefied data) for samples with and without antibiotics treatment was compared with a non-parametric t-test. The difference was considered to be statistically significant if the direction of (abx/control) fold change in both experiments is the same, individual p-value < 2% in each experiment, Fisher’s combined p-value < 0.1% and FDR < 0.1%. Beta diversity was calculated using weighted UniFrac (Lozupone and Knight, 2005) and the distances were used for PCoA (Gower, 1998) and visualized using EMPeror (Vazquez-Baeza et al., 2013). The taxonomic summary bar plots were used to visualize abundance at the phylum and order levels.

The log2 transformed OTU tables were used for limma analysis. Meta-analysis was performed using the same criteria as applied for phenotypes to identify differentially abundant OTUs. A heatmap with row scaling was generated for each experiment using R packages ggfortify v0.2 (Horikoshi and Tang, 2016) and gplots v3.0.1 (Warnes et al., 2016). Hierarchical clustering was used to group OTUs (rows) based on similar abundance patterns across the groups in the first experiment and the same row order was used for the second experiment without row-wise clustering.

Network Reconstruction and Prioritizing Microbe-Phenotype Edges

Spearman rank correlations were calculated between all pairs of genes, microbes, and metabolic parameters across all samples or per-group in an experiment. A combined Fisher’s p-value was calculated for each pair from the p-values for the correlation from each experiment. A FDR was calculated on the combined p-values separately for the following correlations: (i) within genes, (ii) within metabolic parameters, (iii) between genes and metabolic parameters, and (iv) between OTUs and phenotypes (genes or metabolic parameters).

We retained edges that satisfy the following criteria: the sign of correlation coefficients in the two experiments should be consistent, individual p-value of correlation within each experiment is <20%, combined Fisher’s p-value of all experiments <5% and FDR cutoff of 10% for edges without a microbial node (i, ii, and iii), whereas 1% for edges containing at least one microbial node (iv).

Next, the transkingdom network was generated (Dong et al., 2015; Morgun et al., 2015; Greer R.L. et al., 2016; Rodrigues et al., 2017) by keeping the criteria-satisfying phenotypic (i, ii, and iii) and OTU-phenotype (iv) edges, where the OTU has >0.5% median abundance across the two experiments in at least one group.

Finally, an OTU-phenotype edge was retained if it showed consistent sign of per-group Spearman correlation coefficient between the two experiments, principles of causality compliancy (Yambartsev et al., 2016) [i.e., satisfied fold change relationship between the two partners in the appropriate (abx vs. control) comparison] in at least one group, and the same sign of correlation coefficient across different groups. To put this bipartite network in perspective of the phenotypic connections a phenotypic edge was included (only during visualization) if its strength of correlation was stronger than at least one OTU-phenotype edges connecting the phenotypes. Network topology statistics, namely degree and betweenness centrality (BC), were calculated using NetworkAnalyzer (Assenov et al., 2008) in Cytoscape v3.5 (Shannon et al., 2003). These edges were ranked using a score of maximum (per-group OTU abundance) × absolute [median (per-group correlation)] to prioritize OTUs and the phenotypes they potentially affect, where the per-group OTU abundance and correlation are medians across the two experiments. The top hit of BLAST for the Greengenes representative sequence for an OTU was used to obtain species level identification.

Data Availability

Raw reads of 16S rRNA gene sequencing have been deposited at NCBI under BioProject PRJNA394608, Biosamples of SAMN07356206 – SAMN07356264, Sequence Read Archive SRP112596.

Results

Lean, normoglycemic male mice were left untreated, or were treated with ampicillin, metronidazole, neomycin or vancomycin, or a cocktail containing all four antibiotics for 4 weeks to study the effects of antibiotic treatment on glucose tolerance, genes involved in glucose and bile acid metabolism, and the gut microbiota. Antibiotics resulted in different patterns of changes in the metabolic parameters, gene expression, and intestinal microbiome.

Antibiotics Improved Glucose Tolerance in Lean Mice

No metabolic parameter worsened following antibiotics treatment (Figure and Supplementary Figure ). We observed that treatment with individual or cocktail of antibiotics reduced fasting glucose, but did not change body weight. Glucose tolerance improved upon treatment with cocktail, ampicillin, or vancomycin as indicated by reduced AUC of the GTT. Treatment with all antibiotics, including metronidazole or neomycin reduced fasting glucose levels, however, the latter two did not cause changes in systemic glucose tolerance. Fasting insulin was reduced only when the mice were treated with vancomycin. Overall, glucose metabolism was regulated by antibiotic treatment.

Metabolic parameters and gene expression in antibiotic-treated and control animals. (A) Summary table; the red and green colors indicate increase and decrease, respectively, in antibiotic treated group compared to the control. (B) GTT curves for the antibiotics treated and control groups in the two experiments. (C–E) Metabolic parameters and gene (F–O) expression represented as means with standard error bars. The red and blue colors indicate experiments one and two, respectively. Asterisks indicate parameters that show statistically significant differences upon antibiotics treatment compared to untreated control mice [same direction of (abx/control) fold change in both experiments, individual p-value < 20% in each experiment, Fisher’s combined p-value < 5% and FDR < 10%]. Cn, Control; Coc or Cc, cocktail; Amp or A, ampicillin; Met or M, metronidazole; Neo or N, neomycin; Van or V, vancomycin.

Antibiotics Changed Expression of Genes Involved in Glucose and Bile Acid Metabolism

Tissue specific host gene expression is important in many metabolic processes (Thomas et al., 2008; Chiang, 2013) and regulated by gut microbiota (Larsson et al., 2012). These, along with the knowledge that intestinal glucose metabolism can control systemic glucose levels (Saeidi et al., 2013), led us to examine the expression of key glucose and bile acid metabolic genes in the liver and the ileum.

The majority of the tested genes in the ileum showed changes in expression due to antibiotic treatment (Figure ). Ileum Hk2 and G6pase transcripts showed decreased expressions after treatment with cocktail, ampicillin, or vancomycin. Ileum Pck1 and Tgr5 mRNA were increased after treatment with ampicillin or vancomycin, but showed no changes after with cocktail. Ileal Hk1 and Glut1 did not change gene expression after antibiotics, whereas, Fgf15, Fxr, and Shp showed antibiotic-specific patterns in expression.

Only three genes showed differential expression in the liver following antibiotic treatment (Figure ). Fxr and G6pase showed increased and decreased expression, respectively, in ampicillin or vancomycin treated mice. Pck1 showed lower expressions in ampicillin treated samples. Hk2 and Insr genes in the liver did not change following antibiotics treatment.

Despite some variability in tissue specific behavior of genes in response to antibiotics, the improved glucose tolerance upon antibiotic treatment suggests that relationships between gene expression and metabolic parameters are mostly preserved across all groups. Hence, we constructed a correlation network consisting of (differentially expressed) genes and (differentially abundant) metabolic parameters using all samples per experiment (Figure ). Genes from the ileum, including G6pase, Hk2, and Fxr were strongly connected with the GTT-AUC. The Fxr gene in the liver was positively correlated with the ileum Tgr5 but negatively correlated with ileum Fxr and with fasting glucose and GTT. Altogether, this network indicates opposite effects of intestinal and liver Fxr on glucose metabolism. Furthermore, it also suggests that increased glycolytic gene expression program in ileum is connected to worsening of systemic glucose metabolism.

A network consisting of metabolic parameters and gene expression from liver or ileum. An edge indicates the sign of spearman correlation coefficients across all samples in the two experiments are consistent, individual p-value of correlation within each experiment is <20%, Fisher’s combined p-value of all experiments <5% and FDR < 10%. Red and green colors indicate increased and decreased median fold change (abx/control) for nodes, respectively; where the color intensity corresponds to the level of fold change (e.g., dark color indicates fold change ratio is further away from 1); diamond and rectangle shapes indicate genes and metabolic parameters, respectively. Blue and orange colors indicate negative and positive correlated edges, respectively.

Antibiotics Caused Shifts in Microbial Communities

Microbiome composition is known to be affected by antibiotics (De La Cochetiere et al., 2005; Jernberg et al., 2007; Dethlefsen et al., 2008; Jakobsson et al., 2010; Dethlefsen and Relman, 2011; Perez-Cobas et al., 2013; Pallav et al., 2014; Panda et al., 2014; Raymond et al., 2016) and involved in metabolic processes (Larsson et al., 2012; Tremaroli and Backhed, 2012; Sanz et al., 2015; Utzschneider et al., 2016), so we hypothesized that gut microbes might play a mechanistic role in the effect of antibiotics (Morgun et al., 2015; Greer R. et al., 2016; Greer R.L. et al., 2016) on host glucose metabolism (Caesar et al., 2012; Greer R.L. et al., 2016). Sequencing the 16S rRNA gene of the cecal microbiome from the two experiments provided a total of 14,321,948 high quality reads with mean length of 248.50 bases and standard deviation of 9.42. A threshold of 99% cumulative abundance across all samples per experiment retained 734 and 677 OTUs in the two experiments (overlap of 561 OTUs) with 5,450,867 and 5,525,927 assigned sequences to the OTUs. The alpha diversity metrics on the normalized and rarefied OTU tables are provided in Supplementary Table . As expected, the cocktail of antibiotics reduced the diversity of the samples compared to untreated or individual antibiotics (Supplementary Figures , ). Shannon diversity comparisons showed that alpha diversity decreased when treated with cocktail, ampicillin, or vancomycin (Supplementary Figure ).

A PCoA analysis using the weighted UniFrac suggested that the overall community composition from vancomycin and ampicillin treatment was closer to that when treated with antibiotics cocktail (Supplementary Figure ). At the phylum level, Firmicutes decreased in cocktail, ampicillin, and metronidazole treated samples. Bacteroidetes decreased upon cocktail and vancomycin treatment but increased when treated with metronidazole (Figures , and Supplementary Figure ; Sheet in Supplementary Table ). The treatment with antibiotics showed similar patterns of change in the abundant bacteria at the order level, while less abundant bacteria showed antibiotic specific changes (Supplementary Figures , ). Vancomycin treatment increased Verrucomicrobiales in both experiments compared to control, however, increase in the second experiment was extremely high (fold change = 17,480) compared to the first (fold change = 358). Of note, this order was presented by single member (A. muciniphila). Thus, we also analyzed the abundance of this microbe via specific PCR and confirmed differences between two experiments in vancomycin treated groups (0.04 and 9794.4 ng DNA A. muciniphila/ g cecal content in the first and second experiments, respectively).

Taxonomic plots showing mean bacterial abundance across the different groups at the phylum level.

Heatmap of differentially abundant OTUs in the different groups. The represented OTUs passed the same statistical criteria as that used for phenotypes. The rows (OTUs) were clustered based on the abundance in the first experiment and the same row order was used for the second experiment. The rows were scaled, and orange and blue colors indicate increased and decreased abundance, respectively. The ids and taxa names for the rows (OTUs) are provided in Supplementary Table . The boxes indicate OTU clusters with the same taxonomy: Boxes 1 (113/117), 2 (43/46), and 3 (43/45) primarily contained Clostridiales (phylum: Firmicutes), while Box 4 (38/38) had Enterobacteriales (phylum: Proteobacteria).

Cocktail, ampicillin, and vancomycin treated samples showed similar patterns of change at the OTU level as compared to the microbiome of control samples (Figure ), which may be related to the fact that only these antibiotic treatments were able to change GTT-AUC (Figure ). Prevotella sp. (OTU_189721) was the most abundant OTU in control (median abundance across two groups (24%), neomycin (38%), and metronidazole (17.5%) treated samples. Enterobacteriaceae family (OTU_1111294) was the most abundant in cocktail (38%) and vancomycin (28%), and the third most abundant in ampicillin (14%) treated samples. Bacteroides uniformis (OTU_589071) was the most abundant upon ampicillin treatment (22.8%), while A. muciniphila was the second most abundant in vancomycin treated samples (17.5%) (Supplementary Table ).

Microbes Are Associated with Changed Phenotypes

Gut microbiota can control the expression of many genes in the small intestine (Larsson et al., 2012). Therefore, we asked whether the antibiotic-induced changes in the microbiome were potentially connected to the observed changes in gene expression. We constructed a transkingdom network using all groups, consisting of genes, metabolic parameters, and OTUs, to identify candidate interactions whereby microbes can mediate changes in systemic glucose tolerance and found 131 OTU-phenotype edges.

To focus on microbe-phenotype relationships that are not affected by type of antibiotics, we retained the 40 edges (Figure ) that maintained the same sign of correlation coefficient between the various groups of both experiments and consistent with potential causal relations (Dong et al., 2015; Morgun et al., 2015; Greer R.L. et al., 2016; Rodrigues et al., 2017) in at least one group of both experiments. Overall, this means that while a strength of OTU-phenotype interaction may be weak for a particular antibiotic group, this interaction may still be important in mediating effects of antibiotics on the host in general. The abundance of a microbe and its strength of correlation with a phenotype are expected to be crucial in mediating the effects, hence these 40 edges were ranked using a score that takes into account the maximum per-group OTU abundance and the median per-group correlation strength with a phenotype (Figure ; See formula in section “Materials and Methods”).

Transkingdom Network. Red and green colors indicate increased and decreased median fold change (abx/control) for nodes, respectively; diamond, rectangle, and circle shapes indicate genes, metabolic parameters, and OTUs, respectively. Blue and orange colors indicate negative and positive correlated edges, respectively. We indicate a phenotypic edge if its strength of correlation (in the phenotypic network) is stronger than at least one OTU-phenotype edges connecting the phenotypes.

Ranking the transkingdom edges. The score accounts for the bacterial abundance and its correlation strength with a phenotype. Blue and orange colors indicate negative and positive correlated edges, respectively. A broken axis (indicated by “Z”) allows for clear visualization of the other pairs.

NCBI BLAST on an OTU’s Greengenes reference sequence was used to obtain its (closest) species level identification. Interestingly, associations between A. muciniphila with fasting glucose and liver Fxr showed as the top interactions suggesting a possible mechanism through which this bacterium can mediate systemic changes in glucose metabolism. Proteus mirabilis was negatively correlated with GTT-AUC. Bacteroides uniformis was positively and negatively correlated with hepatic Fxr and G6pase, respectively. Importance of phenotypes in the network was also determined by degrees of connectedness (degree) and BC. GTT-AUC (degree = 14, BC score = 0.65) and liver Fxr (degree = 8, BC score = 0.60) were the highly connected metabolic parameter and gene, respectively, as well as the key nodes in the largest connected component of the network. Overall, it suggests that gut microbiota potentially influences the liver metabolic genes and systemic metabolic parameters and mediates the effects of antibiotics on host phenotypes.

Discussion

Germ-free Swiss Webster mice showed improved glucose metabolism (Caesar et al., 2012), suggesting that microbiota regulate metabolism in this strain. Furthermore, Swiss Webster mice are traditionally outbred, so more similar to human population. Therefore, Swiss Webster mice were selected for this study. While research has been done to study the effects of antibiotics on microbiota and glucose tolerance in diseased models (Francino, 2015), these effects in lean, non-diabetic or normoglycemic mice are not well studied. Such a study can provide meaningful insights into the host-microbial interactions and consequences of antibiotics in healthy population, and may allow the prediction of protective mechanisms and risk factors for development of diabetes.

To the best of our knowledge, our study is the first to show the ability of antibiotics to change glucose metabolism in healthy mice. The reduced fasting glucose and GTT-AUC in two experiments, especially in the ampicillin, vancomycin, and cocktail treated samples suggest that antibiotic treatment cause systemic improvements in glucose tolerance. Although our observations of reduced GTT-AUC contradicts with the absence of change observed in healthy humans (Mikkelsen et al., 2015), the cocktail ingredients and time course of antibiotic treatment (1 week vs. 4 weeks) of the two studies may be more critical factors contributing to this disagreement than differences between two species (i.e., mice and humans). Additionally, the unchanged insulin secretion by cocktail treatment in our study is in agreement with the study performed in humans (Mikkelsen et al., 2015). Noteworthy, one study did not observe any changes in fasting glucose and insulin in chow fed C57BL6 mice when treated with broad-spectrum antibiotics (ampicillin, metronidazole, and neomycin) (Pang et al., 2013). This disagreement might be due to differences in the mice strain, gut bacterial communities in different mouse facilities and the antibiotics used in the cocktail. For example, effect of vancomycin on glucose metabolism can be partially attributed (at least for our second experiment) by increased abundance of A. muciniphila which is missing in some mouse colonies. While the two studies (Pang et al., 2013; Mikkelsen et al., 2015) have some discrepancies with our observations, there are numerous supportive studies using germ-free (Caesar et al., 2012) or diet-induced obese (Cani et al., 2008, 2014; Membrez et al., 2008; Carvalho et al., 2012; Hwang et al., 2015; Fujisaka et al., 2016) mice that have shown improved glucose tolerance in the absence of microbiota and with antibiotics usage and the consequently modulated microbiota.

The expression of key genes from the glucose and bile acid metabolism pathways were measured, since bile acid signaling plays an important role in glucose homeostasis (Nguyen and Bouscarel, 2008; Trauner et al., 2010; Prawitt et al., 2011; Chiang, 2013; Nie et al., 2015; Trabelsi et al., 2016). We observed well-known (and therefore expected) relationships between the tissue-specific expression patterns of different genes themselves and with the systemic metabolic parameters in vancomycin and ampicillin. Low hepatic Fxr causes increased gluconeogenesis (Ma et al., 2006) and bile acid synthesis (Duran-Sandoval et al., 2004), while increased liver Fxr (Li and Guo, 2015) and intestinal Fgf15 (Holt et al., 2003) suppress hepatic bile acid synthesis (Kong et al., 2012) and regulate hepatic glucose metabolism (Potthoff et al., 2011). Also, increased liver Fxr represses G6pase (Yamagata et al., 2004; Ma et al., 2006; Zhang et al., 2006), Pck1 (De Fabiani et al., 2003; Yamagata et al., 2004; Ma et al., 2006), and like repressed ileum Fxr (Jiang et al., 2015b), improves glucose tolerance (Ma et al., 2006; Zhang et al., 2013), similar to our results. Along the same lines, mice treated with Fgf15 showed improved glucose metabolism (Zhou et al., 2017) and increased intestinal Fgf15 expression represses liver G6pase and Pck1, key enzymes for liver gluconeogenesis (Potthoff et al., 2011). In line with our observations in vancomycin treatment, Fxr agonist obeticholic acid (OCA; Intercept Pharmaceuticals, New York, NY, United States) increased mRNA levels of Fgf15 and Tgr5 in the ileum of C57BL/6J mice without increase in ileum Fxr (Pathak et al., 2017), and an increase in ileum Fgf15 expression decreased plasma glucose levels even with low insulin levels (Potthoff et al., 2011). Supporting our results from ampicillin treatment, a study showed that treating mice on high fat diet with antibiotics cocktail inhibited Fxr signaling in the ileum but not in the liver, and observed decreased expression of Shp in the ileum (Jiang et al., 2015a). Furthermore, Fxr and Shp mRNA in the ileum were also reduced in germ-free Swiss Webster mice on chow diet compared to the conventionally raised (untreated) group (Sayin et al., 2013). Overall, these studies along with ours support the idea that bile acids repress gluconeogenesis (Modica et al., 2010). However, it is interesting to see that while metronidazole, neomycin, and cocktail do not show the above changes in gene expression there is still improvement in fasting glucose and/or glucose tolerance, suggesting that microbiota might play an even bigger role in mediating the effects of antibiotics on phenotypes through additional mechanisms not explored here.

In fact, microbiota can change the expression of many genes in the ileum (Larsson et al., 2012). Their study observed down-regulated Hk2, G6pase, and Shp in the ileum of germ-free mice when compared to conventionally raised mice on chow diet, supporting our results from cocktail (ampicillin, or vancomycin) treatment. Also, some of the changes in gut microbiota to antibiotics that we observe in our data, e.g., the increased Verrucomicrobiales following vancomycin (Hansen et al., 2012) and Enterobacteriales increase upon treatment with ampicillin or vancomycin (Ubeda et al., 2010), are well documented.

While the effects of microbes on systemic glucose tolerance in lean subjects are rarely studied, their ability to influence glucose metabolism is well-recognized (De Vadder et al., 2016). A good example of well-established causal relations between specific bacteria and glucose metabolism is beneficial effect of A. muciniphila. For example, it was shown that A. muciniphila was able to delay the onset of diabetes in the vancomycin treated mice (Hansen et al., 2012). Furthermore, multiple studies demonstrated that this bacterium can improve glucose metabolism in animal models and in humans (Everard et al., 2013; Zhang et al., 2013; Joyce and Gahan, 2014; Shin et al., 2014; Anhe et al., 2015; Dao et al., 2016; Greer R.L. et al., 2016). It might not be surprising that the negative correlation between A. muciniphila and glucose levels was detected as one of the top ranked edges in our unbiased transkingdom network, thus, providing extra confidence for our results about less investigated bacteria inferred in our analyses.

Our predictions provide insights into host-microbial interactions. For instance, our result of Bacteroides uniformis being correlated with hepatic G6pase and Fxr might indicate potential mechanisms by which this bacterium improves glucose tolerance (Gauffin Cano et al., 2012). Also, it was shown that colonization with Bacteroides thetaiotaomicron makes mice leaner comparing controls despite similar levels of food consumption. The ability and preference of Bacteroides thetaiotaomicron and Bacteroides ovatus to utilize polysaccharide rich diet (McNulty et al., 2013) may explain these effects. However, our result of the negative correlation between abundance of these two bacteria and Pck1 (simple sugar forming gluconeogenic enzyme) may suggest the effect of these bacteria on liver gluconeogenesis. Similarly, P. mirabilis is predicted to have a negative interaction with GTT-AUC in our study, but shown to be positively correlated in rats with and without high fat diet (Lecomte et al., 2015). This disagreement may be explained by different physiological pathways dominating in the same bacterial species in different host that has been clearly shown for other bacteria (Oh et al., 2010). Overall, our study offers testable hypothesis regarding critical microbe-phenotype associations.

Conclusion

We show that antibiotics alter systemic glucose metabolism in lean mice. In addition to reporting changes in the microbiota, expression of key genes from the glucose and bile acid metabolism pathways, and concomitant systemic metabolic measures, we delineate potential mechanisms by which microbes mediate these effects. While there is a general understanding of the different players and mechanisms of microbiome-mediated regulation of the glycemic response (Tremaroli and Backhed, 2012; Devaraj et al., 2013; Cani et al., 2014; Hartstra et al., 2015; Janssen and Kersten, 2015; Parekh et al., 2015; Sanz et al., 2015; Boulange et al., 2016; Marchesi et al., 2016; Stenman et al., 2016; Suez et al., 2016), a lot remains to be understood, especially in terms of identifying the precise pathways operating in host-microbiome interactions. Overall, our data strongly suggests that antibiotics affect systemic glucose metabolism via shaping gut microbial communities and consequently regulating gene expression programs in intestine and liver. Yet, treatment of germfree mice with antibiotics as well as colonization of germ-free mice with antibiotic modified microbiota are required to fully support above statement. Also, it is doubtful that different antibiotics use the same mechanisms of gene expression and microbiota changes to affect systemic glucose tolerance, and the limited number of samples per group makes it difficult to obtain antibiotic-specific mechanisms. Furthermore, while the taxonomical assignments of 16S rRNA sequencing of current study present natural challenges, further studies employing shotgun metagenomics sequencing will allow to overcome this limitation. Finally, our experimental design followed by a data-driven, systems biology approach of network analysis offers consistent and statistically significant interactions that may be integral in mediating the host-microbiome communication. Furthermore, this approach is a useful hypothesis generating strategy and future experimentation can help investigate the distinct mechanisms in the different antibiotics and eventually lead to personalized medicine (Zmora et al., 2016).

Author Contributions

AM and NS conceived the original idea, designed and supervised the experiments, analyses, and writing. RG conceived the original idea, designed and performed the experiments, and supervised the writing. RR designed and performed the analyses, and drafted the manuscript. XD performed the analyses. JW, MG, and KD performed the experiments and analysis, respectively. All authors wrote the manuscript, read and approved the final draft submitted.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.