Integrative analysis of gut microbiota composition, host colonic gene expression and intraluminal metabolites in aging C57BL/6J mice

Physiological differences between mice during aging

Body weight as well as food intake increased during aging (Figure S1A-B). Also liver weight and relative liver weight (% of body weight) were significantly higher in older mice (Figure S1C-D). To determine if metabolic health altered during aging, fasting insulin levels in plasma were assessed. These levels were significantly increased at 24 and 28 months compared to 6 months (p<0.01) (Figure S1E). Death rates of the mice remained low till the age of 20 months and started to increase afterwards (Figure S1F). Most mice that died prematurely suffered from multiple conditions (Data set S1).

Alterations in faecal and colonic microbiota composition during aging

Fresh faeces samples were collected longitudinally at 4, 8, 12, 18, 24 and 28 months of age in a subset of 9 mice. Sequencing of the V3-V4 region of the 16S rRNA gene revealed that the alpha-diversity, in terms of Shannon entropy, did not change during aging (Figure 1A). To explore variance in faecal microbiota composition, a Principal Component Analysis (PCA) on relative abundance of genera was performed. The PCA revealed that the first principal component, explaining 42.7% of the variation, separated the mice aged 4 and 8 months from the other age groups (Figure S2A). Dissimilarities in faecal microbiota composition between age groups were further investigated at two taxonomic levels. At the phylum level, the relative abundance of Proteobacteria was less abundant in young mice (4 months) and increased from 8 months onwards in most mice (Figure 1B). Verrucomicrobia and Actinobacteria were highly abundant at 4 months, but both disappeared in most mice between 4 and 12 months of life (<1%). By exception, the relative abundance of Verrucomicrobia in mouse 9 at 28 months of age was high compared to the other mice at the same age (Figure 1C).

Figure 1. Alterations in gut microbiota composition and metabolites during aging. (A) The alpha-diversity (Shannon entropy) measured by the Shannon entropy in fresh faeces. Boxes extend from the 25th to 75th percentile, line in the middle represents median, and whiskers represent minimum and maximum values. (B) Relative abundance (%) at phylum level in fresh faeces collected longitudinally in a subset of 9 mice. (C) The alpha-diversity (Shannon entropy) measured by the Shannon entropy in colonic luminal content. (D) Principal Coordinate Analysis (PCoA) based on Bray-Curtis distances showing dissimilarities between all individual samples (FF and CLC). (E) The average relative abundance (%) at phylum level in colonic luminal content. Legend corresponds to graph with relative abundance at phylum level in fresh faeces. (F) Redundancy Analysis (RDA) displaying the part of the variation in microbiota composition explained by the age groups and source. Sample plot showing the clustering of the samples. Legend corresponds to PCoA plot. The distance between the ages and sample sources approximates the average dissimilarity of the microbiota composition. The FF analyses were based on a subset of n=9 mice at consecutive ages (4, 8, 12, 18, 24 and 28 months). For the CLC analyses, n=10 mice per age group (6, 12, 24 and 28 months) were included. (G) Levels of fumarate, (H) an unspecified bile acid and (I) methanol detected by 1H-NMR in CLC at 6 months (n=6), 24 months (n=8) and 28 months (n=8) of age. Error bars represent standard error of the mean (S.E.M.). (J) Strongest Spearman correlation coefficients (-0.60.6) for the concentrations of the identified metabolites and the relative abundances of the 50 most abundant genera in colonic luminal content.

At genus level, the higher relative abundance of the Proteobacteria phylum in older mice was largely accounted for by a significantly higher abundance of the pathobionts Desulfovibrio spp. and Bilophila spp. (Table 1). However, the strongest aging-related change within this phylum was the significantly decrease in relative abundance of Parasutterella spp. during aging. The lower relative abundance of the Verrucomicrobia phylum was solely explained by Akkermansia spp. that significantly decreased during aging (Table 1). The relative abundance of Bifidobacterium spp. and Coriobacteriaceae spp. both decreased strongly and significantly explaining the age-dependent decrease in relative abundance of the Actinobacteria phylum. Finally, although the phyla Bacteroidetes and Firmicutes in itself did not significantly change during aging, several genera belonging to these phyla showed strong and significant changes between young and old mice (Table 1 and Table S1).

