Introduction

For most vertebrate species, physical declines with age are inevitable. These changes define the concept of biological aging and contribute to increased disease burden in older individuals (1, 2). While mean patterns of biological aging may be species-typic, the pattern of biological aging differs across individuals within species. Hence, an animal’s age in years—i.e., its chronological age—is not an exact reflection of age-related decline in physical functioning (3–6). Measuring individual differences in biological age is an important first step to understand how socio-environmental conditions influence aging processes and to identify strategies to improve health in old age.

One valuable marker of biological aging may lie in the composition and dynamics of the mammalian gut microbiome (7–10). Age-related changes in gut microbiota are well documented in humans and other animals, and the gut microbiome has the potential to reflect a wide variety of aging processes for individual hosts (11–24). Mammalian gut microbiota interact with the immune, endocrine, nervous, and digestive systems, all of which change with age (25–28). Gut microbiota are also sensitive to host environments and behaviors that change with age, including host diet, living conditions, and social integration (23, 29–32). Finally, gut microbiota may play a causal role in age-related changes in host development and longevity and may, therefore, be directly involved in individual differences in biological age (16, 17, 22, 33, 34). For example, children who experience famine exhibit developmentally immature gut microbiota that, when transplanted into mice, delay growth and alter bone morphology (35–38). Experiments in flies, mice, and killifish find that the gut microbiome can also influence longevity (16, 17, 22, 39).

One strategy for testing if gut microbiota reflect host biological age is to apply supervised machine learning to microbiome compositional data to develop a model for predicting host chronological age, and then to test whether deviations from the clock’s age predictions are explained by socio-environmental drivers of biological age and/or predict host development or mortality. A parallel approach is commonly applied to patterns of DNA methylation, and the resulting epigenetic clocks predict disease and mortality risk more accurately than chronological age alone (40–45). To date, five microbiome age-predicting clocks have been built for humans, which predict sample-specific age with median error of 6 to 11 years (46–49). However, to our knowledge, no clocks have tested whether microbiome age is sensitive to potential socio-environmental drivers of biological age or predicts host development or mortality.

Here we create a microbiome age-predicting clock using 13,476 16S rRNA gene sequencing-based gut microbiome compositional profiles from 479 known-age, wild baboons (Papio sp.) sampled over a 14-year period (Figures 1A and 1B). These microbiome profiles represent a subset of a data set previously described in Grieneisen et al. (50), Björk et al. (51) and Roche et al. (52), filtered to include only the baboon hosts whose ages were known precisely (within a few days’ error). Important to human aging, baboons share many developmental similarities with humans, including an extended juvenile period, followed by sexual maturation and non-seasonal breeding across adulthood (Figures 1C and 1D; [51–54]).

Figure 1.

The baboons in our data set were members of the well-studied Amboseli baboon population in Kenya, which has continuous, individual-based data on early life environments, maturational milestones, social relationships, and mortality (57–62). Relevant to measuring biological age, prior research in Amboseli has identified several demographic, environmental, and social conditions that predict physical condition, the timing of development, or survival: sex, season, social status (i.e., dominance rank), and early life adversity (53–55, 63–69). Consistent with the possibility that these associations arise from causal effects of harsh or stressful conditions on biological aging (e.g. 70–72), and with the idea that microbiomes also serve as a marker of biological age (7–10), we tested whether conditions associated with poor physical condition, delayed development, or higher mortality are also associated with microbiome age estimates. Our predictions about the direction of these effects varied depending on the socio-environmental condition and the developmental stage of the animal (i.e., juvenile or adult). In terms of sex, adult male baboons exhibit higher mortality than adult females. Hence, we expected them to exhibit microbiomes that are old-for-age, compared to females (we expected no sex effects on microbiome age in the juvenile period).

In terms of season, the Amboseli ecosystem is a semi-arid savannah with a 5-month long dry season during which little rain falls, often leading to nutritional hardship (75). We expected that samples from the dry season might appear to be old-for-age, compared to samples from the wet season due to nutritional stress in this difficult season.

In terms of social status, baboons experience linear, sex-specific hierarchies. Female ranks are nepotistic, with little social mobility, and low-rank is linked to low priority of access to food (55, 68, 76–78). In contrast, adult male rank is determined by strength and fighting ability and is dynamic across adulthood (79). High-ranking males experience high energetic costs of mating effort, have altered immune responses, and exhibit old-for-age epigenetic age estimates compared to low-ranking males (45, 64, 65). We expected that individuals who pay the largest energetic costs—low-ranking adult females and high ranking adult males (45)—would appear old-for-age.

In terms of early life adversity, prior research in Amboseli has identified six conditions whose cumulative, and sometimes individual, effects predict adult female mortality, including maternal loss prior to age 4 years, drought in the first year of life, birth into an especially large social group, the presence of a close-in-age competing younger sibling, and having a low-ranking or socially isolated mother (64, 65, 67). For adult female baboons, experiencing multiple sources of adversity in early life is the strongest socio-environmental predictor of baboon mortality in Amboseli; hence, we expected that these individuals would have old-for-age clock estimates in adulthood (58, 59, 64, 66, 67, 80). However, we also expected that some sources of early life adversity might be linked to young-for-age gut microbiota. For instance, maternal social isolation might delay gut microbiome development due to less frequent microbial exposures from conspecifics.

