To investigate the signature of environmental variation on the baboon DNA methylome, we used reduced-representation bisulfite sequencing [RRBS (61, 62)] to profile DNA methylation in blood for 477,270 CpG sites in the baboon genome, in 256 unique individuals (115 males, 141 females; SI Appendix, Fig. S3). This set of sites is a subset of the 1,590,767 CpG sites captured in the full sequencing dataset, filtered to remove invariant, constitutively hypomethylated, and constitutively hypermethylated sites (SI Appendix). For 37 individuals, we profiled repeated, longitudinally collected samples (2-3 samples per individual), for a total of n = 295 samples (SI Appendix, Table S1).

For each CpG site separately, we first modeled DNA methylation levels as a function of habitat quality at birth, cumulative early-life adversity, and age and ordinal dominance rank at the time of sampling, using the binomial mixed effects model implemented in MACAU (we refer to this analysis as Model 1; Fig. 1A; see SI Appendix for model details) (63). We quantified habitat quality at birth as a simple binary variable indicating whether each study subject was born before or after the home range shift described above (N = 57 individuals were born in the low-quality habitat). We treated habitat quality at birth separately from cumulative early adversity because of its nature as a strong cohort effect characterized by two distinct time periods, rather than a set of conditions that vary across individuals living at the same time and place (see SI Appendix for a discussion of the resulting unavoidable correlation with time). We considered five sources of early adversity as components of the cumulative early adversity measure: drought, maternal loss, large group size, the presence of a close-in-age younger sibling, and low maternal dominance rank, which collectively predict both reduced survival and reduced offspring survival in this population (5, 45) (see also Materials and Methods). We estimated dominance rank effects for each sex separately (by nesting rank within sex), as male and female ranks depend on different traits for each sex (i.e., kinship in females and physical condition in males). Further, the hierarchies for each sex are separately estimated, have sex-specific implications, and have sex-specific associations with gene expression (44, 57, 64–67).

In Model 1, the strongest predictors of DNA methylation in adulthood were habitat quality at birth, male dominance rank at sample collection, and age at sample collection. The relationship between habitat quality at birth and DNA methylation was striking, resulting in 3,296 habitat quality-associated sites (10% FDR; SI Appendix, Table S2A). Consistent with the association between dominance rank and other aspects of gene regulation (57, 67), associations between male dominance rank and DNA methylation were also widespread (n = 3,736 sites, 10% FDR), in contrast to a weaker relationship with female dominance rank (57, 67–69) (n = 4 sites; see SI Appendix for a discussion of this sex difference, which is consistent with previous findings in our population). Age strongly predicted DNA methylation across the genome (n = 169,439 age-associated sites), with a bias, as reported in other studies (68, 70), to increases in DNA methylation with age in CpG islands (65%) and decreases in DNA methylation with age in most other regions of the genome (79%). In contrast to these three effects, we observed no significant associations (10% FDR) between DNA methylation and cumulative early adversity.

Our results for Model 1 suggest that habitat quality in early life is particularly important in the lives of baboons and could moderate the association between other sources of early adversity and DNA methylation. To test this possibility, we re-ran our analyses, but in this case tested for the effects of cumulative early adversity experienced in the high-quality habitat and low-quality habitat separately (i.e., by nesting cumulative early adversity within habitat quality; Model 2). To maximize power, we also included individuals for whom early adversity data were available, but dominance rank data were missing because of observational gaps for males. This model not only strengthens the evidence for a main effect of habitat quality (25,509 habitat quality-associated sites; 10% FDR) but reveals an interaction with cumulative adversity: 2,856 sites are associated with cumulative adversity for baboons born in low-quality habitat (10% FDR), while none are significantly associated with cumulative adversity in baboons born in high-quality habitat (Fig. 1 B and C and SI Appendix, Table S2B). Notably, only 64 of 295 samples derive from low-quality habitat individuals, suggesting that the greater number of associations in low-quality habitat is not driven by greater power. Among the significant sites identified in samples from individuals born in low-quality habitat, the effect sizes for cumulative adversity in the low-quality habitat are uncorrelated with the effect sizes for cumulative adversity in high-quality habitat (p = 0.838) but positively correlated with the effect sizes for habitat quality itself (R = 0.508, P < 1 × 10⁻¹⁰). This result suggests that the signature of cumulative adversity is amplified by exposure to ecologically challenging conditions (and vice versa). Importantly, cumulative adversity scores do not differ between animals born in low-quality and high-quality habitats (Wilcoxon rank-sum test P = 0.843).