Next to the analysis of the fresh faeces samples, the microbiota composition of the colonic luminal content of 40 mice sacrificed at the age of 6, 12, 24 and 28 months of age was investigated. Similar to the results obtained from the fresh faeces samples, no significant changes in Shannon entropy were found (Figure 1C) and the PCA revealed that the youngest and oldest age groups explained most of the variation (Figure S2B). These findings were strengthened by the distance-based Principal Coordinate Analysis (PCoA), since the first principal coordinate, explaining 27% of the variation, also separated young and old mice (Figure 1D). Comparison of the relative abundance at phylum and genus level between fresh faeces and colonic luminal content revealed mostly similar trends, although Actinobacteria phylum was hardly abundant in colonic luminal content (Figure 1E, Table 1, Figure S3). To investigate the impact of both age and source (fresh faeces and colonic luminal content) on the total variation, a Redundancy analysis (RDA) was performed. Ordination axis 1 explained 34.0% of the variation and separated the youngest age groups of both sources. Interestingly, ordination axis 2 explained 8.7% of the variation and separated the colonic luminal content and fresh faeces samples, indicating that the explanatory variable age contributed more to the variation in microbiota composition than the source variable (Figure 1F). Taken together, these results reveal pronounced changes in gut microbiota composition during aging, both in fresh faeces and more moderate in colonic luminal content. Additionally, subtle differences between the faecal and colonic microbiota composition were found.

Metabolite profiles in colonic luminal content

To explore whether observed changes in gut microbiota composition may come with potential functional consequences, ¹H-NMR was performed on colonic luminal content at 6, 24 and 28 months of age. Concentrations of the SCFAs butyrate, propionate and succinate did not differ between the age groups, except for acetate that was significantly lower at 24 months compared to 6 months (p<0.05) (data not shown). Intraluminal concentrations of fumarate decreased during aging, as well as an unspecified bile acid (Figure 1G-H). It was not possible to exactly identify this bile acid, since this peak in the NMR spectrum (0.17 ppm) is located in a rather unexplored region. However, it can be ruled out that this is one of the bile acids as described by Kok and colleagues [21], since these bile acids were identified in a different region of the spectrum (0.6-0.8 ppm). The concentration of methanol increased in older mice (Figure 1I). A correlation analysis was performed to investigate whether the levels of these metabolites could be related to the relative abundances of the 50 most abundant genera. The strongest negative correlations were found for fumarate that was negatively correlated with the relative abundance of Desulfovibrio spp. (r=-0.651), an uncultured Ruminococcaceae spp. (r=-0.708) and an uncultured Coriobacteriaceae spp. (r=-0.681), and positively correlated with Parasutterella spp. (r=0.637) (Figure 1J).

Changes in colonic gene expression in 6, 12, 24 and 28-month-old mice

To explore the effect of aging on gene expression in the colonic wall, microarray analysis was performed on mRNA isolated from colonic scrapings from the same mice for which we also had data on the microbiota composition. A PCA carried out on the top 1,000 most variable genes revealed a clear segregation of the 6-months-old mice (Figure 2A). No clear separate clustering was observed for the mice from the other age groups, indicating that only gene expression at 6 months differed substantially from all other age groups. The highest number of significantly differentially expressed genes was found when comparing the youngest (6 months) versus the three older age groups (Figure 2B). There was a large overlap of genes (n=610) that had an increased expression when comparing 6 months with 12, 24 or 28 months (Figure 2C). Together, these results indicate that the expression levels of a substantial number of genes increased between 6 and 12 months, and remained increased also at higher ages.