We began our analyses by identifying microbiome features that change reliably with host age. We then constructed a microbiome clock by comparing the performance of several supervised machine learning algorithms to predict host age from gut microbial composition in each sample from each host. We evaluated the clock’s performance for male and female baboons and tested whether deviations from clock performance were predicted by the baboons’ social and environmental conditions (guided by the predictions outlined above). Lastly, we tested whether baboons with young-for-age gut microbiotas have correspondingly late developmental timelines or longer lifespans. In general, we expected that adult baboons who experienced harsh conditions, especially adversity in early life, which is the strongest socio-environmental predictor of baboon mortality in Amboseli known to date, would appear an old-for-age (58, 59, 64, 66, 67, 80). Alternatively, microbiome age might be better predicted by an individual’s current environmental or social conditions (e.g., season or social status), rather than past events. Such results would support recency models for biological aging (79, 80) and would be consistent with findings from a recent epigenetic clock study in Amboseli (45). In support of this alternative perspective, we find that season and social rank have stronger effects on microbiome age than early life events. Further, microbiome age does not predict host development or mortality. Our work highlights important ways that social and environmental conditions shape microbiome aging in natural systems.

Results

Many microbiome features change with age

Before creating the microbiome clock, we characterized microbiome features that change reliably with the age of individual hosts. Our subjects were 479 known-age baboons (264 females and 215 males) whose microbiota were characterized using 13,476 fecal-derived 16S rRNA gene sequencing profiles over a 14-year period (Figures 1A and 1B; baboon age ranged from 7 months to 26.5 years; 8,245 samples from females; 5,231 samples from males; range = 3 to 135 samples per baboon; mean = 35 samples per female and 26 samples per male).

We tested age associations for 1,440 microbiome features, including: (i) five metrics of alpha diversity; (ii) the top 10 principal components of Bray-Curtis dissimilarity (which collectively explained 57% of the variation in microbiome community composition); and (iii) centered log-ratio transformed abundances of each microbial phylum (n = 30), family (n = 290), genus (n = 747), and amplicon sequence variance (ASV) detected in >25% of samples (n=358) (Table S1). For each feature, we tested its association with host age by running linear mixed effects models that included linear and quadratic effects of host age and three other fixed effects known to explain variation in microbiome composition in our population: the season of sample collection (wet or dry), the average maximum temperature for the month prior to sample collection, and the total rainfall in the month prior to sample collection (50, 51, 62). Baboon identity, social group membership, and the hydrological year of sampling were modeled as random effects.

We found that many aspects of microbiome community composition changed with host age (Figure 2; Figure S1). All alpha diversity metrics, except richness, exhibited U-shaped relationships with age, with high values in early life and old age, and low values in young adulthood. While we should interpret this pattern with caution due to the small sample size beyond age 20 (n=18 females), this U-shaped pattern differs somewhat from patterns in humans and chimpanzees: most human populations exhibit an asymptotic increase in alpha diversity with age (81, 82) while in chimpanzees alpha diversity is highest in early life (23) (FDR < 0.05; Figures 2C and 2E; Table S1). Further, seven of the ten principal components (PCs) of microbiome composition exhibited linear, and in some cases quadratic, relationships with age, with PC1, PC2, PC4 and PC6 exhibiting the strongest age-associations (FDR < 0.05; Figures 2C and 2F; Table S1).

Figure 2:

In terms of individual taxa, 51.6% exhibited significant linear or quadratic relationships with host age (Figure S1 and S2; Table S1; FDR < 0.05). 60% of phyla (18 of 30) decreased proportionally with age, while only three phyla—Kiritimatiellaeota, Firmicutes, and Chlamydiae— increased proportionally with age (FDR < 0.05; Table S1). Similarly, 66% (66 of 100) of age-associated families and 55% (115 of 209) of age-associated genera exhibited declining proportions with age (Table S1). Consistent with the idea that age-related taxa differ across host populations and host taxa, none of the taxa that changed in this baboon population were commonly age-associated in humans (81) The taxa most consistently linked to human aging include Akkermansia, Faecalibacterium, Bacteroidaceae, and Lachnospiraceae (34, 81) while in our sample of baboons, the strongest age-related changes were seen in the families Campylobacteraceae, Clostridiaceae, Elusimicrobiaceae Enterobacteriaceae, Peptostreptococcaceae, and an uncharacterized family within Gastranaerophilales (Figure 2D and 2F; Table S1). The genera that had the strongest relationships with age included Camplylobacter, Catenibacterium, Elusimicrobium, Prevotella, Romboutsia, and Ruminococcaceae UCG-011 (Figure 2D and 2F; Table S1).

Microbiome clock calibration and composition

We next turned our attention to building a microbiome clock using all 9,575 microbiome compositional and taxonomic features present in at least 3 samples (Table S2). We did not restrict the features in the clock to age-associated taxa (Table S1) as the purpose of the clock is to test whether there are signatures of age in the microbiome, and to identify the set of key microbiome features that contribute to age prediction. These could include features that are not strongly age-correlated in isolation. In developing the clock, we compared the performance of three supervised machine learning methods to predict the chronological age of individual hosts at the time each microbiome sample was collected. The three machine learning methods were elastic net regression, Random Forest regression, and Gaussian process regression (see Supplementary Methods and Results). Because gut microbiota are highly personalized in ours and other host populations (51), at least one sample from each host was always present in the training sets for these models (see Methods).