To investigate whether different components of the cumulative adversity score contribute differently to the early adversity–DNA methylation relationship, we then ran a third model (Model 3) to evaluate each of the five individual sources of early adversity, nested within habitat quality (all other biological and technical covariates remained the same as in Model 2). Among the individual sources of adversity we considered, early-life drought most clearly predicted variation in DNA methylation across the genome, especially for individuals born in the low-quality habitat (25,355 sites at a 10% FDR; Fig. 2 A and B; note that drought also affected 11% of baboons born in the high-quality habitat). We also identified detectable, but less common signatures of maternal loss (4,893 sites), large group size (3,124 sites), low maternal rank (730 sites), and the presence of a close-in-age younger sibling (619 sites). In contrast, none of the individual sources of early adversity were robust predictors of DNA methylation for individuals born in the high-quality habitat (≤5 sites associated with any individual predictor at 10% FDR; SI Appendix, Table S2C; see SI Appendix, Fig. S4 for comparisons of effect sizes for each predictor in low- versus high-quality habitat).

Our models indicate that some early-life experiences are linked to more pronounced DNA methylation signatures than others. Drought in particular, which is one of the least predictable environmental exposures in Amboseli, is associated with an order of magnitude more CpG sites than maternal rank or group size, the next most common effects. Notably, early-life rainfall and rainfall at the time of sampling are only weakly correlated in our dataset, supporting the idea that our observations capture the signature of early life, not experience at the time of biological sample collection (SI Appendix, Fig. S5). To investigate whether these signatures are unique to specific early-life experiences or reflect a general signature of stress and adversity (perhaps scaled to the magnitude of the stressor), we therefore tested for overlap between the sets of sites linked to each of the five individual-level predictors and to habitat quality based on results from Model 3.

Our results support a generalized rather than an exposure-specific signature (Fig. 2C). Specifically, among sources of early adversity with a substantial number of associated CpG sites (habitat quality, drought, maternal loss, and group size), sites associated with one early-life exposure are 1.04 – 8.6-fold more likely to be associated with a second early-life exposure (P < 1 × 10⁻¹⁰ for 4 of 6 comparisons). Habitat quality and drought (in samples from individuals born in low-quality habitat) show a particularly striking pattern of overlap: 4,038 CpG sites are significantly associated with both predictors (log₂(OR)=2.23, P < 1 × 10⁻¹⁰), and almost all of these effects (99.8%) are directionally concordant, such that exposure to low habitat quality in early life and exposure to drought predict the same direction of effect. Drought-associated sites were enriched in a broad group of gene sets annotated in the Molecular Signatures Database Hallmark set, including pathways involved in development (e.g., genes expressed in pancreatic B cells:), cell proliferation (e.g., targets of MYC signaling), and cellular metabolism (e.g., genes involved in oxidative phosphorylation; all Bonferroni-corrected P < 0.05; no gene sets pass this threshold for habitat quality-associated sites; see SI Appendix, Table S7).

Comparing these findings to the signature of male dominance rank shows that overlap in sensitivity to the environment is not specific to early-life variables [note that we focused on male rank here for comparison because significant associations with female rank are far less common, consistent with findings from our previous work on rank and gene expression in the Amboseli baboons: (57, 67)]. Male rank-associated sites are 11.21 times more likely to be associated with drought than background expectations and 2.43 times more likely to be associated with habitat quality (both P < 1 × 10⁻¹⁰). In these cases, dominance rank effects tend to have directionally opposite effects to habitat quality and drought (log₂(OR) = −4.06 for overlap with habitat quality; the odds ratio could not be estimated for the overlap with drought because there was no overlap in the direction of effects). Consequently, sites that are more highly methylated in high-ranking males also tend to be more highly methylated for baboons of both sexes who were born in poor-quality habitat and exposed to drought within that habitat.