Figure 2. Changes in colonic gene expression in 6, 12, 24 and 28-month-old mice. (A) Principal Component Analysis (PCA) plot using the top 1,000 most variable genes showing the first 2 principal components. (B) The number of differentially expressed genes for all age comparisons (IBMT q<0.01). (C) Venn diagram showing the number of overlapping and unique differentially expressed genes between 6vs12, 6vs24 and 6vs28 months. (D) Each significantly differentially expressed gene (6vs12, 6vs24, 6vs28 months, q<0.01) was assigned to the model profile that most closely matched their gene expression profile, executed by STEM. The model profiles with the highest number of genes assigned to (as indicated by the number in the box) are shown. (E) Top 10 most significant canonical pathways between 6 and 12 months as determined by Ingenuity Pathway Analysis (R.A. = Rheumatoid Arthritis). Numbers behind bars represent the number of differentially expressed genes (IBMT q<0.01) and numbers behind the pathways represent the total number of genes. (F) The gene expression (Log2) of the epithelial marker E-cadherin (Cdh1), mesenchymal markers OB-cadherin (Cdh11) and Vimentin (Vim), transcription factors Snail Family Transcriptional Repressor 1 and 2 (Snai1, 2) and Zinc Finger E-Box Binding Homeobox 1 and 2 (Zeb1, 2), Transforming growth factor beta 1 (Tgfb1) and Transcription Factor 21 (Tcf21). (G) Length of the colon (cm) at sacrifice. (H) Depth of colonic crypts (μm) measured in H&E stained colon samples. (I) Representative pictures of Haematoxylin and Eosin (H&E) staining and Sirius Red/Fast Green staining of colon tissue at 6, 12, 24 and 28 months (200x magnification). Error bars reflect standard error of mean (S.E.M.). *p<0.05.

Dominant cluster of genes similarly expressed between 6 and 12 months

Visualization of the expression profiles of the 1,990 genes that were differentially expressed in the colonic wall of 6-months-old mice and the three older age groups (Data set S2) revealed a clustering of the majority of genes in four different expression profiles (Figure 2D). The profile in red contained by far the largest number of genes (817 genes) and showed the most prominent expression change between 6 and 12 months, remaining mostly stable afterwards (Figure 2D). The green profiles (602 genes in total) and the purple profile (217 genes) represent genes with decreased and increased expression levels during aging, respectively. Gene Ontology (GO) enrichment analysis of the green and purple profiles did not reveal any significant GO categories. However, analysis of the red profile revealed numerous significant GO categories that were all related to the extracellular matrix (ECM) organization and other structural development processes (Table 2). Ingenuity Pathway Analysis (IPA) and Gene Set Enrichment Analysis (GSEA) carried out on the genes differentially expressed between 6 and 12 months confirmed that similar pathways were affected (Figure 2E and Table S2A). IPA also revealed the ‘Regulation of the Epithelial-to-Mesenchymal Transition Pathway’ as a significantly regulated pathway in the colon between 6 and 12 months of age (Figure 2E and Table S2C). To analyze this aspect in more detail we assessed the expression regulation of genes involved in the epithelial-to-mesenchymal transition (EMT) in our microarray data set (Table S2D). Remarkably, while the expression of the epithelial markers E-cadherin (Cdh1), Occludin (Ocln) and most Claudins was minimally altered, the expression levels of the mesenchymal markers OB-cadherin (Cdh11), Vimentin (Vim), Fibronectin (Fn1) and Vitronectin (Vtn) (Figure 2F and Table S2D) were strongly increased between 6 and 12 months and remained constant afterwards. In addition, the same expression pattern was found for transcription factors associated with EMT, i.e. Snail Family Transcriptional Repressor 1 and 2 (Snai1, 2), Zinc Finger E-Box Binding Homeobox 1 and 2 (Zeb1, 2) (Figure 2F), as well as other EMT-inducers, such as Transforming growth factor beta 1 (Tgfb1) and Transcription Factor 21 (Tcf21). IPA identified the transcription factor Transforming growth factor beta 1 (Tgfb1) as the strongest activated upstream regulator associated with differential gene expression between 6 and 12 months (Table S2E).