We found that the most accurate age predictions were produced by a Gaussian process regression model with a kernel customized to account for heteroscedasticity (Figure 3; Figure S3). This model predicted host chronological age (agec), with an adjusted R² of 48.8% and a median error of 1.96 years across all individuals and samples (Figure 3A, Table 1). As has been observed in some previous age clocks (e.g., 39, 44, 46), microbial age estimates (agem) were compressed relative to the x=y line, leading the model to systematically over-predict the ages of young individuals and under-predict the ages of old individuals (Figure 3A).

Figure 3.

Table 1.

When we subset our agem estimates by sex, we found that the microbiome clock was slightly more accurate for males than for females (Figure 3B, Table 1). The adjusted R² for the correlation between agec and agem for males was 50.0%, with a median prediction error of 1.71 years as compared to an adjusted R² of 48.9% and median error of 2.15 years for female baboons (Table 1). Male baboons also exhibit significantly older gut microbial age than females (Figure 3B, chronological age by sex interaction: β=0.18, p<0.001, Table S3). Across the lifespan, males show a 1.4-fold higher rate of change in agem as a function of agec compared to females (relationship between agec and agem in males: β=0.63, p< 0.001; relationship between agec and agem in females: β = 0.45, p < 0.001; Table S4). Similar to patterns from a recent epigenetic age-predicting clock developed for this population (45), this effect was only present after sexual maturity: when we subset the agem estimates to microbiome samples collected prior to the median age at sexual maturity (5.4 years for testicular enlargement in males and 4.5 years for menarche in females [54]), we found no significant interaction between sex and age (agec by host sex interaction prior to median age of maturity: β= -0.09, p=0.203; agec by host sex interaction after median age of maturity: β=0.15, p<0.001; Table S3). After maturity, we recapitulate the 1.4-fold higher rate of change in males compared to females observed in the full data set (relationship between agec and agem in males only: β=0.53, p<0.001; relationship between agec and agem in females only: β=0.38, p<0.001; Table S4).

Overall, agem estimates performed reasonably well compared to other known predictors of age in the Amboseli baboons (Table S5; [44]). Agem was a stronger predictor of host chronological age than body mass index (BMI; except juvenile male BMI), blood cell composition from flow cytometry, and differential white blood cell counts from blood smears (Table S5). However, agem was a less accurate predictor of chronological age than both dentine exposure (males, females respectively: adjusted R² =73%, 85%; median error = 1.11 years, 1.12 years; Table S5) and an epigenetic clock based on DNA methylation data (males, females respectively: adjusted R² = 74%, 60%; median error = 0.85 years, 1.62 years; Table S5; [44]).

Social and environmental conditions predict variation in microbiome age

To test whether deviations in microbiome age for a given chronological age are correlated with socio-environmental predictors of health and mortality risk, we calculated whether microbiome age estimates from individual samples were older or younger than their hosts’ known chronological ages (⃤age; Figure 3C). We then tested whether several social and environmental variables predicted individual variation in microbiome ⃤age (Table S7; note that whether microbiome ages are old- or young-for-age is correlated with host age, hence our models always included host chronological age as a covariate). Overall, we expected that adult baboons who experienced harsh conditions in early life adversity (the strongest socio-environmental predictor of adult mortality in Amboseli) would tend to look old-for-age based on the microbiome clock (58, 59, 64, 66, 67, 80). Alternatively, microbiome deviations from chronological age might be best predicted by an individual’s current social status or season, rather than past events. These results would support recency models of biological aging (73, 74) and would be consistent with findings from a recent epigenetic clock study in Amboseli (45).

We found that individual baboons varied considerably in gut microbiome △age. For instance, in mixed effects models, individual identity explained ∼25% to ∼50% of the variance in △age for females and males respectively over the course of their lives (Table S6). Further, we found that season, dominance rank, and some aspects of early life adversity (large group size, early life drought, and maternal social isolation) were linked to small deviations from chronological age. In support of our expectation that microbiome samples collected in the dry season are old-for-chronological age, we found that age estimates based on microbiome samples collected from female baboons in the dry season were ∼2 months older than the host’s true chronological age (β= -0.180, p=0.021, Table 2; Table S6C; Figure S6A). However, season did not significantly predict the difference between microbiome age and known age in male baboons.

Table 2.

In terms of social status, we expected to observe that low-ranking females and high-ranking males would be old-for-chronological age (45). In support, we found that estimates from high-ranking males were old-for-age compared to estimates from low-ranking males, but these effects were relatively weak and noisy (rank effect: β=0.033, p<0.001; Figure 4A; Table 2; Tables S6B and D). Specifically, controlling for chronological age, alpha male gut microbiomes (ordinal rank = 1) appeared to be approximately 4 months older than microbiomes sampled from males with an ordinal rank of 10 (Table 2). However, contrary to our expectations, high-ranking female baboons also had old-for-age estimated when compared to low-ranking females (rank effect: β=1.745, p<0.001; Figure 4B; Table 2; Tables S6A and C). Specifically, controlling for chronological age, the microbiome of an alpha female (proportional rank=1) appeared to be approximately 1.75 years older than the lowest-ranking females in the population (proportional rank=0; Table 2).

Figure 4.

Some forms of early life adversity also predicted variation in microbiome △age, but only in males, and in inconsistent directions. For instance, males born into the highest quartile of observed group sizes had old-for-age estimates. Males experiencing this source of early life adversity had gut microbiota that were predicted to be ∼5.4 months older than males not experiencing this source of adversity (β=0.471, p=0.033, Table 2; Table S6D; Figure S6C). However, early life drought and maternal social isolation were linked to young-for-age gut microbiota in males (drought effect: β=-0.451, p=0.021; maternal social isolation effect: β=-0.395, p=0.006, Table 2; Table S6D; Figure S6D). Probably as a result of these conflicting effects, we found no effect of cumulative early adversity on microbiome △age in males (Table S6B).