In contrast to male rank-associated patterns of DNA methylation, age effects only modestly overlap with drought effects and habitat quality (log₂(OR) = 0.16 and 0.43, both P < 10⁻¹⁰) and do not overlap with male rank effects at all (log₂(OR) = 0.045, P = 0.35) (SI Appendix, Fig. S6). These results suggest that despite a shared epigenetic signature of at least some types of early and adult experience (with variation in the magnitude of the effect), the signature of age is distinct and the distribution of differentially methylated sites across the genome is not inevitable for any methylation-associated variable. To test this hypothesis further, we investigated how CpG sites related to age versus socioenvironmental variables are distributed across promoters, gene bodies, CpG islands and shores, putative enhancer elements, and unannotated regions. We focused on the four variables with the strongest DNA methylation signatures: age, habitat quality in early life, drought (in the low-quality habitat), and male dominance rank.

Our results highlight two patterns (Fig. 3 A and B). First, compared to the full reference set of sites we analyzed, drought and male dominance rank-associated sites are systematically enriched in functionally important regions of the genome, especially gene bodies (log₂(OR) = 0.25 and 0.72, respectively) and putative enhancer elements (log₂(OR) = 0.52 and 0.99), but depleted in unannotated regions (log₂(OR) = −0.13 and −0.36) of the genome (all P < 1 × 10⁻⁷; SI Appendix, Table S3). Second, and in contrast, age-associated sites are 1.57-fold more likely to occur in unannotated regions of the genome than expected by chance (i.e., in comparison to the background set of sites included in our study), but are depleted in enhancers (log₂(OR) = −0.27) and gene bodies (log₂(OR) = −0.28, all P < 1 × 10⁻¹⁰). Notably, habitat quality-associated sites, which are much more widely distributed in the genome than rank- or drought-associated sites, follow an intermediate pattern: They are less common in unannotated regions than age-associated sites but are not clearly enriched for gene bodies or enhancers.

To investigate the difference in the functional significance of age- versus socioenvironmental variable-associated sites, we turned to more granular chromatin state annotations available from the Roadmap Epigenomics Consortium, which distinguishes between a class of states associated with “active” versus “repressed” (inactive or quiescent) DNA. We observe clear differences in the distributions of age versus socioenvironmental differentially methylated sites based on enrichment within chromatin state annotations (i.e., predictions of the regulatory state of different regions of the genome based on the presence of epigenetic modifications and chromatin accessibility) (71). Specifically, we lifted over chromatin state coordinates for human peripheral blood mononuclear cells to the baboon genome, Panubis1.0 (27, 72). Relative to the background set of sites assayed in our study, early-life drought and dominance rank are again enriched in regions of the genome marked for regulatory activity, such as enhancer elements (log₂(OR) = 0.91 and 1.75, respectively, both P < 10⁻¹⁰), and transcriptional activity (log₂(OR) = 0.62 and 1.52 respectively, both P < 10⁻¹⁰), but depleted in repressed and silenced regions such as heterochromatin (log₂(OR) = −0.89 and −1.32, P = 4.4 × 10⁻⁵ and 0.055) and weakly repressed, polycomb-marked DNA (log₂(OR) = −0.13 and −0.34, P = 0.03 and 0.016; Fig. 3 C, Top; SI Appendix, Table S3). Age-associated sites show the opposite pattern (Fig. 3 C, Bottom).