Change in colon morphology during aging

To study the potential morphological consequences of the observed changes in ECM- and EMT-related genes we measured the colon length at sacrifice as well as colonic crypt depths in H&E-stained Swiss rolls. The obtained results revealed an increase in colon length during aging (Figure 2G). A marked deepening of the crypts between 6 and 12 months of age was observed, which minimally increased afterwards (Figure 2H). Furthermore, Sirius Red/Fast Green staining of the same Swiss rolls showed that the ECM protein collagen gradually increased with age (Figure 2I).

Integrative analysis of microbiota composition and gene expression in the early phase of aging

To obtain insight into the correlations between differential gene expression during the early phase of aging and the colonic microbiota composition, an integrative analysis was performed based on data of the individual mice. The 50 genera with a relative abundance of more than 0.1% in at least one sample and the 817 genes identified by STEM (Figure 2D, red cluster) showed several strong (r<-0.80 or r>0.80) positive and negative correlations (top 25 strongest correlations shown in Figure 3). Bifidobacterium spp. was negatively correlated with 228 genes, including a cluster of EMT- and ECM-associated genes, e.g. Cdh11, Fn1, Zeb1, Zeb2, Vim,Vtn and Fgfr1, amongst others (Figure 3 and Data set S3). Additionally, Turicibacter spp. and an uncultured member of the Lachnospiraceae family were positively correlated with a large number of genes, including Platelet Derived Growth Factor Receptor Beta (Pdgfrb), Wnt Family Member 2B (Wnt2b) and Gremlin 2 (Grem2), amongst others. Together, these results might imply the potential role of gut microbiota in the regulation of ECM (re)organization and EMT, or vice versa.

Figure 3. Integrative analysis of colonic microbiota composition and gene expression in the early phase of aging. The 817 up-regulated genes as determined by STEM were correlated with the 50 genera with a relative abundance ≥0.1% in at least 1 sample, including the mice aged 6 and 12 months. Relevance network showing the genera (white circles) and genes (grey rectangles) that were correlated (r<-0.80 or r>0.80). Red lines represent positive correlations, blue lines represent negative correlations. Regarding the negative correlations, only the top 50 strongest negative correlations are shown for visualization purposes. Genes with green frames are involved in EMT and/or ECM-related processes. i.s.=incertae sedis; u.c.: uncultured.

Aging-specific changes in colonic gene expression between 12 and 28 months

Next, we investigated alterations in gene expression during the later phase of the aging process by analyzing the changes between 12 and 28 months of age. Since the number of differentially expressed genes (n=170, q<0.01) was too low to perform a functional analysis, a less stringent IBMT p-value<0.01 was applied to include genes displaying more subtle expression changes, resulting in 1,371 significantly differentially expressed genes (Figure 4A). The top 15 highest up- and down-regulated genes presented in Table 3 show that Fatty acid binding protein 6 (Fabp6) had the strongest increased expression. The top 15 down-regulated genes were dominated by immune- and inflammatory-related genes, including 8 immunoglobulin-coding genes, Angiogenin 4 (Ang4), regenerating islet-derived 3 beta (Reg3b) and Resistin like beta (Retnlb). Another interesting immune-related gene displaying age-related reduced expression, but with a smaller fold change ranking them out of the top-15 list, was Indoleamine 2,3-dioxygenase 1 (Ido1) (Data set S4).