Microbiome age does not predict the timing of development or survival

Finally, we tested whether variation in microbiome △age predicted the timing of individual maturational milestones or survival using Cox proportional hazards models (Table S8). Maturational milestones for females were the age at which they attained their first adult dominance rank, reached menarche, or gave birth to their first live offspring (Figure 1C). Male maturational milestones were the age at which they attained testicular enlargement, dispersed from their natal social group, or attained their first adult dominance rank (Figure 1D). We also tested if microbiome △age predicted juvenile survival (in females and males) or adult survival (females only). We did not test adult survival in males because male dispersal makes it difficult to know age at death for most males (83).

Contrary to our expectations, microbiome △age did not predict the timing of any baboon developmental milestone or measure of survival (Tables S9, S10). However, these patterns should be treated with caution, as reflected by the large number of censored animals, large hazard ratios, and small sample sizes for some tests.

Discussion

We report three main findings. First, similar to humans and other animals (11–24), baboon gut microbiota show consistent age-related changes in taxonomic composition that produce a dependable microbiome-based age predictor—a microbiome clock. This clock explains nearly half the variance in true host chronological age, and variation in its age predictions recapitulate well-known patterns of faster male senescence in humans and other primates (84). Second, parallel to a recent epigenetic clock in the Amboseli baboons (45), deviations from microbiome age predictions are predicted by the current socio-environmental conditions experienced by individual hosts, although the effect sizes are relatively small. Notably, recent social competition as reflected in social dominance rank predicts both microbiome and epigenetic age. Third, microbiome age did not seem to predict the timing of individual development or survival (but caution is warranted, given small sample sizes for some tests). Hence, in our host population, gut microbial age reflects current social and environmental conditions, but not necessarily the pace of development or mortality risk.

Our work extends the findings of prior analyses from the longitudinal microbiome data set on the Amboseli baboons (50–52). For instance, Grieneisen et al. (50) investigated microbiome heritability, finding that the heritability of microbiome features rises with host age, suggesting both an increasing role for host genotype and a decreasing role for host environments in shaping microbiome composition with host age. This result is consistent with our finding that social and environmental determinants of microbiome age, while statistically significant, tend to have small effect sizes. Further, Roche et al. (52) found that members of the same age class have more similar ASV-ASV correlation patterns than members of different age classes—a finding consistent with the observation that the most predictive microbiome age predictor was produced by the Gaussian process regression model, which can more explicitly incorporate relationships between features than other machine learning algorithms.

To date, five other microbiome clocks have been built—all in human subjects—that predict host chronological age (46–49). Compared to these clocks, our clock in baboons has comparable or better predictive power, with a median error of 1.96 years, compared to 6 to 11 years in human age-predicting clocks ((46–49); baboon lifespans are approximately one third of a human lifespan). Age prediction may be more successful in baboons than humans for at least three reasons. First, signs of age in the baboon gut microbiome may be more consistent across hosts, perhaps because of the relatively homogeneity in host environments and lifestyles in baboons compared to humans (51).

Second, the baboon clock relies on dense longitudinal sampling for each host, and because the training data set included at least one sample from each host, the microbiome clock may be better able to address personalized microbiome compositions and dynamics than clocks that rely on cross-sectional data (e.g. 46–48). However, because our training data set was not naïve to information from the host being predicted, this approach could leak information between the training and test set. Third, our Gaussian process regression approach allowed us to account for non-linear and interactive relationships between microbes with age, leveraging a wider variety of age-related signatures in the microbiome than other machine learning approaches (e.g., elastic net regression or random forest regression).

Despite the relative accuracy of the baboon microbiome clock, its ability to predict individual age is lower than for age clocks based on patterns of DNA methylation—both for humans and baboons (40–43, 45). One reason for this difference may be that gut microbiota are highly personalized: each host species, population, and even host individuals within populations have distinctive, characteristic microbiota, which likely limits the utility of our clock beyond our study population (51, 60, 62, 85). To make a more generalizable clock, an important next step is to train the clock on data from many more host populations and incorporate features of the gut microbiome that are broader and more universal across host species and populations. Because microbiome microbial taxa are often host-specific, we suspect that making a more generalizable clock will require incorporating microbiome features that are widely shared across host species. Despite the limitations of taxon-based microbiome clocks, one advantage of microbiome clocks over epigenetic clocks is that the required data can be collected non-invasively from host individuals, which may make them more amenable for longitudinal sampling or for use in non-human animal populations where invasive collection of blood or tissue samples is challenging or impossible (86, 87).

We connected our microbiome clock age predictions to the social and environmental conditions experienced by host individuals (45, 69, 80, 84). Among the conditions we tested, the most consistent findings were connections to individual dominance rank. Microbiome samples from high-ranking males and females both appeared old-for-age. These results are interesting considering a growing body of evidence that finds rank-related differences in immunity and metabolism, including costs of high rank, especially for males (45, 64, 65, 88–90). For instance, in Amboseli, high social status in males is linked to old-for-age epigenetic age estimates (45), differences in immune regulation (65, 88), and, for alpha males, elevated glucocorticoid levels (64). These patterns, together with our evidence that high-ranking males tend to look ‘old-for-age’, are consistent with the idea that high-ranking males pursue a “live fast die young” life history strategy (45).

