Early-life social experience affects offspring DNA methylation and later life stress phenotype

Studies in rodents and captive primates suggest that the early-life social environment affects future phenotype, potentially through alterations to DNA methylation. Little is known of these associations in wild animals. In a wild population of spotted hyenas, we test the hypothesis that maternal care during the first year of life and social connectedness during two periods of early development leads to differences in DNA methylation and fecal glucocorticoid metabolites (fGCMs) later in life. Here we report that although maternal care and social connectedness during the den-dependent life stage are not associated with fGCMs, greater social connectedness during the subadult den-independent life stage is associated with lower adult fGCMs. Additionally, more maternal care and social connectedness after den independence correspond with higher global (%CCGG) DNA methylation. We also note differential DNA methylation near 5 genes involved in inflammation, immune response, and aging that may link maternal care with stress phenotype.

obligately social (unlike the other species mentioned) and therefore may be the exception that proves the rule.
-Candidate gene study: What is the window size and other features of the putative promoter for the glucocorticoid receptor? Does it overlap a CpG island or contain more CpG sites in other species? The relative lack of CpG sites and DNA methylation may be explained by the lack of a CpG island in the promoter of this species, the promoter was misidentified, or the bisulfite sequencing was off target. Given the surprising nature of this result, I think more detail may help a reader make sense of what happened or, if unable to disentangle possible causes, maybe it should be omitted to avoid unnecessary confusion.
- Figure 3: The DAG figures are useful to understand the study design, but would it be possible to integrate the results into the same figures such as by weighting the arrows based on the effect sizes and including an annotation for which were significant? Otherwise, upon a quick skim, it may appear as if these are the observed connections rather than all of the links tested.
Reviewer #2 (Remarks to the Author): This paper was a pleasure to read and review. It describes a highly interesting set of studies with relevance to a broad range of disciplines. Of particular interest is the examination of the impacts of maternal behavior and social environment on epigenetic marks and corticosteroid levels in a complex social and ecological environment. As so much of both the epigenetic and stress literature relies on either rodents in highly deprived artificial environments or human studies which have to cope with high degrees of unobserved variability and reliance on self-report, the study should help resolve some long standing questions with regard to development and stress. I'm particularly struck by the fact that maternal behavior has such a limited effect on adult corticosteroid levels, while the richness of social contacts outside of the den appears to be more powerful. I think this may begin to explain why so many early life interventions in humans fail to produce substantial effects across the lifespan, and point the way to new developmental windows in which to intervene. The paper was well written and should be published with a commentary calling attention to it. I have only a few substantive concerns: 1. P-7. Description of the LUMA assay-states: "The majority of CpG sites assessed…occur in gene bodies … as well as in non-coding regions of the genome…" As this description basically covers the entire genome, I'm not sure how meaningful the argument is. I think it might be better to simply cite some references that associate global DNA methylation with either health, epigenetic age or some equivalent variable.
3. For the EWAS analysis (and more generally), why was the domestic cat used when a draft genome of the both the striped and spotted hyenas exist (Yang et al 2020)? Even if these draft assemblies are too spotty to be useful in all cases, it should be noted they do exist and explained why they were not used.
Similarly, it would increase my confidence in the data if the approach of using the cat genome was validated using another well-annotated carnivore such as the domestic dog.
4. With regard to the glucocorticoid receptor data, is it possible the site examined was a pseudogene? If not, it might be useful to give a bit more context in terms of hyena endocrinology, which is relatively unique, and might explain why this species appears to have no plasticity at this gene, which seems to have plastic expression/epigenetic regulation in most other mammalian species thus far examined.
Yang C, Li F, Xiong Z, Koepfli KP, Ryder O, et al. 2020. A draft genome assembly of spotted hyena, Crocuta crocuta. Sci Data 7: 126 Reviewer #3 (Remarks to the Author): Laubach et al.'s interesting manuscript entitled Associations of early social experience with offspring DNA methylation and later life stress phenotype represents an important investigation of the role of social interactions/connectedness on subsequent genetic (DNA methylation) and stress responses (fecalbased glucocorticoid metabolites). The use of a novel model species ---the spotted hyena---in a naturalistic setting brings much needed authenticity to previous investigations in this research area that are based on laboratory bred animals raised in restricted laboratory environments. Further, the importance of social interactions in natural hyena behavior and the relevance of their social hierarchy on survival-related resources make this species ideal for the proposed research in the current study. It's not surprising that the data depart slightly from earlier findings in related studies focusing on maternal behavior and subsequent hypothalamic-pituitary-adrenal (HPA) activity in developing offspring. For example, it was reported in Laubach's manuscript that observations of maternal care early in hyena offspring development seemed to be less impactful than reported in rats [e.g., by Michael Meaney's team's work on low-licking and high-licking maternal rats (e.g., see Caldji et al.,1998, Maternal care during infancy regulates the development of neural systems mediating the expression of fearfulness in the rat. Proceedings of the National Academy of Sciences of the United States of America, 95(9), 5335-5340)]. However, by focusing further in the developmental phase (i.e., Den Independence phase), social interactions were indeed found to be associated with fecal metabolites in a direction indicating that increased social interactions mitigated the glucocorticoid response. Still, Laubach et al.'s findings suggest that the maternal-offspring interactions may be more complex than previously observed in the rodent studies. For example, it is possible that more vulnerable offspring may be viewed to elicit more maternal responses than healthier animals when intra-and inter-litter differences are considered. The limitations of collecting fecal samples at various times throughout the day were addressed by treating the time of collection as a covariate in the statistical analyses.
Looking forward, it would be interesting to also include DHEA assays when processing fecal samples, as considering DHEA/CORT ratios provide another index of adaptive stress responsivity since DHEA has been associated with emotional resilience, among other functions (e.g., see: Sripada RK, Marx CE, King AP, et al. DHEA enhances emotion regulation neurocircuits and modulates memory for emotional stimuli. Neuropsychopharmacology. 2013;38(9):1798-1807. The fact that the samples weren't collected after a stress challenge limits the information about responsivity but this is addressed in the manuscript. Looking back in the data, there may have been isolated times after a natural stressor that may be interesting to try to tease apart from the pooled baseline data.
Focusing on the findings that higher social connectedness during the maternal and later developmental phase affects global DNA methylation corroborate previous work, and enhances the robustness of these earlier findings with ecological relevant data in a unique species model. Concerning the scoring of maternal behavior, it wasn't clear if proximity of the infant to the mother was viewed independently from the nursing and grooming behavioral categories. Was proximity scored when the animals were in proximity but NOT engaging in those behaviors? Was social play considered in the early behavioral observations? The sophistication of the statistical approaches in this manuscript (e.g., linear mixed model) is impressive and appears to consider the most likely causal paths with the social connectedness-stress response scenario. With the data being drawn from a pool of observations/samples extending from 1988-2016, is it necessary to describe other studies published from these data? I see that one study apparently published from this data set is mentioned in the current manuscript (i.e., Laubach et al., 2019, Mol Ecol), but doesn't explicitly convey that the findings were in the same group of animals/data.
To succinctly address the prescribed questions highlighted for reviewers, the findings described in the current manuscript are indeed noteworthy as they provide a naturalistic, ethologically-relevant approach to genetic and endocrine markers of emotional resilience and stress. This approach qualifies specific aspects of past findings with rodent models. The findings support the claims in the study and, although the methodology has limitations in areas (e.g., not being able to control time of fecal sample collections to control for circadian effects), the authors cleverly compensated for these variations with statistical techniques. The methods are clearly described in the manuscript and supplemental materials.
Reviewer #4 (Remarks to the Author): In the manuscript entitled "Associations of early social experience with offspring DNA methylation and later life stress phenotype" the authors chose wild spotted hyenas (Crocuta crocuta) as model species aiming to decipher the impact of early life experiences on faecal glucocorticoid metabolite levels and DNA methylation. The authors also test the hypothesis that DNA methylation may act as physiological mediator among early life experience and stress later in life. The study was inspired by former studies in laboratory model species under controlled experimental set-ups. Early experiences were measured by observation of maternal care which were specified as grooming, nursing and close proximity to the offspring`s mother during two chosen time windows in early development "communal den (CD) phase" and the subsequent "den-independent (DI) phase". Unfortunately, it remains unclear why these specific windows were chosen, as well as their lengths in months. The biological meaning of splitting "maternal care" into the three variables: grooming/licking, nursing and close proximity needs further explanation. Figure 6 shows that only close proximity of offspring to the mother has an effect on global methylation, but neither grooming nor nursing. The biological meaning of this but not the others remain unexplained in the discussion, where the authors refer to this result as maternal care, instead of to close proximity. "We also found that more maternal care in the first year of life, which roughly corresponds with the CD period of development, ….". The same strategy of gathering the results which were before specified holds true for measures of the CD period, stated in the second part of this sentence "…and greater social connectedness during the DI phase were each associated with higher global DNA methylation, a presumed indicator of genomic stability and overall health.", while figure 7 shows only a slight effect for the "degree" of interactions on the global methylation ratio, but not the others. Why were these variables separated into more specific variables, but explained as one? What is the biological relevance of these specifications? Furthermore, the specific maternal care observatory data were measured in 1 to 29 repeated observations of mother-offspring interactions. These variable estimated values were subsequently reduced to a single value (by calculating the distance to the mean of each measurement), erasing its variation, and thus reducing the complexity of the data. This approach has been shown to be problematic as this strategy leads to a decrease in robustness of the model these values are used in (https://academic.oup.com/beheco/article/28/4/948/3059669). Also it is unclear why certain covariates were used in respective models e.g. one includes social rank, but not the others. Please explain how they were chosen. The offsprings' sex seems critical, because females are philopatric while males disperse. As a result it may lead to differences in early behaviour. Additionally, puberty in males occurs earlier (after 18-24 months) than in female spotted hyenas (after 24-30 months), which may also lead to differences in early behaviour, and has been shown to change the overall DNA methylation in human. Please provide the number of males and females in the manuscript.
A comparison of a wild mammal to non-model species seems like the attempt to simplify a highly complex system within an uncontrolled experimental set-up. This simplification is critical. The authors are partly aware of this as they state within the discussion "…, in contrast to controlled rodent experiments, wild hyenas are subject to a multitude of stressors over development, which may hamper our ability to isolate the specific effect of maternal care on stress phenotype. Second, compared to maternal separation common in many experimental studies, natural variation in maternal care is far more subtle and may require a more sensitive measure of physiological stress than average fGCMs, which is a summary baseline stress indictor 51 subject to unmeasured environmental factors and the animal's condition when the sample is collected 52." The measurement of faecal glucocorticoid metabolite is a valid method. But stress levels are prone to outer environmental changes e.g. predatory attacks and human interactions increase stress levels. Accounting for these events is often difficult due to limited time windows of observation. "Behavioral data were collected daily between 0530 -0900 h and 1700 -2000 h." Therefore replicates of fGCMs for each animal are important as the period and day of sampling is critical. Please state if they were applied.
Global methylation analysis using LUMA only provides a methylation ratio over the entire genome. It does not provide the data for functional interpretation. The authors state "The majority of CpG sites assessed via LUMA occur in gene bodies, where they may function in transcription regulation and alternative splicing 27, as well as in non-coding regions of the genome, where they may repress repetitive elements 28 and enhance chromosome stability 29. Therefore, we assumed that lower than average %CCGG methylation is likely disadvantageous to health." This assumption is very farfetched, because LUMA does not provide the details to distinguish between functional important regions and mostly excludes promoters, which are thought to be the main regulatory regions. Therefore methylation changes are difficult to interpret and should be handled with caution. Using LUMA is cost-efficient and does not need a reference genome, which would be two valid reasons to use it. However the usage of three different approaches with different samples numbers: global DNA methylation (n = 186), genome-wide DNA methylation (n = 29), and candidate gene DNA methylation (n = 96) without explanation, is confusing. Please clarify why you chose this strategy.
In order to assess the biological meaning of the established data, please clarify the numbers of animals, mothers, sex of offspring, and siblings, and the pairwise interactions observed. The number of samples collected for CD and DI time windows, stating why these were chosen, the sex, the cofounding factors that can lead to an increase in stress, the robustness of the metabolites measured. Last but not least, providing the data and code used is essential to understand and repeat the analysis.
We want to begin by thanking all of the reviewers for their thoughtful and timely review of our manuscript. We have carefully considered the reviewers' comments below and sought counsel from multiple statisticians (n=2), bioinformaticians (n=1), genetic/genomic epidemiologists (n=2) and evolutionary biologists (n=1) in order to address the issues identified by the reviewers, particularly in relation to the global and genome-wide DNA methylation analyses. We addressed specific comments below by carefully weighing the consensus opinions of the reviewers, the experts with whom we consulted, and drawing on the specializations represented in our author list. In pooling expertise from multiple sources, we acknowledge that multiple approaches are employed in studies like ours and that this is an area of rapid methodological development (e.g. causal inference, statistical analyses of genomics etc.), and so the methodologies employed herein represent what we collectively determined to be the most appropriate options. We feel that the changes made based on the reviewers' comments have improved this work, and we respond to specific concerns below.
Reviewer #1: 1. Comment: This manuscript is seeks to test the hypothesis that DNA methylation represents a mechanistic link between early social experience and a biomarker of stress later in life. The section evaluating whether unfavorable early social experiences affect later less stress physiology is particularly clear and an important contribution to the field. However, the later sections examining global DNA methylation and genome-wide DNA methylation appear inconclusive and are difficult to interpret. Specifically (see following comments):

Response:
We are grateful to the reviewer for their careful review of this manuscript, and we are delighted that they find value in this work as a contribution to the field.  (Woo & Kim, 2012). It is worth noting that higher global DNA methylation has also been associated with some cardiovascular diseases (Sharma et al., 2008). Thus, deviations in either direction, both lower or higher DNA methylation, may be a marker of aberrant physiology and should be interpreted with caution.

Pg 35, Lines 589-977:
We found no evidence of mediation by %CCGG methylation. This could be due to low power given the relatively small sample sizes, and/or related to the fact that %CCGG methylation averaged over the entire genome may be too broad a metric to isolate important regulatory pathways involved in the stress response. For example, we previously found that maternal prenatal socioeconomic status was associated differential methylation at specific regions of the genome but not with global DNA methylation in humans (Laubach, Perng, et al., 2019). Therefore, we complemented our global DNA methylation analysis in hyenas with mERRBS (Meissner et al., 2005), a high throughput sequencing technology that assays single nucleotide DNA methylation resolution at a genome wide scale.

Comment:
The absence of this information makes it difficult to evaluate whether the reader would expect global DNA methylation to be a good biomarker and how to interpret the small effect sizes observed.

Response:
Previous studies have shown larger effect sizes when assessing global DNA methylation in relation to cancer, and as a consequence of inadequate methyl donor intake. However, these scenarios represent extreme cases of genomic instability and nutrient restriction, respectively, and may not be relevant to the majority of individuals in a population. The goal of this study was to understand whether natural variation in early social experience is associated with naturally-occurring differences in global DNA methylation, especially given earlier findings that the hyenas' maternal social rank and food availability correspond with similar magnitudes of effect sizes when %CCGG methylation is the outcome (Laubach et al. 2019. Mol. Eco. 28: 3799-3812). The effect sizes we report in this analysis (β = 1% to 2%) are modest but represent differences detectable across many CCGG sites (~2.1 million) throughout the genome and thus, likely have implications for health and potentially fitness. Moreover, the differences we observed in this analysis are comparable or larger than those reported in observational studies of healthy humans. We also acknowledge that larger differences in %CCGG methylation have been reported in studies that do not focus on cancer -e.g., ~6% lower DNA methylation in brain tissue of wild polar bears in response to mercury exposure (

Comment:
As the authors note near the end of the discussion, cell type variation is an important source of variation that must be accounted for in modeling DNA methylation levels (such as by estimating proportions of major cell types using blood smears).
Response: Unfortunately, we do not have information on the cell type composition, which we discussed as a potential limitation. However, the absence of these data does not affect interpretation of our results given that our goal was to estimate the total effect of each explanatory variable on %CCGG methylation. Cellular composition is likely on the causal pathway given that social stressors have been identified as determinants of cell type composition (Engler et al. 2004. J. Neuroimm. 148: 106-115) which in turn, is a known determinant of DNA methylation. In such a scenario, adjustment for cell type composition would introduce rather than reduce bias into the estimate of interest (Ananth & Schisterman. 2017. Am. J. Obst. Gyn. 217: 167-175;Schisterman et al. 2009. Epid. 20: 488-495). To address this comment and provide further clarity, we have modified the following text in the Discussion section:

Discussion:
Pg 39, Lines 682-687: Fourth, we used archived DNA from whole blood that lacked information on cell type composition, and therefore we were not able to control for cell type heterogeneity as a source of variation in DNA methylation. However, cell type composition may be influenced by factors like social stress (Engler et al., 2004), and therefore could be on the causal pathway between the explanatory variables of interest and DNA methylation. Accordingly, adjusting for such a variable in the analysis would introduce rather than reduce bias (Ananth & Schisterman, 2017;Schisterman et al., 2009).

Comment:
Similarly, genetic covariance between samples is important and can dramatically effect results, especially when only considering a few individuals, and should be accounted for even if only from pedigree records or genotypes called from bisulfite converted DNA (i.e. bisSNP).

Response:
In response to this comment, we ran microsatellite assays to determine paternity, used these data to construct a pedigree and a covariance matrix of relatedness, and included the matrix in our EWAS models to control for genetic correlations. These updates are described in text (Methods: Pages 17-18, Lines 299-301; SI: Pages 16-17, Lines 301-307). For analyses in which global (%CCGG) DNA methylation was the outcome, we included a random intercept for maternal ID to account for familial clustering due to shared genetics. We have noted these changes in the main and supporting texts.
6. Comment: 29 samples is small for a genome-wide association study, especially given the biological noise in estimating DNA methylation and fGCs. I believe a power-analysis in the SI would be extremely useful.
Response: While our sample size is small (which we noted as a limitation, Conclusions: Page 39, Lines 680-682), it is comparable to other recent genome-wide DNA methylation studies in wild animals (c.f. Taff  Further, the small sample was likely not prohibitive given that we were able to detect meaningful differences after false discovery rate correction. Of note, we used strict sample selection criteria to reduce noise in the data (Methods: Pages 16-17, Lines 278-282), thereby reducing the need for covariate adjustment and enhancing statistical power.
Regarding the suggestion to include a power analysis in the Supplemental Material: the EWAS portion of this analysis was exploratory, so we did not conduct a priori power analysis. Additionally, we do not feel it is appropriate to conduct post hoc power analysis given that we were able to detect significant effects after false discovery rate corrections, and because post hoc power analyses are prone to type 2 error (c.f. Bababekov et al. 2018. Ann. Surg. 267: 621-622, and a less-technical piece written by Gelman: http://www.stat.columbia.edu/~gelman/research/unpublished/power_surgery.pdf).
7. Comment: On the topic of power, it is fairly typical to exclude CpG sites with very low or high methylation (e.g. <10% or >90%), which may increase signal by reducing the multiple-testing burden.
Response: Thank you for the suggestion. We are aware of research groups who implement this initial data processing step (Lea et al. 2016. Mol. Eco. 25: 1681-1696. Following this precedent, we have removed 5% of CpG sites which exhibited the lowest variation in DNA methylation, and those that had <10% DNA methylation or >90% DNA methylation in order to reduce the burden of false positives. We then conducted the EWAS using this filtered data set and have made changes describing the data filtering and new analysis in the main and supporting text (Methods: Page 17, Lines 286-288; SI: Page 16, Lines 292-300; SI, Page 23, Lines 394-396).
8. Comment: Furthermore, the distribution shown in Figure S8b is depleted for low p-values and not nearly uniform, suggesting a permutation-based null and FDR correction may be a better fit than the typical Benjamini-Hochberg correction. I think this would be apparent in Figure S8a if the y=x line was added as well, although there may be an enrichment over a permutation-based null.
Response: Thank you for this insightful point. After additional consultation with experts in genomics and statistics, we agree that it is appropriate to control for inflation and bias as indicated by our original P-value distribution and QQ-plot. We used the Bioconductor package, 'bacon' to correct estimates and associated P-values for inflation and bias. We describe these changes in the main and supporting text, and we have provided additional supplemental figures that highlight improvements after correcting for bias and inflation.

Methods: Pg 18, Lines 301-305:
After running the EWAS, we corrected estimates of association and P-values for inflation and bias due to unmeasured confounding using the Bioconductor package, 'bacon' in order to reduce false positive associations (van Iterson et al., 2017). We accounted for multiple comparisons using a Benjamini-Hochberg false discovery rate (FDR) of 5% on the inflation and bias corrected P-values (Hochberg & Benjamini, 1990).

Supporting material:
Pg 24 (SI), Lines 403-408: We graphed the P-value distribution (Supporting Figure 9, a) and a Q-Q plot of the expected vs. the observed P-values (Supporting Figure 9, b) from our uncorrected EWAS. These plots indicated the presence of test statistic bias and inflation, likely due to unmeasured confounding, so we used the package 'bacon' to improve estimation of the empirical null distribution (van Iterson et al., 2017). After applying the 'bacon' correction we no longer observed depletion of low P-values (Supporting Figure 9). Similarly, we annotated each differentially methylated site from our EWAS to each of these respective regions (c.f., Supporting Figure 11). We describe these methods and results and discuss their biological meanings in respective sections of the main text and supporting material. However, the present analysis is not an ideal setting in which to formally compare global vs genome-wide results for a few reasons. First, the overlap in the global and genome-wide DNA methylation samples includes only 12 animals, therefore a direct comparison between data sets is limited. Second, the purpose of the EWAS analysis was exploratory, and to identify potential functional pathways that can be assessed more rigorously in future analyses with larger sample sizes. At this point, we hesitate to add additional analyses on top of multiple existing analyses given the limits of article lengths for the journal and because of the already-dense nature of the paper. However, we mention these ideas in our revised manuscript as a suggestion for future research.

c) d)
Minor points: 10. Comment: -Abstract: it's unclear what "DNA methylation biomarkers" means, especially considering that, at this point, there has been no mention of genome-wide DNA methylation. I would suggest rephrasing as "Finally, using bisulfite sequencing data, we identify XXX CpG sites".

Response: Done
Abstract: Pg 2, Lines 8-10: Using bisulfite sequencing data, we also identified differential DNA methylation near 5 genes involved in inflammation, immune response, and aging that may link maternal care and stress phenotype.
11. Comment: Introduction: It's mentioned that "few studies have measured natural variation…" but gives no mention of those that do. A few citations regarding where this has been done, and therefore what insights can be gained, may be helpful here. Similarly, "few have measured all three elements in the same population" would benefit from examples which do so which could be integrated into the explanation of why it is variable. Currently the only citation is to Weaver et al. 2014 which, while foundational, does not utilize a genome-wide approach.

Response:
We have now cited relevant studies in the Introduction, as suggested:

Introduction: Pg 4, Lines 45-57:
First, there is a lack of studies that measure natural variation in the quantity and quality of early social experience, as much of this research includes experimental studies involving maternal separation and peer isolation, which are routine procedures in studies of captive primates and rodents 14,15 . While informative, studies of extreme deprivation do not capture the range of normative social experiences that are relevant to development 3,16,17 . Second, although studies have examined the relationships among social experiences, DNA methylation, and stress phenotype in piecemeal, few have measured all three components in the same population 10,18-20 . Doing so is necessary to explicitly test the hypothesis that DNA methylation represents a mechanistic link between early social experience and future stress phenotype. Third, few studies explicitly consider that variation in the type and/or timing of social experiences may affect DNA methylation 21 and stress phenotypes 22 , particularly in wild animal populations where development is subject to environmental fluctuation and pressures of natural selection.
12. Comment: Methods: It's mentioned that "the majority of CpG sites assessed via LUMA occur in gene bodies", but is that actually true or are gene bodies just strongly enriched? It would be helpful to extract the location of CCGG sequences and report what percent are contained in promoters, gene exons, gene introns, CpG islands, CpG shores, etc. to know which regions of the genome LUMA is weighted by.
Response: Thank you for noting this point requiring clarification. As suggested, we have annotated the 'CCGG' motif in the hyena genome and addresses this point in response to the reviewer's previous comment #9.

Comment:
Similarly, what was the coverage of the genome-wide DNA methylation data and the distributions for site-specific mean methylation and the genomic context of measured CpG sites?
Response: Thank you, we have added new figures and text to provide added genomic context, (c.f. Supporting Figure 11, and additions to the methods, and results.)

Methods:
Pg 8, Line 139: We filtered read counts to a minimum of 10x coverage.

Results:
Pg 26, Lines 448-455: Based on the CpG density annotation, we observed 62% of CpG sites from our EWAS occurred in CpG islands and CpG shores, while among both hyper-and hypomethylated DMSs 72% of sites were in CpG islands and shores (Figure S11, a, c). According to our genic annotation, 38% of CpG sites from our EWAS occurred in either promoters, exons or introns, and 39% of DMSs were found in the same genic annotations (Figure S11, b, d). There were equal total numbers of hyper-and hypomethylated DMSs, but in both annotations more hypermethylated DMSs were found near CpG islands and gene bodies than in intergenic regions.
14. Comment: What is the NCBI link/published citation for the hyena genome assembly used? If it hasn't been done already, it would be useful to have genes annotated via NCBI or ensemble. UCSC also has publicly available tools to identify CpG islands which are very informative for analyses of DNA methylation and not presented here.

Response:
The DNA methylation calling for the EWAS used the draft spotted hyena genome (at that time unpublished). We have now added the appropriate citation. We also address the identification of CpG islands in response to comment #9.

Methods: Pg 8, Lines 137-139:
After DNA sequencing, we ran a standard bioinformatic pipeline to clean, align, and call DNA methylation reads using the draft spotted hyena genome (Yang et al., 2020).
15. Comment: What ere the parameters used for blat, and how were multiple complete/best hits resolved? Was there a maximum distance for the nearest gene?

Response:
We have added text clarifying the criteria we used for BLAT. We have also noted the approximate distance of each DMS based on the multiple species alignment from our BLAT query in Table 2.

Methods:
Pg 20, Lines 345-350: We used default BLAT parameters (see Supporting Materials) and retained the highest scoring DNA sequence matches in cases where there was concordance between the multiple species alignment when querying the cat and human reference genomes. We do not report the nearest human gene when there was discordance between the cat and human queries and when the BLAT score and % identical sequence matches were low (<500 BLAT score, <85% match). 17. Comment: Following the EWAS analysis between DNAm and fGCMs, it appears as if only 15 sites were tested for associations with maternal care and rank. This type of partial analysis based on arbitrary "significant" sites reduces the ability to identify shared effects (Urbut et al. 2018, Nat Gen) and also limits the ability to test whether such effects are more likely than expect by chance. I would suggest doing the full EWAS for the other predictor variables and testing for overlap/enrichment and whether the subset of sites identified truly survive multiple hypothesis testing.
Response: While carrying out full EWAS on predictor variables (of which there are several) and fGCMs is certainly an option, this approach may be particularly prone to type 1 error for this study. This is because we have multiple early life social experience variables, so we would need to conduct multiple EWAS in addition to the stress related EWAS in order to be able to look for overlapping DMSs. Given the already small sample size for this portion of the analysis, we took the more efficient meet-in-themiddle approach of conducting a single EWAS with FDR correction to identify determinants of the outcome, and then assessing for associations with a key early life variable based on associations detected in the present analysis. We have used this approach in other analyses of high-dimensional data, including an epigenome-wide association study of prenatal SES and DNA methylation at three 18. Comment: Discussion: Overall this was really well written and did a good job discussing the results and the weaknesses. I would just point out that the yellow-bellied marmot example is a species that is not obligately social (unlike the other species mentioned) and therefore may be the exception that proves the rule.

Response:
We have changed the text to make this clarifying point.

Discussion:
Pg 32, Lines 525-528: For instance, a study on wild yellow-bellied marmots found no association of social network metrics with fGCMs measured in a mixed-age, mixed-sex population (Wey & Blumstein, 2012), although it should be noted that this species is not obligately social like those in other above-mentioned studies.
19. Comment: Candidate gene study: What is the window size and other features of the putative promoter for the glucocorticoid receptor? Does it overlap a CpG island or contain more CpG sites in other species? The relative lack of CpG sites and DNA methylation may be explained by the lack of a CpG island in the promoter of this species, the promoter was misidentified, or the bisulfite sequencing was off target. Given the surprising nature of this result, I think more detail may help a reader make sense of what happened or, if unable to disentangle possible causes, maybe it should be omitted to avoid unnecessary confusion.

Response:
We have added information about the size and the %GC content of the putative spotted hyena GR promoter in the figure legend of Supporting Figure 5. In addition, this figure delineates important features including the location of potential NGFI-A transcription factor binding sites as well as consensus sequence based on alignments human and rodent DNA sequences. While we think these null results from hyenas are important to disseminate given that no prior study done this and because it contributes valuable information on the hyena GR promoter region, we are open to removing this finding if the reviewers and editors feel it is appropriate.

Supporting material:
Pg 11 (SI), Lines 190-199: Supporting Figure 5. Spotted hyena putative GR promoter sequence, which has been trimmed to match the aligned sequences published in Oberlander et al. 2008, McGowan et al. 2009, and Perroud et al. 2011. This sequence is 580 bp and has a 73.1% GC content. Underlined is the inferred consensus region (based on alignment) for human exon 1F, and red, italicized is the inferred consensus region (based on alignment) for rat exon 17. GC boxes, the putative binding site for transcription factor NGFI-A, are bold and contained CpG sites capitalized. The yellow highlighted GC box is the potential NGFI-A transcription binding site homologous to the region in rats that was shown to be differentially methylated with respect to maternal licking and grooming (c.f. Weaver et al. 2004).
20. Comment: Figure 3 -The DAG figures are useful to understand the study design, but would it be possible to integrate the results into the same figures such as by weighting the arrows based on the effect sizes and including an annotation for which were significant? Otherwise, upon a quick skim, it may appear as if these are the observed connections rather than all of the links tested.
Response: As discussed in our recent publication (Laubach et al. 2021. Proc. Royal Soc. B. 288: 2020815), the purpose of DAGs is to depict the temporal and causal relationships among the variables of interest. These graphs are conceptual and based on a priori knowledge. To clarify, we have updated all figure legends to reflect this (c.f. the updated legend for Figure 1).

Figure 1. Directed Acyclic Graph (DAG) for Part 1 of the analysis showing the hypothesized relationships between early-life social experience and adult fecal Glucocorticoid Metabolites (fGCMs).
Covariates were included based on a priori biological knowledge and bivariate analysis (see Supporting  Tables S2-S5). Precision covariates (gray text) were assessed when the fecal sample was collected. Variable positions over the timeline roughly corresponded with the timing of assessment.

Dependent variable
Reviewer #2 1. Comment: This paper was a pleasure to read and review. It describes a highly interesting set of studies with relevance to a broad range of disciplines. Of particular interest is the examination of the impacts of maternal behavior and social environment on epigenetic marks and corticosteroid levels in a complex social and ecological environment. As so much of both the epigenetic and stress literature relies on either rodents in highly deprived artificial environments or human studies which have to cope with high degrees of unobserved variability and reliance on self-report, the study should help resolve some long standing questions with regard to development and stress. I'm particularly struck by the fact that maternal behavior has such a limited effect on adult corticosteroid levels, while the richness of social contacts outside of the den appears to be more powerful. I think this may begin to explain why so many early life interventions in humans fail to produce substantial effects across the lifespan, and point the way to new developmental windows in which to intervene. The paper was well written and should be published with a commentary calling attention to it. I have only a few substantive concerns.
Response: Thank you for thoughtfully reviewing this manuscript. We are grateful that the reviewer is excited about this work, and in particular, appreciates the value of the wild study population.
2. Comment: P-7. Description of the LUMA assay-states: "The majority of CpG sites assessed…occur in gene bodies … as well as in non-coding regions of the genome…" As this description basically covers the entire genome, I'm not sure how meaningful the argument is. I think it might be better to simply cite some references that associate global DNA methylation with either health, epigenetic age or some equivalent variable.
Response: This is a salient point and germane to our discussion of the results. Reviewer 1 has also noted this issue, so we refer to our previous responses to Reviewer 1's comments #2, #7, and #10. More specifically, in Rev. 1, comment #2 we discuss associations between global DNA methylation and various health outcomes and provide additional citations in text. In Rev. 1 comments #7 and #10, we present new analyses that annotate the 'CCGG' motif in the recently published draft spotted hyena genome. These new analyses and associated figures provide a much more detailed description of what the LUMA assay measures in this species.
Regarding the lack variation/plasticity in DNA methylation in hyenas: we suspect that this may partially reflect the fact that we used leukocytes instead of brain tissue, the latter of which seems to be more variable in GR DNA methylation. To address this, we exercise caution in interpreting our null results and suggest future studies consider a more complete sequencing and testing of differential CpG methylation across the entire GR gene.

Discussion: Pg 38, Lines 652-655:
Future studies may better replicate results from laboratory rodents by using brain tissue, wherein expression of the GR gene is more tightly coupled with HPA regulation or by exploring DNA methylation of additional CpGs with more extensive coverage of the hyena GR gene.

Reviewer #3
1. Comment: Laubach et al.'s interesting manuscript entitled Associations of early social experience with offspring DNA methylation and later life stress phenotype represents an important investigation of the role of social interactions/connectedness on subsequent genetic (DNA methylation) and stress responses (fecal-based glucocorticoid metabolites).The use of a novel model species ---the spotted hyena---in a naturalistic setting brings much needed authenticity to previous investigations in this research area that are based on laboratory bred animals raised in restricted laboratory environments. Further, the importance of social interactions in natural hyena behavior and the relevance of their social hierarchy on survival-related resources make this species ideal for the proposed research in the current study. It's not surprising that the data depart slightly from earlier findings in related studies focusing on maternal behavior and subsequent hypothalamic-pituitary-adrenal (HPA) activity in developing offspring. For example, it was reported in Laubach's manuscript that observations of maternal care early in hyena offspring development seemed to be less impactful than reported in rats [e.g., by Michael Meaney's team's work on low-licking and high-licking maternal rats (e.g., see Caldji et al.,1998, Maternal care during infancy regulates the development of neural systems mediating the expression of fearfulness in the rat. Proceedings of the National Academy of Sciences of the United States of America, 95(9), 5335-5340)]. However, by focusing further in the developmental phase (i.e., Den Independence phase), social interactions were indeed found to be associated with fecal metabolites in a direction indicating that increased social interactions mitigated the glucocorticoid response. Still, Laubach et al.'s findings suggest that the maternal-offspring interactions may be more complex than previously observed in the rodent studies. For example, it is possible that more vulnerable offspring may be viewed to elicit more maternal responses than healthier animals when intra-and inter-litter differences are considered. The limitations of collecting fecal samples at various times throughout the day were addressed by treating the time of collection as a covariate in the statistical analyses.

Response:
We thank the reviewer for their time and effort spent vetting this work. We are pleased that they find this analysis to yield interesting findings.
2. Comment: Looking forward, it would be interesting to also include DHEA assays when processing fecal samples, as considering DHEA/CORT ratios provide another index of adaptive stress responsivity since DHEA has been associated with emotional resilience, among other functions (e.g., see: Sripada RK, Marx CE, King AP, et al. DHEA enhances emotion regulation neurocircuits and modulates memory for emotional stimuli. Neuropsychopharmacology. 2013;38(9):1798-1807). The fact that the samples weren't collected after a stress challenge limits the information about responsivity but this is addressed in the manuscript. Looking back in the data, there may have been isolated times after a natural stressor that may be interesting to try to tease apart from the pooled baseline data.
Response: Thank you for this suggestion. We agree that there are a number of additional metrics of stress responsivity that would provide valuable insight and potentially help disentangle the adaptive function of developmentally plastic stress physiology. We noted that we were not able to measure change in stress hormone levels in response to a standardized challenge as a limitation, and we suggest that future studies should consider quantifying this difference in response to the darting event as has been done by Dr. Robert Sapolsky in baboons. Quantifying DHEA and modeling the ratio of CORT/DHEA is another good idea. We have added this suggestion in the discussion of future directions.

Discussion:
Pg 39, Lines 677-680: However, a more nuanced measure of stress, like reaction and recovery to a standardized stressor or the ratio of cortisol to dehydroepiandrosterone (CORT/DHEA), may provide a more direct assessment of acute HPA function (Kamin & Kertes, 2017;Romero, Dickens, & Cyr, 2009) over time.
3. Comment: Focusing on the findings that higher social connectedness during the maternal and later developmental phase affects global DNA methylation corroborate previous work, and enhances the robustness of these earlier findings with ecological relevant data in a unique species model.

Response:
We agree and thank the reviewer for this comment. We have slightly modified the text to clarify this point.

Discussion:
Pg 33, Lines 547-551: Greater maternal care (specifically, close proximity received during the first year of life) and higher social connectivity (namely, degree during the DI life stage) were both associated with higher global DNA methylation later in life. These results corroborate those of previous studies that report associations of early social experiences with differential patterns of DNA methylation in rodents, Rhesus macaques, and humans (Anier et al., 2014;Provencal et al., 2012;Unternaehrer et al., 2015).

Comment:
Concerning the scoring of maternal behavior, it wasn't clear if proximity of the infant to the mother was viewed independently from the nursing and grooming behavioral categories. Was proximity scored when the animals were in proximity but NOT engaging in those behaviors? Was social play considered in the early behavioral observations?
Response: Thank you for the opportunity to clarify. Close proximity includes time spent nursing or grooming. We have updated the text in the supplemental information to reflect this detail and discussed this as a potential reason for the strong associations of close proximity with global DNA methylation. Social play was not considered in these behaviors.

Discussion:
Pg 34, Lines 560-565: One potential explanation for the robust effect of close proximity but not the other metrics of maternal care (nursing and grooming) on global DNA methylation may reflect the fact that close proximity includes the time mother hyenas spent nursing and grooming. Thus, close proximity is a composite metric of maternal investment in offspring that may capture more variation in offspring global DNA methylation.

Supporting material: Pg 1 (SI), Lines 9-15:
We focused on three specific maternal care behaviors: time spent in close proximity, nursing, and grooming. We chose these specific metrics of maternal care because comparable behaviors in primates and rodents have been shown to influence offspring behavior and physiology, albeit in captive settings. Additionally, these behaviors capture the effect of nutrition as well as physical contact on offspring development. Nota bene, maternal care behaviors were not mutually exclusive in their occurrence; for example, a portion of time spent in close proximity includes time spent nursing or grooming.

Comment:
The sophistication of the statistical approaches in this manuscript (e.g., linear mixed model) is impressive and appears to consider the most likely causal paths with the social connectedness-stress response scenario. With the data being drawn from a pool of observations/samples extending from 1988-2016, is it necessary to describe other studies published from these data? I see that one study apparently published from this data set is mentioned in the current manuscript (i.e., Laubach et al., 2019. Mol. Ecol.), but doesn't explicitly convey that the findings were in the same group of animals/data.
Response: Thank you. The reviewer is correct that this analysis is based on a subsample of animals who are part of a long-term study of wild hyenas. Some of the individuals in the present analysis were also in previous analyses and are published in manuscripts, including Laubach et al. 2019. Mol. Ecol. 11912: 1413-1427.We cite these previous publications in text.

Methods:
Pg 7, Lines 104-107: For each hyena, we constructed association networks during the CD and DI periods (two life stages during which social interactions were previously identified as key determinants of fitness 29 ) as previously described for this population 26 .

Pg 7, Lines 119-121:
Descriptions of the LUMA assay, laboratory procedures, and data cleaning protocol are available in (Laubach, Faulk, et al., 2019), from which the current global DNA methylation samples were drawn.
6. Comments: To succinctly address the prescribed questions highlighted for reviewers, the findings described in the current manuscript are indeed noteworthy as they provide a naturalistic, ethologically-relevant approach to genetic and endocrine markers of emotional resilience and stress. This approach qualifies specific aspects of past findings with rodent models. The findings support the claims in the study and, although the methodology has limitations in areas (e.g., not being able to control time of fecal sample collections to control for circadian effects), the authors cleverly compensated for these variations with statistical techniques. The methods are clearly described in the manuscript and supplemental materials.
Response: Thank you, again.
Reviewer #4: 1. Comment: In the manuscript entitled "Associations of early social experience with offspring DNA methylation and later life stress phenotype" the authors chose wild spotted hyenas (Crocuta crocuta) as model species aiming to decipher the impact of early life experiences on faecal glucocorticoid metabolite levels and DNA methylation. The authors also test the hypothesis that DNA methylation may act as physiological mediator among early life experience and stress later in life. The study was inspired by former studies in laboratory model species under controlled experimental set-ups. Early experiences were measured by observation of maternal care which were specified as grooming, nursing and close proximity to the offspring`s mother during two chosen time windows in early development "communal den (CD) phase" and the subsequent "den-independent (DI) phase".
Response: Thank you for the detailed and timely review of our manuscript. These comments have been insightful, and no doubt have improved the manuscript.
2. Comment: Unfortunately, it remains unclear why these specific windows were chosen, as well as their lengths in months.

Response:
The time-windows of communal den (CD) and den independent (DI) were chosen to reflect developmental stages that correspond with key life history events in spotted hyenas and are based on previous studies in our population that linked early life social experiences to long-term health and fitness. More specifically, we defined the CD period based on data from daily observations of individual hyenas and their activity patterns, which were confined to the communal den. We also recorded when young hyenas began forays away from the communal den out into their clan's territory, which we noted as the start of the DI period. The duration of the DI period was matched to the duration of the CD period in order to balance our sampling design, as previously done in this study population (Turner et al. 2018. Behav.Eco. Soc. 72). Furthermore, social networks metrics assessed within the CD and DI phases using this method are associated with fitness in this population, supporting the biological relevance of these metrics (Turner et al. 2021. J. Anim. Eco. 90: 183-196). We have updated the main text and the text in the Supporting Material to clarify our justification of the time-windows used in our analyses.

Methods:
Pg 7, Lines 104-107: For each hyena, we constructed association networks during the CD and DI periods (two life stages during which social interactions were previously identified as key determinants of fitness 29 ) as previously described for this population 26 .

Supporting material: Pg 1 (SI), Lines 4-7:
Our maternal care data included 1533 FAS totaling approximately 779 hours of observations of 258 mother-infant pairs when offspring were ≤1 year old, which is the approximate age at weaning in this sample (11.59 months) and among spotted hyenas more generally (Holekamp & Smale, 1998).
3. Comment: The biological meaning of splitting "maternal care" into the three variables: grooming/licking, nursing and close proximity needs further explanation.

Response:
We have added additional explanations in the supporting material to address this comment.

Supporting material: Pg 1 (SI), Lines 9-15:
We focused on three specific maternal care behaviors: time spent in close proximity, nursing, and grooming. We chose these specific metrics of maternal care because comparable behaviors in primates and rodents have been shown to influence offspring behavior and physiology, albeit in captive settings. Additionally, these behaviors capture the effect of nutrition as well as physical contact on offspring development. Nota bene, maternal care behaviors were not mutually exclusive in their occurrence; for example, a portion of time spent in close proximity includes time spent nursing or grooming.
4. Comment: Figure 6 shows that only close proximity of offspring to the mother has an effect on global methylation, but neither grooming nor nursing. The biological meaning of this but not the others remain unexplained in the discussion, where the authors refer to this result as maternal care, instead of to close proximity. "We also found that more maternal care in the first year of life, which roughly corresponds with the CD period of development, ….". The same strategy of gathering the results which were before specified holds true for measures of the CD period, stated in the second part of this sentence "…and greater social connectedness during the DI phase were each associated with higher global DNA methylation, a presumed indicator of genomic stability and overall health.", while figure 7 shows only a slight effect for the "degree" of interactions on the global methylation ratio, but not the others. Why were these variables separated into more specific variables, but explained as one? What is the biological relevance of these specifications?
Response: Thank you for bringing this point to our attention. We have revised our text in the discussion so that our explanations reflect findings with respect to specific early life social experiences, and to address the reason why we may have observed an association with close proximity but not the other maternal care variables.

Discussion:
Pg 33, Lines 547-551: Greater maternal care (specifically, close proximity, received during the first year of life) and social connectivity (namely degree during the DI life stage) were both associated with higher global DNA methylation later in life, enhancing the robustness of previous findings that suggest early social experiences influence patterns of DNA methylation (Anier et al., 2014;Provencal et al., 2012;Unternaehrer et al., 2015).

Pg 34, Lines 560-565:
One potential explanation for the robust effect of close proximity but not the other metrics of maternal care (nursing and grooming) on global DNA methylation may reflect the fact that close proximity includes the time mother hyenas spent nursing and grooming. Thus, close proximity is a composite metric of maternal investment in offspring that may capture more variation in offspring global DNA methylation.
5. Comment: Furthermore, the specific maternal care observatory data were measured in 1 to 29 repeated observations of mother-offspring interactions. These variable estimated values were subsequently reduced to a single value (by calculating the distance to the mean of each measurement), erasing its variation, and thus reducing the complexity of the data. This approach has been shown to be problematic as this strategy leads to a decrease in robustness of the model these values are used in (https://academic.oup.com/beheco/article/28/4/948/3059669).
Response: Indeed, the paper by (Houslay & Wilson. 2017. Behav. Eco. 28: 948-952) posits that use of BLUPs as a data reduction technique for the dependent variable of interest may artificially shrink confidence intervals to confer a larger degree of precision than use of use of mixed effects models with random effects to account for correlations between repeated measurements (or clusters of data points). However, a key distinction is that we use BLUPs to summarize repeated measures of the explanatory variables (maternal care behaviors) -not the dependent variable of interest. This data reduction step is necessary for modeling repeated measures data as the independent variable, does not confer any sort of shrinkage effect on the estimate of interest, and has been used in other similar settings (e.g., repeated measures of four LINE-1 DNA methylation sites in an analysis of LINE-1 DNA methylation as a predictor of subsequent growth: Perng et al. 2013. PLoS One. 8: e62587). More specifically, we used this data reduction technique because of the varying number of times (1 to 29 times) that each mother-offspring pair was observed, and because it would be very challenging to assess and interpret associations of behaviors at each individual time-point in relation to the dependent variables of interest.
6. Comment: Also it is unclear why certain covariates were used in respective models e.g. one includes social rank, but not the others. Please explain how they were chosen.

Response:
We selected covariates for multivariable models based on prior knowledge on determinants of stress phenotype (including use of the directed acyclic graphs to visualize the hypothesized temporal and causal associations among variables of interest) and bivariate associations (Supporting Tables 2-5). This approach is widely-used by epidemiologists to efficiently and methodically make causal inference from observational data, as discussed in a recent review that our group led (Laubach et al. 2021. Proc. Royal Soc. B. 288: 2020815). We now mention this in the Methods section (Methods: Page 10, Lines 167-168) 7. Comment: The offsprings' sex seems critical, because females are philopatric while males disperse. As a result it may lead to differences in early behaviour. Additionally, puberty in males occurs earlier (after 18-24 months) than in female spotted hyenas (after 24-30 months), which may also lead to differences in early behaviour, and has been shown to change the overall DNA methylation in human. Please provide the number of males and females in the manuscript.
Response: Done (Page 5, Lines 75-77; and in table and figure legends) 8. Comment: A comparison of a wild mammal to non-model species seems like the attempt to simplify a highly complex system within an uncontrolled experimental set-up. This simplification is critical. The authors are partly aware of this as they state within the discussion "…, in contrast to controlled rodent experiments, wild hyenas are subject to a multitude of stressors over development, which may hamper our ability to isolate the specific effect of maternal care on stress phenotype. Second, compared to maternal separation common in many experimental studies, natural variation in maternal care is far more subtle and may require a more sensitive measure of physiological stress than average fGCMs, which is a summary baseline stress indictor 51 subject to unmeasured environmental factors and the animal's condition when the sample is collected 52." Response: Thanks for the suggestion. We have added the following text:

Discussion:
Pg 33, Lines 540-543: Considering not only the differences in how early-life experience and stress outcomes are measured, but also differences between controlled laboratory settings vs. studies of wild animals, we urge caution when comparing our results to previous studies involving captive primates and rodents.
9. Comment: The measurement of faecal glucocorticoid metabolite is a valid method. But stress levels are prone to outer environmental changes e.g. predatory attacks and human interactions increase stress levels. Accounting for these events is often difficult due to limited time windows of observation. "Behavioral data were collected daily between 0530 -0900 h and 1700 -2000 h." Therefore replicates of fGCMs for each animal are important as the period and day of sampling is critical. Please state if they were applied.
Response: Unfortunately, we do not know if hyenas experienced highly stressful events (e.g. encounter's with lions or Masai pastoralists) prior to defecating due to limited time windows of observation, and due to hyenas' fission-fusion social structure, which means that observers are never observing all the hyenas at once. However, we measured fGCMs, which represent stress levels over the course of hours, from multiple samples in order to calculate an average level of the stress hormone that minimizes the effect of any particular acute experience, such as predatory attacks or out-of-theordinary interactions with humans. Accordingly, we feel that our lack of data on such events does not pose a major issue.
Regarding variation in fGCMs that result from natural cycles in steroid hormones over the course of the day, we include time of day at which the fecal sample was collected in our models. This is noted in the results and the footnote of Table 1.
10. Comment: Global methylation analysis using LUMA only provides a methylation ratio over the entire genome. It does not provide the data for functional interpretation. The authors state "The majority of CpG sites assessed via LUMA occur in gene bodies, where they may function in transcription regulation and alternative splicing 27, as well as in non-coding regions of the genome, where they may repress repetitive elements 28 and enhance chromosome stability 29. Therefore, we assumed that lower than average %CCGG methylation is likely disadvantageous to health." This assumption is very farfetched, because LUMA does not provide the details to distinguish between functional important regions and mostly excludes promoters, which are thought to be the main regulatory regions. Therefore methylation changes are difficult to interpret and should be handled with caution. Using LUMA is costefficient and does not need a reference genome, which would be two valid reasons to use it. However the usage of three different approaches with different samples numbers: global DNA methylation (n = 186), genome-wide DNA methylation (n = 29), and candidate gene DNA methylation (n = 96) without explanation, is confusing. Please clarify why you chose this strategy.
Response: Thank you for raising these important points. First regarding the issue of interpreting our global DNA methylation results we kindly refer the reviewer to our responses to Reviewer 1, Comments #2 and 3 Regarding the use of multiple assays and differing sample sizes for the assays: each assay presents a tradeoff, some of which are logistical, including the lack of an available reference genome, and cost efficiency as noted by the reviewer. Other reasons for using multiple assays are biological and include interpretation of results that reflect broad scale genomic stability versus functional pathways vs potential candidate genes and their influence on phenotype. We included additional analyses and results that more explicitly highlight what each assay measures and provide an overview of the genomic context that is targeted by each assay. The revised text and new results and figures are in our response to Reviewer 1, Comment 7.
11. Comment: In order to assess the biological meaning of the established data, please clarify the numbers of animals, mothers, sex of offspring, and siblings, and the pairwise interactions observed. The number of samples collected for CD and DI time windows, stating why these were chosen, the sex, the cofounding factors that can lead to an increase in stress, the robustness of the metabolites measured. Last but not least, providing the data and code used is essential to understand and repeat the analysis.

Response:
Regarding "the numbers of animals, mothers, sex of offspring, and siblings, and the pairwise interactions observed": we report the requested information for each data set in the main text and provide additional information regarding sample sizes and the overlap between different data sets (i.e. the numbers of animals included in different models) in the supplemental information and table/figure legends.
Regarding "the number of samples collected for CD and DI time windows, stating why these were chosen, the sex, the cofounding factors that can lead to an increase in stress, the robustness of the metabolites measured": We have previously addressed this issue in our response to Reviewer 4 Comment #2.
Regarding "Last but not least, providing the data and code used is essential to understand and repeat the analysis": We will submit our code and related data sets to a public GitHub repository such that anyone can reproduce these analyses or implement a similar approach using their own data. We agree that making code available is important for reproducibility, and we have done this with previously published projects that use these data (c.f. the code repository for the analyses in Laubach et al.

<b>REVIEWER COMMENTS</B>
Reviewer #1 (Remarks to the Author): Overall, the authors succeed in addressing all but two of my concerns. I feel their revisions (as well as the strong original manuscript) will make an important contribution to our understanding of stress biology and early life adversity.
The following two points have not been adequately addressed: • I'm not convinced by the authors' response to concerns regarding cell type composition. The problem with changes in cell type composition isn't that they produce local effects at individual genes, but rather systemic biases. They therefore obscure true signal of specific genes or biological pathways, with signatures of cell-type differences. New citations to Schisterman and colleagues represent a minority and outdated opinion in the field.
If the authors still have access to the primary samples, I would suggest they use blood smears to estimate cell type composition for at least the EWAS analysis. If these samples are exhausted, I would remove the references to Schisterman et al. and leave the revised text as "Fourth, we used archived DNA from whole blood that lacked information on cell type composition, and therefore we were not able to control for cell type heterogeneity as a source of variation in DNA methylation. However, cell type composition may be influenced by factors like social stress (Engler et al., 2004), and therefore could be on the causal pathway between the explanatory variables of interest and DNA methylation." • Second, I repeat my suggestion that the authors consider systematic biases in the effect sizes across genomic contexts to add insight into the global methylation patterns. This was partially addressed for significantly differentiated sites [lines 453-455], but not for all sites which may reveal subtle, systematic effects that underlie the global methylation effects.
In response to my original comment, the authors stated "the present analysis is not an ideal setting in which to formally compare global vs genome-wide results for a few reasons. First, the overlap in the global and genome-wide DNA methylation samples includes only 12 animals, therefore a direct comparison between data sets is limited. Second, the purpose of the EWAS analysis was exploratory, and to identify potential functional pathways that can be assessed more rigorously in future analyses with larger sample sizes." The first point seems to be an advantage rather than a disadvantage. If the patterns observed in the larger, global DNA methylation dataset are robust we would expect them to be consistent with the genome-wide methylation dataset, even for completely independent datasets. Second, the exploratory nature of these EWAS results does not prevent them from being used for this purpose. In fact, directional biases in effect size and analyses by genomic context are in many ways less problematic than site-specific analyses using a small sample size. Finally, while I appreciate concerns regarding manuscript length, I think such an analysis could be a concise addition to the main text, such as: "Furthermore, our EWAS results reveal a bias in reduced methylation levels, consistent with our global analysis of DNA methylation. This effect is largely driven by reduced methylation in XXX genomic contexts, where XX contexts have no directional bias".
Reviewer #2 (Remarks to the Author): I thank the authors for their thoughtful revisions. I look forward to seeing the paper in print.
All the Best, Richard Hunter, Ph.D.
Reviewer #3 (Remarks to the Author): As conveyed in my initial review of this manuscript, I feel that it represents a valuable body of work that will be extremely informative as we reconcile laboratory and field studies targeting genetic and stress response perspectives. After reviewing the authors' thoughtful and thorough responses to the reviewer comments, I am satisfied that all points were appropriately considered and included in the revised manuscript.
Reviewer #4 (Remarks to the Author): Reviewers comment to Response: Thank you for the authors' additional work and replies to the reviewers` comments.
After re-reading the publication of Houslay & Wilson (2017. Behav. Eco. 28: 948-952) on the misuse of BLUPs, I remain doubtful that the method used by the authors is appropriate, as it hides uncertainty. The authors` response does not clarify the statistical or logical reason justifying that the problem highlighted by Houslay & Wilson doesn't apply, when BLUPs are used in subsequent analysis as explanatory variable. In the referred article I could not find any instance implying such a statement. Rather I read: "[…], it has become common practice to extract predictions of individual random effects from fitted mixed models and to use these in subsequent analyses, such as correlation tests or linear regression models (Table 1). Problems arise from this approach because individual point estimates from random effects in mixed models (sometimes known as conditional modes, or best linear unbiased predictors, BLUPs) are predicted with large amounts of error. Their use in secondary analyses can therefore lead to highly anticonservative tests of biological hypotheses, because the error inherent in their prediction is excluded from these further tests (Hadfield et al. 2010)." Houslay & Wilson further state "We stress that BLUP is an incredibly useful technique that should not be dismissed in any way as inherently "bad" (Robinson 1991). Indeed, it is entirely appropriate to use individual-level predictions to say something about individuals (or genotypes, or specific levels of some other random effect).", and "Nonetheless, use of these "stats on stats" approaches that are known to be inappropriate for hypothesis testing (see Brommer 2013b for further discussion) continues unabated." Using stats on stats should be avoided, as data might be highly fragile to changes e.g. after an increase in sample number. Thank you for adding the sample numbers, and including numbers of males and females. "Using blood samples from immobilized hyenas, we 74 constructed three primary datasets for analyses: global DNA methylation (n = 186 total; n = 99 75 females and n = 87 males), genome-wide DNA methylation (n = 29 total; n =29 females and n = 76 0 males), and candidate gene DNA methylation (n = 78 total; n = 43 females and n = 35 males)." With respect to these numbers, the experimental set-up seems not fully thought through beforehand and remains confusing. For example, for the analysis of DNA methylation three methods were applied using different individuals, amounts and different sex composition. Sex is, as pointed out before, an essential determiner within the hyena system, which as such has an effect on reproduction, behaviour and thus on DNA methylation. Furthermore, the usage of whole blood samples is tricky, blood-cell composition quickly changes e.g. in case of an infection numbers of leucocytes are increasing. The cell composition change is accompanied by a change in methylation. The authors cannot control for this -as stated in their response to reviewer 1.
Despite the authors` best effort, the results look fragile and dependent on the method used. I suggest applying a measure of sensitivity or even more holistic approach, such as a multiverse analysis proposed by Gelman (https://doi.org/10.1177/1745691616658637), to increase the study's significance.
While the idea and species` system is exciting, at this stage I am not convinced the results provide deeper knowledge into the species system.

Comment:
Overall, the authors succeed in addressing all but two of my concerns. I feel their revisions (as well as the strong original manuscript) will make an important contribution to our understanding of stress biology and early life adversity.

Response:
We are pleased that we have addressed the majority of this reviewer's concerns and look forward to making a formal contribution to the literature.
2. Comment: I'm not convinced by the authors' response to concerns regarding cell type composition. The problem with changes in cell type composition isn't that they produce local effects at individual genes, but rather systemic biases. They therefore obscure true signal of specific genes or biological pathways, with signatures of cell-type differences. New citations to Schisterman and colleagues represent a minority and outdated opinion in the field. If the authors still have access to the primary samples, I would suggest they use blood smears to estimate cell type composition for at least the EWAS analysis. If these samples are exhausted, I would remove the references to Schisterman et al. and leave the revised text as "Fourth, we used archived DNA from whole blood that lacked information on cell type composition, and therefore we were not able to control for cell type heterogeneity as a source of variation in DNA methylation. However, cell type composition may be influenced by factors like social stress (Engler et al., 2004), and therefore might be on the causal pathway between the explanatory variables of interest and DNA methylation." Response: Thank you. We have included the suggested language in the discussion of limitations.
Pg 39-40, Lines 761-765: Fourth, we used archived DNA from whole blood that lacked information on cell type composition, and therefore we were not able to control for cell type heterogeneity as a source of variation in DNA methylation. However, cell type composition may be influenced by factors like social stress 102 , and therefore might be on the causal pathway between the explanatory variables of interest and DNA methylation.
3. Comment: Second, I repeat my suggestion that the authors consider systematic biases in the effect sizes across genomic contexts to add insight into the global methylation patterns. This was partially addressed for significantly differentiated sites [lines 453-455], but not for all sites which may reveal subtle, systematic effects that underlie the global methylation effects. In response to my original comment, the authors stated "The present analysis is not an ideal setting in which to formally compare global vs genome-wide results for a few reasons. First, the overlap in the global and genome-wide DNA methylation samples includes only 12 animals, therefore a direct comparison between data sets is limited. Second, the purpose of the EWAS analysis was exploratory, and to identify potential functional pathways that can be assessed more rigorously in future analyses with larger sample sizes." The first point seems to be an advantage rather than a disadvantage. If the patterns observed in the larger, global DNA methylation dataset are robust we would expect them to be consistent with the genome-wide methylation dataset, even for completely independent datasets. Second, the exploratory nature of these EWAS results does not prevent them from being used for this purpose. In fact, directional biases in effect size and analyses by genomic context are in many ways less problematic than site-specific analyses using a small sample size. Finally, while I appreciate concerns regarding manuscript length, I think such an analysis could be a concise addition to the main text, such as: "Furthermore, our EWAS results reveal a bias in reduced methylation levels, consistent with our global analysis of DNA methylation. This effect is largely driven by reduced methylation in XXX genomic contexts, where XX contexts have no directional bias".
Response: Thanks to the reviewer for taking the time to clarify what they would like to see here. We now show patterns of positive and negative associations between DNA methylation and fecal corticosterone, comparing distributions of beta estimates within genomic regions vs across the whole genome (SI Figure 11; pasted below). We have also added text in the Methods & Results sections describing this figure.

Methods:
Pg 18, Lines 350-355: Finally, among the 25 individuals for which we assayed mERRBS, we annotated CpG sites from the EWAS in genomic regions of the draft hyena genome as described in the Supporting Material. We then assessed patterns of positive and negative associations of DNA methylation with fGCMs and compared the distribution of Beta estimates within genomic regions vs. across the entire genome to assess for bias in the direction of associations from the EWAS.

Results:
Pg 26, Lines 499-503: While we did not observe systemic bias in enrichment of positive or negative associations from our EWAS, within particular genomic regions, including CpG islands, promoters and introns there appeared to be fewer Beta estimates near zero than compared to all CpG sites across the whole genome (Supporting Figure 11). the CpG density annotation. b) Beta values from stress hormone EWAS based on genic annotation.

Supporting
In addition, we agree that in a scenario where maternal care is positively associated with global DNA methylation, then it would be interesting to examine systematic associations/biases within genomic contexts identified by the EWAS. In doing this, we expect that a comparison of counts of loci with positive vs. negative associations between maternal care and (site-specific) DNA methylation would reveal genomic contexts driving the association captured by the global DNA methylation assay. This approach would be informative if we were modeling the same predictor of global DNA methylation as in the EWAS analysis. However, this is not the case for the present study. In this analysis, we examined associations of multiple early-life social experiences and global DNA methylation as part of a broader analysis that sought to test global DNA methylation as a possible mediating mechanism linking early-life social experiences to stress physiology. Because we did not find any associations of global DNA methylation and stress physiology (fecal glucocorticoids), we did not continue the assessment of mediation. Instead, we implemented an outcome-based EWAS (i.e., an EWAS that identified hits with respect to stress physiology) to identify genomic sites/regions associated with fecal glucocorticoids, and then used a meet-in-the-middle approach to identify those that are also associated with maternal care as an indicator of the early-life social experience (NB: we selected maternal care based on prior knowledge and findings in this study population). Since we did not conduct an EWAS in which early-life social experience is the explanatory variable, a direct comparison will not be informative.