Figure 4. Aging-specific changes in colonic gene expression between 12 and 28 months and integration with colonic microbiota composition. (A) Number of differentially expressed genes between 12vs28, 12vs24 and 24vs28 months (IBMT p<0.01). (B) Top 10 most significant canonical pathways between 12 and 28 months identified by Ingenuity Pathway Analysis (IPA). Numbers behind bars represent the number of differentially expressed genes (IBMT p<0.01) and numbers behind the pathways represent the total number of genes. (C) Correlations between colonic microbiota composition and gene expression during aging. The 1371 genes differentially expressed between 12 and 28 months were correlated with the 50 genera with a relative abundance ≥0.1% in at least 1 sample, including the mice aged 12 and 28 months. Relevance network showing the genera (white circles) and genes (grey rectangles) that were correlated (r<-0.80 or r>0.80). Red lines represent positive correlations, blue lines represent negative correlations. For visualization purposes, only the top 25 strongest positive top 25 negative correlations are shown. i.s.=incertae sedis; u.c.: uncultured.

The most significantly enriched down-regulated gene sets identified by GSEA were related to DNA replication and DNA synthesis (Table S3B). In addition, IPA identified ‘Protein Ubiquitination Pathway’, ‘Mismatch Repair in Eukaryotes’ and ‘Cell Cycle Control of Chromosomal Replication’ as most significantly regulated canonical pathways during the later phase of the aging process (Figure 4B) and X-box binding protein 1 (Xbp1) as the strongest inhibited upstream regulator (Table S3C). Taken together, these results indicate that during the later phase of aging (between 12 and 28 months), gene expression profiles pointed towards a diminished antimicrobial defense, an altered protein degradation response and aberrations in DNA repair mechanisms.

Integrative analysis of microbiota composition and gene expression in the later phase of aging

A second integrative analysis was carried out to investigate the correlation between gene expression and colonic microbiota composition during the late phase of life. For that purpose, the data of the 50 genera with a relative abundance of more than 0.1% in at least one sample and the 1,371 genes that were differentially expressed between 12 and 28 months (IBMT p<0.01) for each individual mouse were integrated. The highest number of most pronounced correlations (r<-0.80 or r>0.80) was found for an uncultured member of Gastranaerophilales spp. (Figure 4C). Pcdhb16, Pcdhb18 and Abcg5, all in the top 15 of strongest up-regulated genes (Table 3), as well as Ubiquitin Specific Peptidase 20 (Usp20), were negatively correlated with the relative abundance of this uncultured member. Interestingly, strong positive correlations (R≥0.75) were found for the uncultured member of Gastranaerophilales spp. with the down-regulated immune-related genes Ang4, Retnlb and Ido1 (Data set S5).

Markers of gut permeability increased during aging

To determine whether there were changes in markers of gut permeability, concentrations of lipopolysaccharide-binding protein (LBP), a marker for bacterial endotoxins, were assessed in plasma of mice aged 6, 12, 24 and 28 months. LBP concentrations were significantly higher (p<0.05) in plasma from mice aged 28 months compared to 6-months-old mice (Figure 5A). Besides, the microarray data of the colonic wall revealed increased expression levels of the Lbp gene during aging (Figure 5B). LBP is also synthesized as an acute-phase protein in the liver and therefore we assessed expression levels of Lbp in liver tissue, which were increased during aging (Figure 5C). Interestingly, although we found increased plasma markers of gut permeability at higher age, there was no collective increase in expression of tight junction genes (Data set S2).

Figure 5. Markers of gut permeability during aging. (A) Concentrations of Lipopolysaccharide-binding protein (LBP) (ng/ml) in plasma. (B) Expression (Log2) of Lbp in the colon and (C) liver. *p<0.05.

Ethics statement

The institutional and national guidelines for the care and use of animals were followed and the Local Committee for Care and Use of Laboratory Animals at Wageningen University approved the experiment (code number: drs-2010151b).

Mice and study design