Interestingly, however, we also found some evidence for old-for-age microbiome age estimates in samples from high-ranking females who do not seem to be “living fast” in the same sense as high-ranking males (indeed alpha females have lower glucocorticoid hormones than other females (91)). This outcome points towards a shared driver of high social status in shaping gut microbiome age in both males and females. One candidate is priority of access to food, which is a benefit experienced by both high-ranked male and female baboons. Prior research in this population suggests that as animals age, their diets become more canalized and less variable (50). Priority of access to food and fewer foraging disruptions may result in a higher quality, more stable diet, which may be reflected in an old-looking gut microbiome However, this explanation is speculative and more work is needed to understand the relationship between rank and microbiome age.

While some social and environmental conditions associated with baboon development or mortality predict microbiome age, our microbiome clock predictions do not themselves predict baboon development or mortality. This finding supports the idea that microbiome age is sensitive to transient social and environmental conditions; However, these patterns do not have long-term consequences for development and mortality. One reason for this may be that the biological drivers of development and mortality are be too diverse to be well reflected in gut microbial communities.

For instance, animals in Amboseli die for many reasons, including interactions with predators and humans, conflict with conspecifics (92), and disease, and the biological predictors of these events in the gut microbiome are likely weaker or more diverse than the biological signals that predict developmental milestones (i.e., sex steroids, growth hormones, metabolic status, and physical condition). Despite this variation, three important next steps will be (i) to test whether microbiome age is correlated with other hallmarks of biological age in this population, (ii) to test whether it is possible to build a microbiome-based predictor of individual lifespan (i.e., remaining lifespan as opposed to years already lived), and (iii) to test the relationships between microbiome compositional features and individual survival. These future directions are important for connecting the microbiome, including individual features of the microbiome, to aging processes, as opposed to a simple measure of chronological age.

In sum, our findings support the hypothesis that the gut microbiome serves as a biomarker of some aspects of host age. By leveraging microbial, social, environmental, and life history data on host individuals followed from birth to death, we bolster the validity of microbiome clock studies in humans and find some socio-environmental predictors of microbiome age. Future work may also benefit from searching for more universal aspects of the microbiome that may predict host aging across populations and even host species.

Materials and methods

Study population and subjects

Our study subjects were 479 wild baboons (215 males and 264 females) living in the Amboseli ecosystem in Kenya between April 2000 to September 2013. The Amboseli baboon population is primarily composed of yellow baboons (Papio cynocephalus) with some admixture from nearby anubis baboon (Papio anubis) populations (93–95). Prior research in our population finds no link between host hybrid ancestry and microbiome composition (85). Since 1971, the Amboseli Baboon Research Project (ABRP) has been collecting continuous observations of the baboons’ demography, behavior, and environment (75). The baboons are individually identified by expert observers who visit and collect data on each social group 3 to 4 times per week (the subjects lived in up to 12 different social groups over the study period). During each monitoring visit, the observers conduct group censuses and record all demographic events, including births, maturation events, and deaths, allowing us to calculate age at maturity and lifespan with precision. This research was approved by the IACUC at Duke University, University of Notre Dame, and Princeton University and the Ethics Council of the Max Planck Society and adhered to all the laws and guidelines of Kenya.

Sample collection, DNA extraction, and 16S data generation

The 13,476 gut microbiome compositional profiles in this analysis represent a subset of 17,277 profiles, which were previously described in (50, 51). The 13,476 samples in our analyses include those from baboons whose birthdates, and hence individual ages, were known with just a few days’ error. Each baboon had on average 33 samples collected across 6 years of their life (Figure 1A and 1B; range = 3 to 135 samples per baboon; median days between samples = 44 days).

Samples were collected within 15 minutes of defecation, homogenized, and preserved in 95% ethanol. Samples were freeze-dried and sifted to remove plant matter prior to long term storage at -80°C. DNA from 0.05 g of fecal powder was manually extracted using the MoBio (Catalog No. 12955-12) and QIAGEN (Catalog No. 12955-4) PowerSoil HTP kits for 96-well plates using a modified version of the MoBio PowerSoil-HTP kit. Specifically, we followed the manufacturers’ instructions but increased the amount of PowerBead solution to 950 μL/well and incubated the plates at 60°C for 10 minutes after the addition of PowerBead solution and lysis buffer C1.

Following DNA extraction, a ∼390 bp segment of the V4 region of the 16S rRNA gene was amplified and libraries prepared following standard protocols from the Earth Microbiome Project (96). Libraries were sequenced on the Illumina HiSeq 2500 using the Rapid Run mode (2 lanes per run). Sequences were single indexed on the forward primer and 12 bp Golay barcoded. The resulting sequencing reads were processed following a DADA2 pipeline (97), with the following additional quality filters: we removed samples with low DNA extraction concentrations (< 4x the plate’s blank DNA extraction concentration), samples with <1000 reads, and amplicon sequence variants that appeared in one sample (see (50) for details). ASVs were assigned to microbial taxa using the IdTaxa(…) function in the DECIPHER package, against the Silva reference database SILVA_SSU_r132_March2018.RData (98, 99). The final set of samples had 1,017 to 427,454 reads (median =51,839 reads), with 8,492 total ASVs.

Identifying microbiome features that contribute to age predictions and that change with age