Although the individuals in our dataset were predominantly adults, individuals exposed to poor habitat quality were sampled at a range of ages (range = 2.5 to 26.3 y). We took advantage of this variation to test whether the signature of early-life adversity attenuates over time, resulting in weaker signatures with longer times from exposure. To do so, we focused on habitat quality, the strongest association with early life we observed in our data. We first built an elastic net model to ask whether early-life exposure to low-quality habitat (a binary variable indicating whether the subject was born before or after the home range shifts) is predictable based on DNA methylation levels sampled in adulthood (73).

We found that an elastic net model achieves high accuracy in our sample (AUC = 0.92 based on leave-one-out cross-validation; Fig. 4 A and B). However, among animals born in low-quality habitats, the ability of the model to correctly and confidently predict habitat quality in early life depends on the time elapsed between the habitat shift and blood sample collection (linear model P = 0.0084; Fig. 4C), but not on the cumulative amount of time spent in low-quality habitat (linear model P = 0.279). Animal age does not predict overall habitat quality, so we are not indirectly capturing habitat quality through the age of individuals in our sample (linear model P = 0.80). Specifically, animals who had spent more time in high-quality habitat prior to sampling were less confidently predicted to be born in low-quality habitat than those who experienced it more recently. This result suggests that, although DNA methylation signatures of early adversity can persist for years in baboons, they also decay over time or are overwritten by later life experience.

The distribution of environment-associated CpG sites in loci related to transcription and active gene regulation suggests that some subset of these sites have the capacity to causally influence gene expression. To formally test this hypothesis, we performed a massively parallel reporter assay (MPRA), mSTARR-seq, designed to both identify loci capable of regulatory activity in vitro and quantify the effects of differential methylation on the magnitude of this activity (Fig. 5A and SI Appendix, Table S4) (29). mSTARR-seq tests a sequence fragment’s ability to drive gene expression in a self-transcribing plasmid, in hundreds of thousands of genomic fragments simultaneously. Because active test regions drive transcription by looping to interact with their own promoter, fragments capable of driving their own transcription are considered to have enhancer-like activity in vitro (74). Since the plasmid backbone is devoid of CpG sites, inserted fragments containing CpG sites in their sequence can be tested in either a fully CpG methylated or fully unmethylated state to investigate whether enhancer activity can be modified by changes in DNA methylation alone.

We performed mSTARR-seq using a mechanically fragmented and restriction enzyme-digested library of baboon DNA fragments, transfected into the human K562 cell line (SI Appendix, Table S4). K562s are a myelogenous leukemia line that shares properties with several types of peripheral blood mononuclear cells and are therefore often used in studies of immune variation. Importantly, mSTARR-seq has also been extensively optimized in K562 cells (29). Following quality control, we were able to test for regulatory activity in 252,463 500-base pair windows across the baboon genome (4.4% of the genome), of which 32,634 contained tested CpG sites in the Amboseli baboon dataset. Among these 32,634 windows, we identified 492 windows (1.5% of those tested, using a 10% FDR threshold; SI Appendix, Table S5) capable of enhancer activity in either an unmethylated state, a fully methylated state, or both (similar to estimates from ref. 29).

As expected, experimentally identified regulatory regions were strongly enriched in predicted strong enhancers (based on chromHMM annotations: log₂(OR) = 2.50. P < 10⁻¹⁰) and depleted in insulators (log₂(OR) = −2.32; P = 2 × 10⁻⁹) and repressed regions (log₂(OR) = −0.94, P = 5.8 × 10⁻¹⁸; Fig. 5B and SI Appendix, Table S6). Among the 492 regulatory windows overlapping our tested sites, 86% also exhibited methylation-dependent activity, where the capacity to drive transcription differs depending upon whether CpG sites are methylated or not. Of note, 94% of methylation-dependent regions exhibited reduced activity in the methylated state compared to the unmethylated state, consistent with a general role for DNA methylation in repressing gene regulation.