In this study, 9-week old male C57BL/6J mice were housed individually and fed an ad libitum semi-synthetic (AIN-93W) diet. The study design was described previously in more detail [52]. In order to study the aging process, mice were randomly distributed into four groups that were sacrificed at 6, 12, 24 and 28 months (n=11-16 mice per group). For 4 mice in the 24-months group and 5 mice in the 28-months group, fresh faeces were consecutively collected (directly after defecation) at the age of 4, 8, 12, 18, 24 and 28 months. Bodyweight and food intake were monitored bi-weekly. During sacrifice, the colon was opened and the colonic luminal content was collected first. Then, the middle part of the colon (~1 cm) was turned inside out, rolled and embedded in paraffin (‘Swiss rolls’) for histology. From the remaining colon parts, the mucosa and submucosa were scraped for RNA isolation. The mice included in the present study were also used as control group for diet-intervention studies in the liver and colon of which the data have been reported in previous publications [21,52–54].

Gut microbiota composition analysis

For both fresh faeces (42 samples) and colonic luminal content (40 samples) the microbiota composition was determined. The detailed methods have been described previously [21]. Briefly, DNA was isolated from faecal and colonic luminal content samples using the ZR Fecal DNA MicroPrep kit (ZYMO Research, Irvine, CA, USA) according to the manufacturer’s instructions. Lysis of the samples was performed by bead beating using ZR BashingBead™ Lysis Tubes (ZYMO Research, Irvine, CA, USA). Next, the V3-V4 region of the 16S rRNA gene was amplified. After purification of the amplicon, DNA quality was checked and a second PCR was performed using sample-specific barcoded primers (Nextera XT index kit, Illumina, San Diego, CA, USA). Purified PCR products were sent to BaseClear service laboratory (BaseClear BV, Leiden, The Netherlands) for sequencing on the MiSeq platform (Illumina, San Diego, CA, USA). The Casava pipeline (version 1.8.3, Illumina, San Diego, CA, USA) was used for de-multiplexing of the FASTQ files. After initial quality control, three fresh faeces samples were excluded for further analysis, as the coverage was below threshold (50% of the median number of reads). CLC Microbial Genomics Module version 1.2.1 (CLC Bio, Qiagen, Aarhus, Denmark) was used for further analysis of the sequencing data. FASTQ files of the fresh faeces and colonic luminal content were imported and processed simultaneously. The identified OTUs were aligned against the Silva database (version 119) at 97% similarity [55]. A number of 1,034,303 reads were detected, for which 1208 OTUs were identified. 25 of these OTUs (265 reads) were unavailable. Shannon entropy was considered as a metric for alpha-diversity and the beta-diversity was determined using Bray-Curtis distances. Principal Component Analyses and Redundancy analysis were carried out using Canoco 5 [56]. The input consisted of the 303 OTUs with a relative abundance of ≥0.1% in at least one sample. With regard to the RDA, the numbers of reads were centered and standardized using the Hellinger transformation [57]. The 16S rRNA gene data described in this study have been deposited in the sequence read archive at the NCBI with accession number SRP145060.

Determination of metabolites in colonic luminal content

The colonic luminal content of mice aged 6 months (n=6), 24 months (n=8) and 28 months (n=8) was prepared as described previously [21] and levels of metabolites were determined using ¹H-NMR. From the aligned spectra, integrals for SCFAs were identified. Besides, an untargeted approach was taken to identify metabolites that differed between age groups. After baseline correction, the concentrations of the metabolites were calculated by taking into account the dilution factor and number of hydrogen atoms. Correlations between concentrations of identified metabolites and the relative abundances of the 50 most abundant genera in colonic luminal content were determined using Spearman’s correlation.

Gene expression analysis