To identify microbiome features that change with host age, we ran linear mixed models on 1,440 microbiome features (Table S1). Models were run using the R package lme4, with p-value estimates from lmerTest (100, 101). These features included: (i) five metrics of alpha diversity; (ii) the top 10 principal components of microbiome compositional variation; (iii) centered log ratio transformed abundances of each microbial phyla (n = 30), family (n = 290), genus (n = 747), and amplicon sequence variance (ASV) detected in >25% of samples (n=358). Alpha diversity metrics were calculated using the R package vegan and principal components of microbiome compositional variation were calculated using the R package labdsv (102, 103).

For each feature, we modeled its relationship to host chronological age using both linear and quadratic terms. To make our quadratic terms more interpretable, we centered our age estimates on zero by subtracting the average age in the dataset from each age value. Specifically, when a quadratic term is negative, the curve is concave, whereas when the term is positive, the curve is convex. We also included season (wet or dry) and z-scored rainfall and temperature as fixed effects, and individual identity, social group at time of collection, hydrological year, and the DNA extraction/PCR plate identity were modeled as random effects. All community features (i.e., alpha diversity and principal components) and all taxa present in >25% of samples were modeled using a Gaussian error distribution. We extracted the coefficient, standard error, and p-value for the age term, then corrected for multiple testing using the false discovery rate approach of Benjamini and Hochberg (104).

Building the gut microbiome clock

We created a microbiome clock by fitting a Gaussian process (GP) regression model (with a kernel customized to account for heteroskedasticity) to predict each baboon’s chronological age at the time of sample collection using 9,575 microbiome compositional and taxonomic features present in at least 3 samples (Table S2; i.e., we did not restrict the features in the clock to the 1,440 most abundant features used in the age-association analyses described above). The GP regression model with heteroskedasticity correction was the best performing of four supervised machine learning approaches we considered, including elastic net, random forest, and Gaussian process regression with and without the heteroskedasticity kernel (Figure S3; Table S11; See Supplemental Methods for a comparison of other algorithms). Pearson’s correlations between age predictions across the four methods ranged from 0.69 between the random forests and the GP regression model without the heteroskedasticity kernel to 0.96 between the two GP regression models (with and without the heteroskedasticity kernel; Table S11).

Gaussian process regressions were conducted in Python 3 using scikit-learn (105, 106). As a nonparametric, Bayesian approach that infers a probability distribution over all the potential functions that fit the data, the Gaussian process regression does not assume a linear relationship between chronological and predicted age (107). For the prior distribution in the Gaussian process regression, we used a radial basis function as our kernel and set the scale parameter to the mean Euclidean distance of the dataset, as calculated in the R package vegan (102). Because initial, exploratory models exhibited heteroskedasticity (Figure S7), we multiplied the variance in the training data by the radial basis function, which distributed the higher variance in later life more evenly across lifespan.

To calculate a microbial age estimate for every sample, and to estimate generalization error, we used nested five-fold cross validation. In each of the five model runs, we used 80% of the data to train the model, and the remaining 20% of the dataset as the test data. Because host identity has a strong effect on microbiome composition in our population (51), we distributed samples from each host across the five test/training data sets by randomly assigning each sample a test set without replacement. Hence, the training data set was not naïve to information from the given host being predicted as some samples for that host were included in the training set. For each model run, four of the test datasets were treated altogether as training data and the 5th set was the validation test set. We then took the estimates from all five model runs and estimated global model accuracy on the aggregated estimates.

We assessed the accuracy of our microbiome clock by regressing each sample’s chronological age (agec) against the model’s predicted microbial age (agem) and determining the R² value and Pearson’s correlation between agec and agem. We also calculated the median error of the model fit as the median absolute difference between agec and agem across all samples (40).

Calculating microbiome △age estimates

To characterize patterns of microbiome age from our microbiome clock, we calculated sample-specific microbiome △age as the difference between a sample’s microbial age estimates, agem, and the host’s chronological age at the time of sample collection, agec. Higher microbiome △ages indicate old-for-age microbiomes, as agem > agec, and lower values (which are often negative) indicate a young-for-age microbiome, where agec > agem. Because the microbiome clock systematically over predicted the ages of young animals and under predicted the ages of old animals), we also calculated a “corrected microbiome △age” as the residuals of agem correcting for host chronological age, season, monthly temperature, monthly rainfall, and social group and hydrological year at the time of collection. This measure is used for visualizations of the predictors of microbiome age and for testing whether average microbiome △age predicts developmental milestones or survival.

Testing sources of variation in microbiome △age

Many social and environmental factors have been shown to predict fertility and survival in the Amboseli baboons (54, 59, 67, 68, 77). To test if some of the most important known factors also predict patterns of microbiome age, we used linear mixed models to test predictors of microbiome age in individual samples separately for males and females.

In these models, the response variable was the sample-specific measure of △age (agem - agec). All models included the following fixed effects: individual chronological age at the time of sample collection, to correct for model compression; the average maximum temperature during the 30 days before the sample was collected, total rainfall during the 30 days before the sample was collected, and the season (wet or dry) during sample collection. Every model also included, as fixed effects, measures of early life adversity the individual experienced prior to 4 years of age (Table S7). These were modeled as either as six, individual, binary variables, reflecting the presence or absence of each source of adversity in the first four years of life, or as a cumulative sum of the number of sources of adversity the individual experienced, also in the first four years of life (67). Social rank at the time of sampling was also modeled as a fixed effect. For males we used ordinal rank, and for females we used proportional rank (108). To make model interpretation more intuitive (high rank corresponds to higher values), we multiplied the coefficients for ordinal rank and maternal rank by -1. Random effects included individual identity, the social group the individual lived in at the time of collection, and hydrological year. In models of microbiome age in females, the number of adult females in the group at the time of sample collection was included as female-specific measure of resource competition.