For most discovery thresholds between 5% and 20% FDR, male dominance rank-associated CpG sites are found in mSTARR-seq-identified regulatory windows more often than chance, such that most discovery threshold-predictor variable combinations reach statistical significance (Fig. 5C). For example, 40 dominance rank-associated sites (FDR = 20%; 2.2% of significant sites tested) fall in regions of the genome capable of behaving like enhancer elements [log₂(OR) = 0.59 FET, P = 0.014; note that promoter regions often also exhibit enhancer-like activity in massively parallel reporter assays (29)]. These sites are also enriched (relative to all sites we analyzed within mSTARR-seq analysis windows) in windows where modifying the DNA methylation level of the tested sequence alters its capacity to drive gene expression (log₂(OR) = 0.71 FET, P = 0.006). The pattern for drought-associated sites is less clear: while they are not more likely to occur in mSTARR-seq regulatory elements than expected by chance (drought-associated sites at FDR = 20%, n = 196 sites, (log₂(OR)=0.18 FET P = 0.1), they exhibit a modest enrichment for DNA methylation-dependent activity (log₂(OR) = 0.23, FET P = 0.06), and 88% of early drought-associated sites that overlap with mSTARR-seq regulatory windows are in methylation-dependent windows. Age-associated sites again provide a clear contrast, with no evidence of enrichment of such sites among mSTARR-seq-identified regulatory windows (log₂(OR) = 0.065, P = 0.45).

These results suggest that the associations detected in our field-based sample partially reflect targeted, functionally important changes in the response to the environment, some of which are detectable years to decades post-exposure. If so, environmental effects on DNA methylation should also colocalize with environmental effects on gene expression in the Amboseli baboons. To test this possibility, we drew on RNA-seq gene expression data from white blood cells collected from 2013 to 2018, in which several thousand associations between male rank and gene expression have previously been identified (57, 67) (note that individuals born in low-quality habitat are not well-represented in this dataset because the start of collection for gene expression analysis long post-dated the habitat shift, so we were not able to perform a parallel analysis for habitat quality or early-life drought). Male rank-associated CpG sites fall closer to, and more often within, genes associated with male dominance rank than they do for the background set of tested genes (Kolmogorov–Smirnov test P = 1.81 × 10⁻⁵). Blood-expressed genes that contain a male rank-associated CpG site are also 1.22-fold more likely to exhibit male rank-associated gene expression levels (P = 6.60 × 10⁻⁴), even though the individuals represented in the gene expression dataset and the DNA methylation dataset are largely distinct (34 of 115 males in the DNA methylation dataset were included in the gene expression dataset; 34 of 52 males in the gene expression dataset were included in the DNA methylation dataset). Finally, rank effects on gene expression are negatively correlated with rank effects on DNA methylation for CpG sites in the same gene (SI Appendix, Fig. S7). Thus, if DNA methylation levels are higher in high-ranking males, gene expression levels tend to be lower in high-ranking males, and vice versa (FET for sign: log₂(OR) = −1.25, P = 8.6 × 10⁻⁹). As a result, multiple pathways enriched among rank-associated genes based on gene expression [in comparison to the background set of expressed and analyzed genes (57)] are also enriched among genes linked to rank-associated DNA methylation patterns, including interferon alpha signaling, NFkB signaling, and the inflammatory response (all P < 0.05; SI Appendix, Table S7).

J.A.A., A.J.L., E.A.A., S.C.A., and J.T. designed research; J.A.A., D.L., A.J.L., R.A.J., T.V., and M.Y.A. performed research; J.A.A., D.L., A.J.L., and R.A.J. analyzed data; and J.A.A., D.L., A.J.L., R.A.J., T.V., M.Y.A., E.A.A., S.C.A., and J.T. wrote the paper.

The authors declare no competing interest.

Raw FASTA files data have been deposited in NCBI Short Read Archive (PRJNA970398 and PRJNA871297) (88, 89). All statistical analyses in this work were performed in R (Version #4.1.2) (90), with code available at https://github.com/janderson94/Anderson_et_al_socioecological_methylation_predictors (91). Previously published data were used for this work, and are available via NCBI as referenced in SI Appendix, Table S1.