For a number of 4 to 8 mice per age group, RNA was isolated from scrapings of the colonic wall and liver using TRIzol reagent (Invitrogen, Breda, The Netherlands). Isolated RNA was purified using RNeasy Micro columns (Qiagen, Venlo, The Netherlands) and total RNA yield (Nanodrop ND-1000, Nanodrop Products, Maarssen, The Netherlands) and RNA integrity (Agilent 2100 Bioanalyzer, Agilent Technologies, Amsterdam, The Netherlands) were measured. Only RNA was used that had a RNA integrity number (RIN) above 8.0 (Table S4). Purified RNA (100 ng per sample) was converted to cDNA and labelled using an Ambion WT expression kit (Life Technologies, Bleiswijk, The Netherlands). Microarray hybridization and analysis was performed as previously described [52]. Differences in colonic gene expression between the age groups were analyzed using the Intensity Based Moderated T statistics (IBMT), using q-values <0.01 as threshold for the comparison 6 versus 12 months and p-values <0.01 for the comparison 12 versus 28 months. Microarray data has been submitted to the Gene Expression Omnibus (GEO) at the NCBI, and is accessible under number GSE113257. The Principal Component Analysis plot was generated using MultiExperiment Viewer version 4.9.0. Ingenuity Pathway Analysis (IPA) was used for identification of canonical pathways and upstream regulators [58]. Gene Set Enrichment Analysis was performed to functionally interpret gene regulation [59]. A false discovery rate (FDR) q-value of <0.1 was considered significantly enriched. Short Time-series Expression Miner (STEM) version 1.3.11 [60] was used to identify gene expression profiles during aging, using the mean Robust Multichip Average (RMA) values as input.

Histology

Swiss rolls of 6 mice per age group were cut in 5 micrometer sections and mounted on microscope slides (VWR Superfrost Plus Micro Slide). The sections were dewaxed multiple times with xylene, rehydrated in alcohol and stained with Haematoxylin/Eosin (H&E) or Sirius red/Fast green. Pictures of the stained sections were taken using a microscope with Olympus CellSense Entry software (Olympus Europe, Hamburg, Germany). For the H&E stained sections, the depths of the colonic crypts were measured using Olympus Cell^B software.

Multivariate analyses for gut microbiota and gene expression

In order to integrate colonic microbiota composition and colonic gene expression, we performed multivariate analyses in R 3.3.1 using the MixOmics package [61]. Partial Least Squares (PLS) regression in the canonical mode was used to determine the bi-directional relationship between the two datasets. We used data from individual mice for whom both gene expression as well as the microbiota composition from colonic luminal content were available. For the MixOmics analysis early during aging (i.e. between 6 and 12 months of age), a number of 11 mice were included (6 months: n=4; 12 months: n=7). A number of 14 mice were included for the analysis using the genes differentially expressed between 12 and 28 months (p<0.01) (12 months: n=7; 28 months: n=7). The genera with a relative abundance of ≥0.1% in at least 1 sample were included, resulting in a number of 50 genera. First, the sequencing reads of the genera were transformed and zeros were imputed based on the centered log ratio (clr) using the aldex2 package with 1024 Monte Carlo permutations [62]. Gene expression data were included after transformation into log2 RMA intensities. Relevance networks with the strongest positive and negative correlations (R<-0.8 or R>0.8) were visualized using Cytoscape [63].

Plasma measurements

Plasma levels of LBP were determined using the LBP ELISA Kit Mouse (kit OKBB00573, Aviva Systems Biology, San Diego, USA), following the manufacturer’s protocol. Plasma samples were diluted 100 times. Plasma levels of insulin were determined using the Mouse Adipokine kit (MADKMAG-71K, Merck Millipore, Darmstadt, Germany), according to the manufacturer’s protocol.

Statistical analyses

Statistical analyses were, unless stated otherwise, carried out using GraphPad Prism 5.04 (GraphPad Software, San Diego, California, USA). The Kolmogorov-Smirnov test was used to test if data were normally distributed. In case of normal distribution, one-way analysis of variance (ANOVA) in combination with Tukey’s multiple comparisons test was used to test for differences between ages. In case the data were not normally distributed, the non-parametric Kruskal-Wallis one-way analysis of variance following a Dunn’s post hoc test and a Mann-Whitney test were used for independent samples. A Friedman test or Wilcoxon Signed Rank test were used in case samples were dependent. The Bonferroni or Benjamini-Hochberg approach were used to adjust for multiple testing, if appropriate. For all statistical tests, unless stated otherwise, a p-value of p<0.05 was considered as statistically significant.