Testing whether microbiome △age predicts baboon maturation and survival

We used Cox proportional hazards models to test whether microbiome △age predicted the age at which females and males attained maturational milestones and the age at death for juveniles and adult females (Table S8). We only measured adult survival in females because males disperse between social groups, often repeatedly across adulthood, making it is difficult to know if male disappearances are due to dispersal or death (83). For females, the maturational milestones of interest were the age at adult rank attainment (median age 2.24 in Amboseli), age at menarche (median age 4.51 in Amboseli), and the age at first live birth (median age 5.97 in Amboseli). For males, these milestones were the age of testicular enlargement (median age 5.38 in Amboseli), the age of dispersal from natal group (median age 7.47 in Amboseli), and the age of first adult rank attainment (i.e., when a male first outranks another adult male in his group’s dominance hierarchy; median age 7.38 in Amboseli) (55, 56). See full descriptions of each milestone in Table S8. To be included in these analyses, animals must have reached the milestone after the onset of sampling (April 2000) and had at least three samples available in the timeframe of interest. We verified that none of our models violated the proportional hazards assumption of a Cox regression.

The variables we modeled differed based on the event of interest. However, all models included as fixed effects corrected △age as the residuals of gut microbiome △age averaged over the timeframe. All models of developmental milestones also included variables tested in Charpentier et al. 2008 (55, 56): (iii) maternal presence at the time of the milestone, (iv) the number of maternal sisters in the social group, averaged over the timeframe, (v) rainfall averaged over the timeframe, and (vi) whether the subject’s mother was low ranked (was in the lowest quartile for female ordinal rank). For female-specific milestones, we also included (vii) the average number of adult females in the group averaged over the timeframe, and for male-specific milestones we included the number of excess cycling females in the group averaged over the timeframe, or the difference between the number of cycling females and the number of mature males within a subject’s social group. Last, we included (viii) the subject’s hybrid score, which is an estimation of the proportion of an individual’s genetic ancestry attributable to anubis or yellow baboon ancestry (95).

All juvenile survival models included as fixed effects the residuals of microbiome △age averaged over the timeframe and measures the cumulative number of sources of early life adversity each individual experienced (67). Additionally, we ran three versions of the juvenile survival analysis: two subset to each sex, and one version that included both sexes. In the model including both sexes, we included sex as a predictor.

Adult female survival models included the same variables as for juvenile survival, but additionally included average lifetime dyadic social connectedness to adult females, average lifetime dyadic social connectedness to adult males, and average lifetime proportional rank. Full descriptions of all predictors are available in Table S12.

Acknowledgements

We thank Jeanne Altmann for her essential role in stewarding the Amboseli Baboon Project, and in collecting many of the fecal samples used in this manuscript. In Kenya, we thank the Kenya Wildlife Service, the Wildlife Training and Research Institute, the National Council for Science, Technology, and Innovation, and the National Environment Management Authority for permission to conduct research and collect biological samples. We also thank the University of Nairobi, Institute of Primate Research, National Museums of Kenya, the Amboseli-Longido pastoralist communities, the Enduimet Wildlife Management Area, Ker & Downey Safaris, Air Kenya, and Safarilink for their cooperation and assistance in the field. We thank Karl Pinc for managing and designing the database. We thank Raphael Mututua, Serah Sayialel, Kinyua Warutere, and Long’ida Siodi for collecting the field data and fecal samples. We also thank Tawni Voyles, Anne Dumaine, Yingying Zhang, Meghana Rao, Tauras Vilgalys, Amanda Lea, Noah Snyder-Mackler, Paul Durst, Jay Zussman, Garrett Chavez, and Reena Debray for contributing to fecal sample processing. Complete acknowledgments for the ABRP can be found online at https://amboselibaboons.nd.edu/acknowledgements/.

Additional information

Funding

This work was supported by the National Institutes of Health and National Science Foundation, especially the National Institute on Aging for R01 AG071684 (EAA), R21 AG055777 (EAA, RB), NIH R01 AG053330 (EAA), NIH R35 GM128716 (RB), NSF DBI 2109624 (MRD), and NSF DEB 1840223 (EAA, JAG). We also thank the Duke University Population Research Institute P2C-HD065563 (pilot to JT), the University of Notre Dame’s Eck Institute for Global Health (EAA), and the Notre Dame Environmental Change Initiative (EAA). We also thank Duke University, Princeton University, the University of Notre Dame, the Chicago Zoological Society, the Max Planck Institute for Demographic Research, the L.S.B. Leakey Foundation and the National Geographic Society for support at various times over the years.

Author contributions

Conceptualization: MRD, JT, EAA

Methodology: MRD, KER, DJ, JA, SM, JT, EAA

Supervision: SCA, LBB, JAG, RB, SM, JT, EAA

Writing—original draft: MRD, EAA

Writing—review & editing: MRD, KER, DJ, JA, SCA, LBB, JAG, RB, SM, JT, EAA

Competing interests

The authors declare that they have no competing interests.

Data and materials availability

All data is available on Dryad at https://doi.org/10.5061/dryad.b2rbnzspv. Code is available at the following GitHub repository: https://github.com/maunadasari/Dasari_etal-GutMicrobiomeAge

Supplementary Text

Creating and assessing age-predictive machine learning models

