Effects of storage and technical variation
We first validated our methods by assessing the effect of storage and technical variation on microbiome composition. To quantify the effect of the two storage methods on bacterial composition in fresh samples, we performed a separate pilot study with nine faecal samples sourced from nine captive meerkats at Zurich University. Samples were immediately frozen after collection, and then either freeze-dried or kept frozen at −80 °C for seven days. Microbiome composition clustered strongly by sample identity in their beta diversity (Supplementary Fig. 1b), and storage did not significantly affect composition (Weighted Unifrac: F = 0.7, p = 0.52; Unweighted Unifrac: F = 1.0, p = 0.37). Across samples analysed in this study, storage had significant yet small effects on estimated bacterial load, with frozen samples overall having slightly lower estimated abundance (t = 7.2, p < 0.001, R2 = 0.04; Supplementary Fig. 1c). Observed ASV richness did not differ between storage types (t = 0.7, p = 0.48; Supplementary Fig. 1d). Storage had weak but significant effects across four measures of beta diversity (p < 0.001, R2 = 0.01–0.02; Supplementary Fig. 1e). Because storage had small effects on the measured composition, we account for storage in all models, and we only consider associations robust if they exhibit significant trends across both frozen and freeze-dried samples.
We tested how micro-variation in weighing and other sources of technical variation affected estimated load and diversity measures by including 16 extraction replicates, which were treated separately at every stage of processing. Technical variation accounted for 10% of variation in estimated bacterial load (Supplementary Fig. 2a), and 1–2% of variation across measures of alpha diversity (Supplementary Fig. 2b) and the first axis of variation of four measures of beta diversity (Supplementary Fig. 2c–f). Whilst technical variation was non-negligible, sample ID accounted for 90–98% of variation across measures.
Time of day is the strongest predictor for bacterial load and diversity
We first aimed to investigate how estimated bacterial load, alpha diversity, and beta diversity change over the diurnal, seasonal, and lifetime scales. We modelled bacterial load and alpha diversity across the three time scales by fitting generalised additive mixed models (GAMMs) to log-transformed bacterial load and (untransformed) ASV richness, applying non-linear smoothing functions to time of day, month, and meerkat age, whilst controlling for sampling depth, sequencing run, storage method, and time in field conditions prior to being frozen as fixed effects, and including individual ID and social group as random effects. Definitions of all terms included in models are outlined in Supplementary Table 1.
Mean bacterial load underwent the largest shifts across the day, in comparison to seasonal and lifetime scales, which were both much weaker (Hours after sunrise: F = 54.4, p < 0.0001; Month: F = 1.1, p = 0.007; Age: F = 9.1, p = 0.003; model R2 = 0.47; Supplementary Table 2). Bacterial load tended to be highest early in the morning and lowest approximately 10 h after sunrise (Fig. 2a), although it should be noted there is considerably uncertainly regarding estimates for the middle of the day when sampling is sparse. Bacterial load fluctuated only weakly with season (Fig. 2b) and age (Fig. 2c). Whilst seasonal and lifetime shifts in bacterial load were weak but significant across the full dataset, they were not replicable across both frozen and freeze-dried samples (Fig. S3a). Sequentially excluding terms to assess the drop in model explanatory power (R2) indicated that diurnal, seasonal, and lifetime dynamics accounted for 12%, 1%, and 1% variation in bacterial load, respectively. Individual ID and social group together accounted for 2% of variation, whilst methodological variables accounted for 31%.
a–f Smoothed estimates and partial residuals from two GAMMs predicting (a–c) bacterial load and (d–f) observed ASV richness across the day, year, and life. Shaded area represents 95% CIs. g–i Beta diversity MDS ordination coloured by (g) hours after sunrise (orange = morning; afternoon = purple), (h) season (green = wet season; yellow = dry season), and (i) age (blue = young; grey = adult; red = old). For clarity, samples were grouped into discrete categories to generate 95% CI ellipses and group centroids (large circles), which represent samples collected in the morning field session (< 7 h after sunrise) and afternoon (> 7 h after sunrise), wet (October–April) and dry (May–September) season, and age categories (< 1 years, 1–5 years, and >5 years). Source data are provided in the source data file.
ASV richness also demonstrated the strongest fluctuations across the day (Hours after sunrise: F = 20.8, p < 0.0001; Month: F = 3.4, p < 0.0001; Age: F = 3.8, p = 0.072; model R2 = 0.3; Supplementary Table 3). In contrast to bacterial load, observed richness was lowest in the early morning and evening, and highest in the middle of the day (Fig. 2d). Across the year, mean ASV richness increased over the wet summer, and declined over the dry winter (Fig. 2e). ASV richness did not exhibit significant changes across life (Fig. 2f). These results were consistent across frozen and freeze-dried samples (Supplementary Fig. 3b). Overall, our model accounted for 30% of variation in ASV richness. However, diurnal, seasonal, and lifetime dynamics accounted for only 6%, 3% and 1%, respectively; meerkat ID and social group accounted for 4%, whilst methodological variables accounted for 16%. We tested whether these trends were also evident when applying the abundance-weighted Shannon diversity, and found weaker and inconsistent trends across all temporal scales (Supplementary Fig. 3c). Thus, changes to observed ASV richness across the day were driven by rare taxa.
We explored shifts in beta diversity across the three temporal scales by ordinating taxa composition using Multi-Dimensional Scaling (MDS) analysis of Weighted Unifrac distances, which accounts for both abundance and phylogeny. Community composition along the first two ordination axes clustered most strongly by the time of day (Fig. 2g), and only very weakly by season (Fig. 2h) or age (Fig. 2i). A PERMANOVA test on the weighted Unifrac distance matrix indicated that time of day had the strongest effect size yet explained relatively little variation (Hours after sunrise: F = 49.7, R2 = 0.038, p < 0.001; Month: F = 3.1, R2 = 0.002, p = 0.003; Age: 5.28, R2 = 0.004, p < 0.001; Supplementary Table 4a). Meerkat ID explained a large proportion of variation yet had a small effect size (i.e. centroids are close together; F = 1.1, R2 = 0.19, p = 0.032), whilst social group was not significant (F = 0.9, R2 = 0.03, p = 0.88). Similar results were generated when using Unweighted Unifrac (Supplementary Table 4b) and across storage types (Supplementary Fig. 3d, e).
Genus-level dynamics
We next aimed to identify which genera were influencing shifts in beta diversity and to model the dynamics of genus-level abundances across temporal scales. Changes to beta diversity across the day reflected a decrease in the relative abundance of the genus Clostridium between morning and afternoon (Fig. 3a). Axes 1 and 2 of the Weighted Unifrac ordination shown in Fig. 2g–i largely represented the continuum of Clostridium and Bacteroides abundances, respectively, and sample placement along these two axes was strongly influenced by time of day (R2 = 0.27, p < 0.001), but only very weakly by season (R2 = 0.004, p = 0.02) and not by age (R2 = 0.001, p = 0.48). Samples taken in the morning were more likely to be dominated by Clostridium than those taken in the afternoon, which tended to be dominated by a more diverse suite of Raoultibacter, Cellulomonas, Bacillicaceae, Enterococcus and Lactococcus (Fig. 3b). In contrast, axes 3 and 4 largely represented the continuum of Paeniclostridium and Blautia abundances, respectively, and sample placement along these axes was more associated with age (R2 = 0.043, p < 0.001) and weak seasonal effects (R2 = 0.007, p = 0.003), rather than time of day (R2 = 0.003, p = 0.22; Fig. 3c). The microbiomes of old meerkats (>5 years old) more likely to be dominated by Romboutsia, and Paeniclostridium, than the microbiomes of adults (1–5 years old) and young meerkats (< 1 year old; Fig. 3c).
a Summary of taxonomic shifts in relative abundance at the genus-level per 30-min interval from sunrise. The number of samples that were summarised per 30-min slot are indicated by the histogram. b, c Weighted Unifrac ordination plots of (b) axes one and two and (c) three and four, coloured and grouped by the most abundant genus in each sample. Large circles represent group centroids for samples sharing the same most abundant genus. Arrows indicate the direction and influence of significant temporal variables when categorised into morning/afternoon, wet/dry seasons, and young/adult/old. Statistics for temporal variables (arrows) are shown. Source data are provided in the source data file.
Changes to beta diversity are driven by shifts in relative abundance rather than absolute abundance. We therefore modelled genus-level dynamics in absolute abundance over daily, seasonal and lifetime scales. We first performed simple differential abundance non-parametric tests across all genera with over 15% prevalence across samples (n = 117) to identify genera that were differentially abundant in the morning compared to afternoon, in the dry season compared to the wet season, young meerkats versus adults, and adult meerkats versus old meerkats (Supplementary Fig. 4). Almost all genera were significantly associated with time of day (Supplementary Fig. 4a), suggesting that diurnal oscillations are widespread across gut microbiome members. Only a few genera significantly differed between dry and wet seasons (Supplementary Fig. 4b). A small number of genera were differentially abundant in adults compared to young meerkats (Supplementary Fig. 4c), whilst none were differentially abundant in old meerkats compared to adults (Supplementary Fig. 4d).
We next focused on 16 notable genera in order to model their temporal dynamics using GAMMs whilst controlling for potentially confounding methodological variables. We focused on the most prevalent and abundant genera (n = 13) which all had at least 60% prevalence across samples and together accounted for 75% relative abundance. However, we used the results from the differential abundance analysis to select three additional rarer genera that exhibited notable trends for additional analysis, including Raoultibacter (43% prevalence), and Callulomonas (38% prevalence). We also include a particularly rare genus, Eubacterium (18% prevalence), which was only present in young individuals.
Thirteen of the 16 focus genera significantly shifted over the day (Fig. 4a, b), whilst five genera fluctuated with season (Fig. 4c, d), and six were associated with age (Fig. 4e, f). Relative to their mean abundance, most genera followed the same diurnal pattern, with the highest abundance in the early morning, and decreasing until the late afternoon approximately 10 h after sunrise, before increasing slightly again prior to sunset (Fig. 4a, b). Clostridium exhibited the strongest diurnal oscillations, and made up a large proportion of bacterial load in the morning (Fig. 4b). Raoultibacter, Cellulomonas, and a Bacillaceae genus increased in the afternoon as the abundances of most genera fell, yet in general still maintained relatively low abundances (Fig. 4b). Across the year, Clostridium and Blautia increased over the dry (winter) season, whilst Raoultibacter, Cellulomonas, and the Bacillaceae genus increased over the wet (summer) season (Fig. 4c, d). Over meerkat life, six genera, including Peptococcus, Paeniclostridium, Romboutsia, and Lachnoclostridium, increased over the first 1–2 years of life before levelling off (Fig. 4e, f). In contrast, Eubacterium decreased over the first 6 months of age. There was no evidence for changes to any taxa in old age (Fig. 4c), although we note that there was a tendency for old individuals to have reduced abundance of Christensenellaceae (Supplementary Fig. 4, Supplementary Fig. 5). Trends for each genus individually, including 95% confidence internals and split by storage, are presented in Supplementary Fig. 6 (diurnal trends), Supplementary Fig. 7 (seasonal trends), and Supplementary Fig. 8 (lifetime trends).
Comparison of temporal dynamics of 16 focus genera across diurnal (a, b), seasonal (c, d), and lifetime (e, f) scales. Top panel shows GAMM abundance estimates across the three temporal scales compared to the mean, where zero (indicated by the dashed line) represents the mean (log10) abundance of each genus. Estimates have been back-transformed in the bottom panel to represent absolute abundance, and bacterial load (grey) is shown for comparison. Only genera that significantly shift across both frozen and freeze-dried samples are shown (see Supplementary Figs. 6–8 for 95% confidence intervals and trends split by storage). Source data are provided in the source data file.
Overall, a summary of effect sizes across diurnal, seasonal, and lifetime scales for all models presented thus far is presented in Fig. 5a, showing that diurnal oscillations are generally stronger, more prevalent, and more robust than seasonal and lifetime trends for both community and single-genus phenotypes.
a Effect sizes for non-linear temporal variables on bacterial load, diversity measures, and abundances of 16 focus genera. b Effect sizes and 95% confidence intervals for eight environmental, biological, and behavioural variables. Dark blue/red denotes significant associations that are consistent across storage treatments, whilst light blue/red denotes significant yet inconsistent across storage. c Partitioning of model R2 into variance explained by temporal, mechanistic, and methodological variables. Source data are provided in the source data file.
Microbiome dynamics is associated with temperature and foraging schedules
Our third aim was to investigate the biological mechanisms underpinning gut microbial dynamics. We were particularly interested in exploring how well climate and foraging schedules explained diurnal dynamics, which climatic variables best explain seasonal changes to the microbiome, and whether body condition and social status underpin changes to genera associated with ageing. We therefore examined the potential mechanisms underpinning temporal dynamics in five measures of bacterial load and diversity and abundance of the 16 focus genera by building GAMMs incorporating eight environmental, biological, and behavioural variables. These mechanistic variables were cumulative rainfall over the previous month, maximum and minimum temperature on the day of sample collection, the temperature at the time of sample collection, meerkat sex, body condition at the time of sampling, social status (dominant/subordinate), and the number of hours spent foraging at the time of sample collection (see Supplementary Table 1 for definitions). Foraging schedules were estimated from long-term observational data across the year, and consisted of a foraging period in the morning and a foraging period in the afternoon which shifted across the year (Supplementary Fig. 9). To summarise results, we present effect sizes of (non-linear) temporal trends with mechanistic variables excluded (blue points, Fig. 5a), and effect sizes for just mechanistic variables with temporal smooths excluded (red points, Fig. 5b), categorising associations by how robust they are to methodology.
As demonstrated in Figs. 2–4, diurnal dynamics were much stronger than seasonal and lifetime scales for both diversity measures and genus-level abundances, and were also more robust to splitting the dataset by sample storage (Fig. 5a). When considering mechanistic variables, diurnal oscillations of bacterial load, alpha and beta diversity, as well as the abundance of most genera, were largely associated with temperature and foraging schedule (Fig. 5b). Most genera, and notably Clostridium, increased over the foraging period, and decreased as diurnal temperatures climbed. Conversely, most diurnally oscillating taxa also increased with maximum temperature. Together, these indicate that bacterial load and abundances of Clostridium, Paeniclostridium, and Romboutsia, as well as others, peak early in the morning on hot days, likely tracking foraging schedule. Rainfall had little effect on microbiome composition, and neither did minimum temperature. Bacillaceae, Roultibacter, and Cellulomonas, in contrast, were associated with low maximum temperatures and high sampling temperature, and decreased with foraging.
We found little robust evidence for sex, body condition, or social status playing a strong role, although there were a number of weak signals across the whole dataset, albeit not replicable between frozen and freeze-dried samples. Abundances of Paeniclostridiam and Romboutsia, both of which increase with age, tended to be more abundant in dominant individuals than subordinates. Yet, these results were not replicable, potentially due to reduced statistical power in split datasets combined with the relatively few samples from dominant individuals (201 dominant vs 908 subordinate). Moreover, we found little evidence for associations between the gut microbiome and body condition at the time of sampling, although, again, a number of genera such as Bacteroides were weakly and negatively associated with the condition when considering all data together.
We partitioned model variance by excluding temporal and mechanistic terms from the global model that included all terms (Fig. 5c). Models explained 15–44% in variation in abundance, including methodological terms. Whilst mechanistic variables accounted for part of the variation in microbiome dynamics, a substantial proportion was only explained by temporal variables, suggesting that additional process (e.g. light-dark cycles) are at least partially responsible for diurnal oscillations. Overall, temporal and mechanistic variables accounted for approximately 5–20% in variation across measures.
Diurnal oscillations do not decay with age
Finally, we tested whether genus-level oscillations decrease in magnitude in older individuals by splitting samples into three age categories: samples taken from meerkats under 1 year old (n = 385), from adults between 1- and 5 years old (n = 627), and samples from old meerkats that were over 5 years of age (n = 97). Patterns in diurnal oscillations were very similar across all age categories, with Clostridium sensu stricto 1 demonstrating the strongest oscillations across all groups (Fig. 6a–c). We looked in closer detail at five genera (Clostridium, Bacteroides, Paeniclostridium, Cellulomonas, and Raoultibacter) that demonstrated the strongest non-linear diurnal trends and found that confidence intervals around the mean overlapped across the age categories (Fig. 6d–h). Therefore, we found no evidence for a reduction in circadian rhythms with age.
a–c Diurnal oscillations of 16 focus genera in a) young meerkats (<1 year old; n = 385). b adult meerkats (1–5 years old; n = 627); and (c) particularly old meerkats (>5 years old; n = 97). Zero represents the taxa mean and sample distributions are indicated by the histograms. d–h Model estimates and 95% confidence intervals for (d) Clostridium (sensu stricto 1); (e) Bacteroides; (f) Paeniclostridium, (g) Cellulomonas, and (h) Raoultibacter, split by age category (blue = young; grey = adult; red = old). Source data are provided in the source data file.
Source: Ecology - nature.com