Combining STAT with high-fat diet increases body weight
We first sought to confirm and extend our prior studies of the effect of STAT on murine development [12, 13], in both males and females (Fig. 1). Analysis of the whole-life growth curves shows that STAT mice were heavier than controls from the very first weights obtained after weaning at week 4 (males only), with differences continuing to the end of the experiment (Fig. 2a–c). Both male and female STAT-exposed mice had increased body weight over time compared to controls, with the major differences occurring after HFD initiation at week 13 (Fig. 2a). After introduction of HFD at week 13, weight gain of STAT mice was greater than in controls (males, 20.0 ± 2.5 g vs. 13.1 ± 3.7 g; p <0.001; females, 13.7 ± 5.8 g vs. 5.1 ± 2.4 g; p <0.001), showing that the antibiotic exposure potentiated the effects of the HFD. At 32 weeks, both STAT males and females remained significantly larger than controls (Fig. 2c). These studies confirm our prior findings of enhanced growth of mice in the STAT model [12, 13], with acceleration of the growth differences in the presence of HFD.
STAT with a high-fat diet increases body fat
Beginning at weaning, body composition of all mice was measured by DEXA. Although STAT mice tended to have slightly higher lean mass (Fig. 2d, e), the significant weight differences observed largely reflected fat mass (Fig. 2f, g), which were enhanced by HFD in both sexes. Measurements of bone composition (mineral density, mineral content, and area) were not significantly different in relation to sex, treatment, or diet throughout the experiment (Additional file 1: Figure S1, Panels E, F, and G respectively). Taken together, these data indicate that STAT led to weight gain predominantly in fat mass, beginning early in life, exacerbated by HFD, with little or no effect on lean mass or on bone development, under the conditions studied.
STAT does not markedly perturb host energy balance
To determine whether STAT was altering food intake or energy harvest, 21-week-old mice were studied in metabolic cages. For individually housed control and STAT male and female mice, we measured food and water intake and waste production for 5 days. Food intake in STAT males was not different compared to controls, but STAT females consumed fewer total calories daily than control females (Additional file 1: Figure S1A). Fecal calorie content (per gram) measured using bomb calorimetry did not vary by sex or exposure group (Additional file 1: Figure S1B). Neither net calories (Additional file 1: Figure S1C; calories IN minus OUT), nor the proportion of calories retained (Additional file 1: Figure S1D; IN minus OUT/IN) was altered by STAT exposure. These data provide evidence that STAT-related adiposity did not result from either increased appetite or enhanced energy harvest.
STAT affects glucose and insulin homeostasis
Based on the increased weight and adiposity phenotypes, we hypothesized that STAT would increase the incidence and severity of metabolic diseases, including type 2 diabetes (T2DM) and NAFLD. To address this hypothesis, we conducted several studies in STAT and control mice in the weeks prior to sacrifice, including tests of glucose and insulin tolerance. While there was no significant difference in recovery of glucose levels in the STAT and control mice, both groups had markedly impaired glucose tolerance and incomplete recovery (Fig. 3a, b); the obese mice in this study were highly glucose intolerant, regardless of STAT exposure.
In insulin tolerance tests, there was significant insulin resistance in both STAT males and females compared to controls (Fig. 3c, d), in the earliest time period after the insulin provocation. Because many of the control animals experienced severe hypoglycemic shock and had to be withdrawn from the test prior to 120 min, we lacked sufficient power for assessment across the usual course of the ITT. At the relatively high level of insulin used, the STAT mice were less sensitive to hypoglycemia than were the controls, due to their relative insulin insensitivity (resistance).
To further quantify the metabolic impact of STAT, we calculated the HOMA-IR index [14]. By this index, based on fasting glucose and insulin values, STAT was found to significantly increase insulin resistance (Fig. 3e) in both males and females. Alternatively, using a pre-defined threshold for elevated HOMA-IR scores, STAT males had a significantly higher incidence of elevated HOMA-IR (Fig. 3f; p <0.05) compared to controls. Although not statistically significant, only STAT females, and not control females, had elevated HOMA-IR scores (Fig. 3f). These results, consistent with the ITT results, point to substantial alterations in glucose regulation in the STAT/HFD model.
STAT affects metabolic hormones and inflammatory markers
Based on the altered glucose homeostasis observed in STAT mice, we measured six other hormones and inflammatory markers involved in metabolism, which we hypothesized would be differentially affected by the STAT exposure. As expected from the IPGTT and IPITT, fasting serum insulin (p <0.05; Fig. 3g) and C-peptide (p <0.001; Fig. 3h) were significantly elevated in STAT males, although not significantly affected in females. Concordant with the increased adiposity, serum leptin was increased in both STAT males and females (p <0.001 in both; Fig. 3i). In contrast, serum ghrelin levels were significantly lower in STAT males and females compared to control mice (p <0.05 in both; Fig. 3j). Since metabolic and sex differences may be related to levels of the pro-inflammatory cytokines-TNFα and IL-6, respectively [42], we examined these in the context of the experiment. Circulating IL-6 was significantly elevated in females (p <0.05; Additional file 5: Figure S4A) but not in males, and circulating TNFα (Additional file 5: Figure S4B) was not significantly elevated in either sex. These data reflect the enhanced obesity in mice exposed to both STAT and HFD, and provide further definition of the sex differences observed.
STAT affects hepatic steatosis
Upon sacrifice, fatty infiltrates in the liver were visible in 13 of 37 mice (Fig. 4a; 10/18 in STAT, 3/19 in control; p = 0.017). Based on these observations, we performed microscopic examination of the liver, grading histology using the NAFLD Activity Score (NAS) [21] (Fig. 4b). Scores for all STAT males were above the diagnostic level for NAFLD with values significantly higher than for controls (p <0.001) (Fig. 4c). Although hepatic injury was not as advanced in STAT females, values were significantly higher than in controls (p <0.01) (Fig. 4c). Fibrosis (Fig. 4d), evaluated using the same scoring system, was more severe and extensive in STAT than controls (p <0.05) in both males and females. Significantly more STAT mice had scores diagnostic for NAFLD (score >5, with fibrosis) compared to controls (Fig. 4e; p <0.001 males; p <0.01 females). These findings indicate that the combination of STAT and HFD increased the incidence and severity of NAFLD-like histologic lesions compared to HFD alone.
STAT alters hepatic lipid storage and metabolism
Based on the hepatic histology at week 32, we quantified the hepatic lipid content. Total cholesterol was increased in STAT mice to a greater extent than in controls (p <0.05) (Fig. 4f) in females, but not males. In both sexes, STAT livers had nearly twice the triglyceride concentration of controls (p <0.001 for each comparison) (Fig. 4g). Free fatty acids were elevated in STAT compared to controls (p <0.05) (Fig. 4h) in males, but not females, representing another instance of sex differences in responses to STAT.
Next, we assessed expression of several genes relevant to steatosis. Cd36 and Vldlr expression were increased in STAT to a greater extent than in control (Additional file 5: Figure S4C; p <0.05 in both), consistent with the increased lipid infiltration of the liver. However, genes involved in fatty acid metabolism, lipid droplet formation, fatty acid oxidation, and related transcription factors did not differ. When the data were analyzed by outcome rather than treatment group, expression of the cholesterol efflux regulator Abca1 was lower in mice that had more insulin resistance (Additional file 5: Figure S4D; p <0.05). Conversely, Cd36 was borderline elevated in mice that had increased insulin resistance (Additional file 5: Figure S4E; p = 0.055) and in those that had met criteria for NAFLD (Additional file 5: Figure S4F; p = 0.055).
STAT effects on phylogenetic diversity of the intestinal microbiota
To assess the relationship of the phenotypic changes to gut microbial composition, we first addressed parameters of community ecology, beginning with markers of community richness. Although α-diversity values were generally higher for the STAT mice compared to controls early in the experiment, the only significant difference (p <0.05, Mann–Whitney U test) was at week 8 (Fig. 5a).
Microbiota community structures are distinct between groups
To assess the microbial community structure determined by unweighted UniFrac analysis of the studied samples, we visualized selected time points in relation to dietary transitions by principal coordinates analysis (PCoA) (Fig. 5c). The weeks shown represent the last sample before weaning onto normal chow (week 4), prior to the transition from normal chow to HFD (week 11), shortly after the transition (week 16), and study end (week 30), respectively.
Three female mice received STAT but did not show changes in any phenotypic changes specific to the STAT exposure (see Fig. 2c; these mice were termed female non-responders (FnR)). Based on UniFrac distances, at week 4, community structure of two of the FnR mice were STAT-like, while the third was control-like (p >0.05), but by week 11, all three FnR communities were indistinguishable from those in other STAT-exposed female mice, continuing through week 30 (p <0.05, compared with controls at each week; data not shown). These findings suggest that microbiota differences linked to differential outcomes in the FnR mice may have occurred prior to week 11.
When samples were grouped by treatment, the UniFrac distances were significantly different between intra- and inter-group measurements at each week (Additional file 6: Figure S5; p <0.005 for all weeks), indicating that the community structures of the control and STAT groups were distinct. When specimens were grouped by clinical phenotype (NAFLD, insulin resistance) or by not showing the pre-defined disease definitions (healthy), there were distinct differences at weeks 4, 11, and 30 (Additional file 6: Figure S5; p <0.005), but not at week 16. These results provide evidence that before HFD initiation, the intestinal microbial communities in mice that developed disease were distinct from those that did not. Although the addition of HFD diminished this distinction, the communities again were separate, long after the transition (week 30).
Differentiating taxa
On the day of weaning (week 4), control mice were enriched in Firmicutes and Candidatus Arthromitus (“Savagella”) (SFB), while STAT mice were enriched in Bifidobacterium, S24-7, and Prevotella, as determined by LEfSe [24]. While some individual taxa differed, that pattern was unchanged while the mice were receiving normal chow (week 11). When the mice were receiving HFD (week 16), the patterns continued similarly, except that SFB and Prevotella differences disappeared, and Allobaculum and Actinobacteria, enrichment was seen in control and STAT mice, respectively (Fig. 6a).
In controls, the week 4 to 11 transition showed a shift from Firmicutes-dominance, whereas the STAT transition was accompanied by a bloom in Proteobacteria (Fig. 6b). The transition between weeks 11 and 16 differed from the earlier transition, with the selective power of HFD having similar effects on control and STAT mice. The further transitions between weeks 16 and 30 similarly affected control and STAT communities (Fig. 6b). Thus, the effects of HFD on individual taxa appear to overwhelm the continuing effects of STAT.
Based on the LEfSe results, we sought to determine if any taxa could predict whether a host would develop metabolic disease (defined as insulin resistance or NAFLD). To accomplish this, a Random Forest classification model was built to predict disease outcome (class) based on bacterial OTU relative abundances (features) for each week of life. Of particular interest was the observation that for early weeks (prior to week 6); six OTUs were predictive (Prevotella, Lactobacillus, Erysipelotrichaceae, SFB, and two different S24-7 OTUs). The model had substantially (more than twofold) better than random predictive power at nearly all time points (Fig. 6c).
To understand the developmental differences in microbial communities, we calculated microbiota-by-age z-scores (MAZ) [11, 25] to compare the communities observed in control and STAT that did not develop disease, with the STAT mice that did (Fig. 6d). Intestinal microbiota follow reproducible patterns of community succession during early life, allowing “microbiota age” to be used as a benchmark of normal intestinal development, as described in studies of humans [25]. In this model, a maturity difference from control indicates either accelerated or delayed development of an age-appropriate microbial community. At week 4, samples from the STAT mice that would later develop insulin resistance or NAFLD had significantly lower MAZ scores than controls (p <0.001), but differences were lost at weeks 11, 16, and 30. These data provide evidence that STAT can delay the normal development of the early-life microbiome, and that this delay is associated with elevated risk for metabolic diseases in later life.
Associations between host phenotypes and microbial taxa
We applied multi-level, sparse PLS models to fecal microbiota data to assess linear relationships between OTUs and seven host phenotypes (Fat, Lean, BMC, DMI, Weight, Weight + 1, and NAFLD). We verified the efficacy of a multi-level linear model by visualizing the within-subject portion of the clr-transformed data. ISOMDS indicates clear separation between subjects of differing groups (Fig. 7a compared to Additional file 4: Figure S3A). We also computed biplots for the sPLS model (Fig. 7b and Additional file 7: Figure S6B), with sample scores colored by (scaled and centered) response variable and significant OTUs, represented by a loading vector colored by phylum.
Overall, 29 taxa (about 4 % of the total) were selected by the fully specified sPLS model, and three additional OTUs (two Clostridiaceae, and Odoribacter) were found to not be significant at α = 10–2. However, we found a large number of significant associations between taxa and body composition phenotypes (Additional file 8: Table S2). With the exception of two S24-7 families, all other Bacteroidetes OTU abundance levels were positively associated with body mass phenotypes, while Firmicutes associations were mixed.
More specifically, we found that Lactobacillus (n = 2) OTUs to be significantly associated with Lean, BMC, DMI, and Weight and one other Lactobacillus directly associated with Fat, DMI, Weight, and Weight + 1. This is consistent with prior findings that Lactobacillus reuteri reduces abdominal fat and age-associated weight gain [43]. Turicibacter genera (n = 2) were found to be negatively associated with DMI and Fat, but were not significantly associated with other body composition measurements, consistent with prior studies of low-dose antibiotic exposure [12] and HFD feeding [44] in mice. A single Anaeroplasma genus was negatively associated with BMI, but not NAFLD, which is consistent with HFD administration in C57BL/6 J mice [45] and abundance enrichment in low-weight rabbits [46].
Notably, we found a negative association between an Allobaculum OTU and NAFLD, accompanied by significant positive associations to other body composition measurements. Our findings are consistent with the previous observations that Allobaculum has been directly correlated with adiposity after switch to a HFD [12], yet negatively correlated with the development of the metabolic syndrome and total cholesterol levels [47, 48]. Finally, we also find Candidatus Arthromitus (SFB, n = 3 OTUs) to be negatively associated with body composition phenotypes, primarily Weight, BMC and Lean (consistent with elevated levels of SFB in control vs. STAT mice reported in [12]) with one particular SFB OTU predicted to have additional associations with NAFLD, Fat, and Weight + 1.
Microbial network topology corresponds to host physiology
We next sought to develop a network model that would permit insights about microbial relationships with the physiology of the hosts studied. The PLS model that we used transforms the primary microbiota population data into a subspace that maximally co-varies with the host responses. Using a Gaussian mixture model with bootstrap stability validation of cluster assignment, we performed unsupervised clustering of these transformed data. These studies revealed that samples are best grouped into six clusters, each of which has a distinct phenotype profile (Fig. 8a). Clusters 1, 3, and 5 were primarily associated with STAT mice. The switch from normal chow to high fat diet largely corresponds to the transition from Cluster 3 to Cluster 5. Clusters 2 and 4 were associated with Control mice receiving normal chow or HFD, respectively. As such, Cluster 4 comprises the microbiota in fecal samples primarily from 18–30-week-old mice.
To identify whether changes in host physiology are also reflected in the global rewiring of the gut microbial community structure that we observed, next we inferred microbial association networks from each of the six sample groups and analyzed their global topological network properties. Using graphlet correlation distance as a global distance measure between networks, and using isometric MDS as an analytic tool, we inferred a low-dimensional embedding of the microbial association networks (Fig. 8b). Importantly, these largely recapitulate the transitions seen in the subspace clustering described above (Fig. 8a). Networks 2 and 3 are closest to network 1. Networks 3 and 5, representing the gut microbiome community in samples from mice that received STAT are distinct from networks 2 and 4, which represent the microbial communities in samples from control mice Network 6, which is inferred from samples of older mice, is distinct from all the other networks.
Since clusters are dominated by samples that were obtained from mice under specific experimental perturbations, we classified the networks as being dominated by STAT (clusters 1, 2, and 3) or Control (2 and 4) or by normal chow (NC) (clusters 1, 2, and 3) or by HFD (4, 5, and 6). Then we calculated several graph topology statistics to assess trends as a function of sample type (Additional file 9: Figure S7). Overall, NC and STAT networks comprise more taxa, have larger network diameters, and show lower average betweenness and degree centrality. These are ecological terms indicating a node’s centrality in a network and the number of cross-ties, and low values are consistent with greater dispersion within the network. The HFD and Control networks tended to be more modular. Finally, the NC and Control networks had higher assortativity at the phylum level; this means that under normal conditions in the absence of antibiotics or an abnormal diet, particular OTUs are more likely to be directly associated with common phyla than under antibiotic and HFD exposure.
We also analyzed OTUs that could potentially serve as keystone species in the different association networks. For each network, we identified the top two taxa that serve both as hub species (having high node degree) and as bottleneck species (as characterized by the highest betweenness centrality) (Additional file 10: Figure S8). Across all networks, these potential keystone taxa are largely represented by Lactobacillus, Lachnospiraceae, and S24-7 families. For instance, in network 1, the top two taxa are Eubacterium dolichum and Lactobacillus reuteri. While these OTUs are not predicted to be directly associated with host physiological changes, both species are known to have strains that are resistant to penicillin [49, 50] and require sugar and amino acid import for survival in the host GI tract. In particular, L. reuteri has been shown to be a key mediator in host and microbe interactions for processing carbohydrate metabolites [51].
In addition to changes in microbial compositions, we also analyzed whether overall network robustness correlates with host health, since microbial ecological networks should have evolved to be resilient to disturbances. One example of this concept would be redundancy in network wiring that may ensure access to a vital metabolite. Thus, we hypothesized that Western-style interventions would promote network fragility by disrupting a critical threshold of keystone taxa or by changing the flux of normal metabolic exchange.
To test this hypothesis, we used natural connectivity as a general stability metric of the inferred networks after simulated network “attacks”. We found that only the network from cluster 2 – control mice receiving normal chow – was reliably robust, independent of the specific node removal strategy (Fig. 8c). Network 4, representing the microbial community after the switch to HFD, showed a decrease in network robustness, yet remained more stable than most STAT networks. Interestingly, the natural connectivity of network 4 decreased at a slower rate when bottleneck taxa were removed compared to hub taxa. This property suggests an increased redundancy of bottleneck taxa in the absence of antibiotic exposure. Importantly, microbial networks inferred in the communities in the STAT-exposed mice were found to be particularly fragile under targeted attacks, independent of the diet.