Introduction to the approaches

To create our final microbiome aging clock, we tested three supervised machine learning algorithms: elastic net regression, Random Forest regression, and Gaussian process regression (107, 109, 110). Below we summarize the strengths and weaknesses of each machine learning algorithm.

Elastic net regression is a regression algorithm that produces a linear model. It improves upon the predictions from simple linear regressions by incorporating coefficient penalties from the L1 regularization (LASSO regression) and L2 regularization (ridge regression) (109). Elastic net regression is infrequently used in microbiome studies but has produced promising results in epigenetic aging clocks due to its flexibility in choosing which features to keep and which to remove (40–43, 45). However, elastic net regressions produce linear relationships between the input chronological age and the predicted age, which may not accurately affect the true relationship between chronological age and the microbiome.

Random Forest regression is an ensemble learning method that creates a number of parallel decision trees, each producing its own prediction (110). The prediction is then averaged among all trees to create the final estimate. A key advantage of Random Forest regression over elastic net regression is that it does not assume a linear relationship between the predicted estimate and the input chronological age, but the model may be biased by correlated features. Random Forest regression is commonly used in microbiome research, including other microbiome clocks (36, 111–113).

Gaussian process regression is a nonparametric, Bayesian approach that infers a probability distribution over all the potential functions that fit the data (107). Like random forest, Gaussian process regressions do not assume a linear relationship between chronological age and predicted age but has the additional advantage of kernel customization. As such, Gaussian process regressions may be able to better handle heteroskedasticity in the data (an issue in our clock; see below). As an increase in chronological age is often associated with the breakdown of physiological processes (e.g. aging), heteroskedasticity in microbial age estimates may indicate a breakdown of the host’s processes that regulate the gut microbiome.

Methods and optimization of machine learning algorithms

Prior to running each algorithm, all features were center log ratio transformed within sample. We then chose a ratio of training to test dataset. To do this, we first compared the model fit of different ratios of training to test sets. These included the following training:test splits: 50:50, 60:40, 75:25, 80:20, and 90:10. We found that an 80:20 data split provided the best balance between model performance and the risk of overfitting.

In order to calculate a microbial age estimate for every sample and estimate generalization error, we used a nested cross-validation framework. Each of the three algorithms has its own internal cross-validation where a subset of the training data is held apart and used to internally validate the model. We added an additional, external layer of cross-validation with our 80:20 training:test data split. We classified samples into five different test sets where individual was as evenly represented as possible in all training and test sets. As the number of samples varied between individuals, we randomly assigned each sample a test set without replacement if an individual’s sample count was less than five, or with replacement if an individual’s sample count was greater than five. For each model run, four of the test datasets were treated altogether as training data and the fifth set was the validation test set.

Elastic net regressions were run in R using function cv.glmnet() from package glmnet (114). The two main parameters for this model are λ, which is the penalty from the LASSO regression that penalizes extra predictors by shrinking coefficients to zero, and α, the parameter that balances between minimizing between the residual sum of squares and minimizing the magnitude of the coefficients. cv.glmnet() automatically fits 100 values of λ by default and names the λ that produces the minimum cross-validated error “lambda.min”. We used lambda.min as our value of λ. For α, we manually ran the model with 200 values of alpha (from 0 to 1 in increasing increments of 0.005) and picked a value of alpha that would minimize the mean absolute error and maximize the adjusted R².

Random Forest regressions were conducted in Python 3 using scikit-learn (105). The main parameter was the number of decision trees being used, which defaults to 100. Too many trees could result in overfitting so in order to minimize overfitting and optimize R², we ran a series of Random Forest regressions with different numbers of trees: we increased the number of trees in increments of 50, stopping at 400 because of minimal changes in R² relative to 200 trees.

Gaussian process regressions were also conducted in Python 3 using scikit-learn (105, 106). In both the non-heteroskedastic-kernel model and heteroskedastic-kernel model, the main parameters we used to modify the kernel function included the scale and bounds. These parameters moderate the level of overfitting in the algorithm: the scale parameter specifies a starting point for which the algorithm optimizes within the confines of the bounds parameters. As with the other models, we incrementally changed both the scale parameter within a wide range of bounds and checked the output model’s R² and median error. Our final model retained a wide range of bounds (1 to 100) and set the scale parameter to the median euclidian distance of the dataset as calculated in R using function vegdist() from R package vegan (102).

Due to the heteroskedasticity exhibited by the models above (Figure S7), we modified the Gaussian process regression’s kernel function further to account for the variance within the dataset. Specifically, we multiplied the variance in the training data by the radial basis function, which distributed the higher variance in later life more evenly across lifespan.

Comparison of machine learning algorithms

To assess model accuracy, we used the predicted age estimates from all 5 runs of the nested cross-validation procedure to assess model fit and accuracy. As in Horvath (2013), we regressed the sample’s predicted microbial age (agem) against the host’s known chronological age (agec) and calculated: (1) the R² between agec and agem; (2) the Pearson’s correlation coefficient between agec and agem; and (3) the median error as the median absolute difference between agec and agem (Table S13 and Figure S3). Across all algorithms, we observed that males always aged faster than females, which is consistent with well-known patterns of sex-specific senescence in humans and other primates (84) (Figure S3). The Gaussian process regression with the heteroskedastic kernel was the best model for every metric assessed - it maximized R² and Pearson’s R to 0.488 and 0.698 (respectively) while minimizing median error. It also was the only model with which we were able to alleviate any heteroskedasticity.

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.