Neighborhood disadvantage is a chronic form of adversity that is often characterized by high rates of poverty, limited physical (e.g., green space) and built (e.g., grocery stores, pharmacies) resources, community violence, high exposure to toxicants, and low social cohesion (Jutte et al., Reference Jutte, Miller and Erickson2015; Wodtke et al., Reference Wodtke, Harding and Elwert2011). This form of adversity has been demonstrated to have a robust effect on long-term physical (e.g., cardiovascular disease, cancer, obesity) (Cubbin et al., Reference Cubbin, Sundquist, Ahlén, Johansson, Winkleby and Sundquist2006; Jutte et al., Reference Jutte, Miller and Erickson2015) and mental health (e.g., depression, substance use) (Diez Roux & Mair, Reference Diez Roux and Mair2010; Jutte et al., Reference Jutte, Miller and Erickson2015) outcomes. Indeed, neighborhoods have been demonstrated to be a robust predictor of life expectancy discrepancies, future health, and life chances (Evans et al., Reference Evans, Zimmerman, Woolf and Haley2012; Haley et al., Reference Haley, Zimmerman, Woolf and Evans2012; Jutte et al., Reference Jutte, Miller and Erickson2015; Lavizzo-Mourey, Reference Lavizzo-Moureyn.d.). The effects of neighborhood on health outcomes have also been shown to persist when accounting for individual deprivation and characteristics (Steptoe & Feldman, Reference Steptoe and Feldman2001). Even so, positive adjustment and competent functioning within the context of such adversity or resilience (Luthar et al., Reference Luthar, Cicchetti and Becker2000; Masten, Reference Masten2001), is quite common (40%–62% of exposed youth) (Luthar et al., Reference Luthar, Grossman and Small2015; Masten, Reference Masten2001; Vanderbilt-Adriance & Shaw, Reference Vanderbilt-Adriance and Shaw2008). Resilient youth thus provide a model of successful adaptation to adversity as understanding how environmental and biological factors may enable these positive outcomes is of great importance to informing prevention and intervention efforts for youth in disadvantaged neighborhood contexts.
While much of the early literature in the field conceptualized resilience as a static individual trait, contemporary work has explicitly reconceptualized resilience as a dynamic outcome that is influenced by the individual’s attributes, as well as their familial and community-level contexts (Luthar et al., Reference Luthar, Cicchetti and Becker2000; Masten, Reference Masten2001; Rutter, Reference Rutter2006). The extant empirical literature on resilience has, in turn, largely focused on socioecological factors, identifying several factors (e.g., parenting behavior) that promote or constrain resilience (Curtis & Cicchetti, Reference Curtis and Cicchetti2003). In recent decades, however, a growing number of studies have begun to examine the role of biological mechanisms in the development of resilience (Burt, Reference Burt2017; Curtis & Cicchetti, Reference Curtis and Cicchetti2003; Karatsoreos & McEwen, Reference Karatsoreos and McEwen2013; Luthar et al., Reference Luthar, Cicchetti and Becker2000; McEwen et al., Reference McEwen, Gray and Nasca2015; Panter-Brick & Leckman, Reference Panter-Brick and Leckman2013). Recent theoretical work (e.g., biopsychosocial model) (Feder et al., Reference Feder, Fred-Torres, Southwick and Charney2019) specifically highlights the transactional relationship between socioecological and biological influences on youth resilience. One potential mechanism undergirding these transactions relates to epigenetics and the biological embedding of stress via DNA methylation (e.g., silencing or activation of genes). Several epigenetic studies have found evidence of DNA methylation that results from environmental stressors, predicting outcomes ranging from stress response (Smith et al., Reference Smith, Zhao, Wang, Ratliff, Mukherjee, Kardia, Liu, Roux and Needham2017) to physical health (Notterman & Mitchell, Reference Notterman and Mitchell2015) and depression (Sun et al., Reference Sun, Kennedy and Nestler2013).
Given the growing literature examining the role of DNA methylation in response to stressors, it is somewhat surprising to note that the literature examining the role of DNA methylation in resilience to stressors remains scarce. Three published studies have examined DNA methylomic biomarkers of resilience in human samples (Milaniak et al., Reference Milaniak, Cecil, Barker, Relton, Gaunt, McArdle and Jaffee2017; Miller et al., Reference Miller, Shakespeare-Finch and Bruenig2020), two of which examined DNA methylation in only one or two specific gene regions. Milaniak and colleagues (Milaniak et al., Reference Milaniak, Cecil, Barker, Relton, Gaunt, McArdle and Jaffee2017) found that DNA methylation in the oxytocin receptor gene at birth predicted psychological resilience (i.e., a lack of conduct problems) to prenatal environmental stressors in middle childhood (N = 321). Similarly, Miller and colleagues (Miller et al., Reference Miller, Shakespeare-Finch and Bruenig2020) found that DNA methylation of sites located on the NR3C1 and FKBP5 genes predicted psychological resilience (i.e., measured using the Brief Resilience Scale) among emerging to middle-aged adults (N = 49). Although these studies begin to provide proof of concept for the idea that DNA methylation is a mechanism supporting resilience to adversity, they were notably limited by their focus on specific gene regions despite the availability of methylome-wide arrays. Indeed, Lu and colleagues (Lu et al., Reference Lu, Hsieh, Yang, Wang and Lin2023) appear to have conducted the only methylome-wide association study (MWAS) on psychological resilience (N = 78; discovery sample N = 16, validation sample N = 62). While their study identified three differentially methylated probes (DMPs), there are a number of critical limitations of this study. First and foremost, the authors failed to account for multiple testing; given the p-values reported, the DMPs they identified would likely not remain statistically significant if appropriate correction methods (e.g., false discovery rate; Benjamini & Hochberg, Reference Benjamini and Hochberg1995) were applied. In addition, the study did not control for blood cell-type proportions in their analyses, which is known to lead to inflated test statistics. Also, despite their conceptualization of resilience as adaptation in the context of exposure to substantial stress, adversity, or trauma, the authors do not report whether their sample was restricted to individuals with such exposures (Lu et al., Reference Lu, Hsieh, Yang, Wang and Lin2023), thus limiting the generalizability of their findings. Of note, all three of these studies are further limited by their focus on only a single form of resilience despite the multidimensional nature of resilience (i.e., individuals may be resilient in one area but not another). Also, Lu et al. (Reference Lu, Hsieh, Yang, Wang and Lin2023) and Miller et al. (Reference Miller, Shakespeare-Finch and Bruenig2020) employed notably small samples, calling into question the robustness and generalizability of their findings, particularly given that the MWAS conducted by Lu et al. (Reference Lu, Hsieh, Yang, Wang and Lin2023) requires large samples to adequately detect effects. Thus, there is a clear and compelling need for studies to examine DNA methylomic biomarkers of multiple dimensions of resilience (academic, social, psychological, overall) across the entire methylome in sizable samples exposed to moderate to severe levels of adversity.
That said, there are a handful of relevant empirical studies using animal models. For example, Weaver et al., (Weaver et al., Reference Weaver, Cervoni, Champagne, D’Alessio, Sharma, Seckl, Dymov, Szyf and Meaney2004), revealed that high levels of maternal care altered DNA methylation at the GR exon 17 promoter site (accompanied by negative effects on the stress response system) (Szyf et al., Reference Szyf, Weaver, Champagne, Diorio and Meaney2005; Weaver et al., Reference Weaver, Cervoni, Champagne, D’Alessio, Sharma, Seckl, Dymov, Szyf and Meaney2004) in the first week of life and persisting into adulthood. What’s more, DNA methylation at this site appeared to be directly programmed by maternal behavior and reversible through the use of a histone deacetylase inhibitor trichostatin A. Elliot and colleagues (Elliott et al.,, Reference Elliott, Ezra-Nevo, Regev, Neufeld-Cohen and Chen2010) also assessed changes in DNA methylation in rats exposed to a social defeat protocol (Krishnan et al., Reference Krishnan, Han, Graham, Berton, Renthal, Russo, LaPlant, Graham, Lutter, Lagace, Ghose, Reister, Tannous, Green, Neve, Chakravarty, Kumar, Eisch, Self, Lee, Tamminga, Cooper, Gershenfeld and Nestler2007) and found, that while most mice avoided their neighbor following the protocol, a subset of mice with significantly increased DNA methylation of the Crf promoter exhibited behavioral resiliency to the social defeat and interacted with the neighbor. These findings collectively bolster conclusions that both promotive and stressful life events may alter DNA methylation with downstream developmental consequences.
In sum, although research is still limited, there is reason to expect that DNA methylation may be an important component of resilience to adversity. Meaningfully extending this line of work to understand resilience in living humans is trickier than it might seem, however. Although usually discussed as a product of the environment only, DNA methylation is also genetically influenced (Grundberg et al., Reference Grundberg, Meduri, Sandling, Hedman, Keildson, Buil, Busche, Yuan, Nisbet, Sekowska, Wilk, Barrett, Small, Ge, Caron, Shin, Lathrop, Dermitzakis, McCarthy and Zondervan2013; Zhang et al., Reference Zhang, Moen, Liu, Mu, Gamazon, Delaney, Wing, Godley, Dolan and Zhang2014; van Dongen et al., Reference van Dongen, Nivard, Willemsen, Hottenga, Helmer, Dolan, Ehli, Davies, van Iterson, Breeze, Beck, Pool, van Greevenbroek, Stehouwer, Kallen, Schalkwijk, Wijmenga, Zhernakova, Tigchelaar, Beekman and Boomsma2016). As such, what may appear to be environmentally induced DNA methylation for a given outcome could in fact reflect genetic effects, a potential confound that undercuts the conclusions of human DNA methylation studies, including the previously discussed MWAS on psychological resilience (Lu et al., Reference Lu, Hsieh, Yang, Wang and Lin2023). Monozygotic (MZ) twin difference designs are considered the gold standard for overcoming this uncertainty in living humans (Burt et al., Reference Burt, McGue, Iacono and Krueger2006). MZ twins are genetically identical and yet can and do have different DNA methylomes as a result of their unique environmental experiences (Fraga et al., Reference Fraga, Ballestar, Paz, Ropero, Setien, Ballestar, Heine-Suñer, Cigudosa, Urioste, Benitez, Boix-Chornet, Sanchez-Aguilera, Ling, Carlsson, Poulsen, Vaag, Stephan, Spector, Wu, Plass and Esteller2005). Unfortunately, most twin studies are population-based and include relatively few youths exposed to adversity and even fewer who demonstrate resilience to that adversity. The utilization of a sample enriched for disadvantage to study the role of DNA methylation in resilience would thus offer significant promise for our understanding of differences in adaptability to adversity.
Current study
The current study aimed to identify DNA methylation biomarkers of resilience in a unique sample of twins enriched for disadvantage. We identified DNA methylation sites associated with academic resilience, social resilience, psychological resilience, and resilience across domains. Analyses were conducted using the entire sample of twins, allowing us to identify general methylomic biomarkers of resilience. We then conducted twin difference analyses of the significant and suggestive CpG sites in only MZ pairs, allowing us to narrow in on those sites that are specifically environmental in origin. We hypothesized that we would find evidence of methylated sites that are associated with resilience (i.e., academic, social, psychological, and across domains) to disadvantage and that differences in DNA methylation between MZ twins will predict differences in their resilience, strengthening causal inferences.
Methods
Participants
Participants were recruited as part of the Twin Study of Behavioral and Emotional Development in Children (TBED-C), a study within the population-based Michigan State University Twin Registry (Burt & Klump, Reference Burt and Klump2019). The TBED-C sample encompasses two arms of participants assessed between 2008 and 2015: a population-based arm of 1,054 twins from 528 families recruited from across lower Michigan and an under-resourced arm of 1,000 twins from 502 families residing in modestly to severely disadvantaged neighborhoods in the same recruitment area. Participating twins were screened for cognitive and physical conditions that would impede completion of the assessment (e.g., a significant developmental delay). Children provided informed assent, and informed consent was obtained from parents. Zygosity was determined using physical similarity questionnaires administered to the twins’ primary caregiver (Peeters et al., Reference Peeters, Van Gestel, Vlietinck, Derom and Derom1998).
Recruitment procedures are detailed at length in prior work (Burt & Klump, Reference Burt and Klump2019). In brief, families were recruited directly from birth records, or from a population-based registry that was itself recruited via birth records, via anonymous recruitment mailings in conjunction with the Michigan Department of Health and Human Services. Recruitment procedures for the under-resourced sample were restricted to those families residing in neighborhoods where neighborhood poverty was 10.5% (the median for Michigan neighborhoods in 2008) or greater, meaning that 10.5% or more of households were living below the poverty line according to census-level data. The response rate for the population-based and under-resourced arms of the sample was 62% and 57%, respectively. The under-resourced arm of the sample was significantly more racially diverse (15% Black, 75% White) than the population-based arm of the sample, reported lower family income (the means were $72,027 and $57,281, respectively; Cohen’s d = –0.38), and had higher paternal felony convictions (d = 0.30). The final under-resourced arm of the sample appears representative of the full sample of families we attempted to recruit as indexed via a brief questionnaire administered to approximately 85% of nonparticipating families (Burt & Klump, Reference Burt and Klump2019).
Participants in the current study represent a subsample of the under-resourced arm of the sample, as well a subsample of families from the population-based arm of the sample who would have met the criteria for the under-resourced arm (i.e., those living in neighborhoods with above median poverty). This totaled a possible sample of 768 families residing in disadvantaged neighborhood contexts, of which saliva assays have been completed for 240 twin pairs (the majority of whom were MZ pairs). Following assay quality control procedures and exclusion of participants with insufficient informant data to compute outcomes of interest, 270 participants from 135 full twin pairs (115 MZ; 20 dizygotic [DZ]) and six singletons (i.e., twins without a pair) formed the primary analytic sample for the current study (total N = 276 individuals). All 20 DZ pairs were male-male, whereas among MZ pairs, 69 were male-male, and 46 were female-female. The remaining singletons included five males and one female. All twins ranged in age from 6 to 11 years old at the time their questionnaires, and saliva samples were collected. The majority of participants in our final analytic sample identified as White (77.8%), 10.6% identified as Black, 2.1% as Native American, 2.1% as Pacific Islander, 1.4% as Latinx, and 6% identified as “Other” or a race prominent in less than 1% of the sample (i.e., Asian). Finally, the mean level of neighborhood poverty was 23%, while the mean family income was approximately $40,000 for a family of four.
Measures
As resilience is inherently a conditional construct – in that youth cannot demonstrate resilience without having first been exposed to adversity – it must be inferred through demonstrated competency and positive mental health despite exposure to adversity. In our case, we focused on resilience to moderate to severe neighborhood disadvantage, a form of chronic adversity. Competency and mental health were assessed via maternal reports on the Child Behavior Checklist (CBCL) (Achenbach & Rescorla, Reference Achenbach and Rescorla2001). The CBCL is one of the most commonly used and well-validated instruments for assessing academic and social competence, as well as mental health (internalizing and externalizing) problems prior to adulthood (Nakamura et al., Reference Nakamura, Ebesutani, Bernstein and Chorpita2009).
Academic resilience
The School Competency subscale of the CBCL served as our continuous measure of academic resilience (α = .64). This subscale includes items that assess school performance across subject domains, special education services received, repeated classes, and academic or other school-related problems (e.g., Does your child receive special education or remedial services or attend a special class or special school?). Mothers responded to a four-part question about academic performance on a 4-point scale ranging from “failing” to “above average,” as well as 3 binary (yes/no) questions. Of note, this score was kurtotic due to the narrow range of the subscale and was thus transformed by taking the natural log of each score to remove kurtosis prior to analyses.
Social resilience
The Social Competency subscale of the CBCL served as our continuous measure of social resilience (α = .49). Mothers responded to six questions assessing the child’s involvement in organizations, number of friends, contact with friends, behavior with others, and behavior alone (e.g., About how many times a week does your child do things with any friends outside of regular school hours?). Of note, the lower reliability evidenced in the school and social competence subscales is not uncommon given that they are multidimensional in nature such that most items capture different aspects of social and school competence.
Psychological resilience
An absence of psychopathology count variable served as our measure of psychological resilience (α = .78). Mothers rated the extent to which a series of statements described their child’s behavior during the past 6 months; responses were made on a 3-point scale ranging from 0 (never) to 2 (often/mostly true). We examined all eight psychopathology scales in the CBCL: anxious/depressed (e.g., fears certain animals, situations, or places, other than school), withdrawn/depressed (e.g., there is very little he/she enjoys), somatic complaints (e.g., constipated, doesn’t move bowels), social problems (e.g., complains of loneliness), thought problems (e.g., hears sounds or voices that aren’t there), attention problems (e.g., can’t concentrate, can’t pay attention for long), rule-breaking (e.g., breaks rules at home, school, or elsewhere), and aggressive behavior (e.g., destroys things belonging to his/her family or others). For the current study, we recoded each of these eight subscales as binary variables that indicate whether the child was at or above (0) or below (1) the CBCL’s empirically established borderline clinical significance cut point for that scale (Achenbach & Rescorla, Reference Achenbach and Rescorla2001). The eight dichotomous variables were then summed to form an absence of psychopathology score ranging from 0 to 8, where a higher score reflects less psychopathology and greater psychological resilience. Of note, this score was negatively skewed due to lower levels of psychopathology in our nonclinical sample and was thus transformed by taking the natural log of each score to reduce the skew prior to analyses.
Resilience across domains
Consistent with state-of-the-science studies of socio-emotional resilience, we are defining overarching resilience in the face of disadvantage as both the absence of psychopathology and the presence of social and academic competencies (Luthar et al., Reference Luthar, Cicchetti and Becker2000; Masten, Reference Masten2001; Rutter, Reference Rutter2006). Therefore, a dichotomous indicator of resilience across domains was computed with individuals above the CBCL social and academic competency subscale cut points (t-score = 40) (Achenbach & Rescorla, Reference Achenbach and Rescorla2001) and below the CBCL internalizing and externalizing score borderline cut points (t-score = 60) (Achenbach & Rescorla, Reference Achenbach and Rescorla2001) considered “resilient” (N = 135), whereas all others were considered “non-resilient” (N = 141) in at least one domain. Seventy-five twin pairs were concordant for resilience across domains, while 60 pairs were discordant for resilience across domains.
Assaying the methylome
Saliva samples were collected during the twin family’s assessment using Oragene collection kits (DNA Genotek). DNA was extracted using the Oragene Laboratory Protocol Manual Purification of DNA. Extracted DNA was then sodium bisulfite converted, and methylation was assessed in the converted DNA using the Infinium Human Methylation EPIC Bead Chip (Illumina). DNA conversion and methylation measurement were performed by the University of Michigan Sequencing Core.
Thorough quality control and intra-sample normalization procedures were employed using the Chip Analysis Methylation Pipeline for Illumina HumanMethylation450 and EPIC (ChAMP) Bioconductor package (Butcher & Beck, Reference Butcher and Beck2015; Morris et al., Reference Morris, Butcher, Feber, Teschendorff, Chakravarthy, Wojdacz and Beck2014). Samples with a high proportion of failed probes (≥10%) were removed (n = 1). Probes were removed if their detection p-value was above 0.01 (n = 86,415 probes), if the bead count was less than 3 in at least 5% of samples (n = 3608 probes), if probes aligned to multiple locations (cross-hybridizing probes) (Nordlund et al., Reference Nordlund, Bäcklin, Wahlberg, Busche, Berglund, Eloranta, Flaegstad, Forestier, Frost, Harila-Saari, Heyman, Jónsson, Larsson, Palle, Rönnblom, Schmiegelow, Sinnett, Söderhäll, Pastinen, Gustafsson, Lönnerholm and Syvänen2013), if probes were not located at CpG sites (n = 2242), if probes overlapped with single nucleotide polymorphisms, or if probes were located on sex chromosomes (n = 12,610). In order to correct for probe design bias, we used the champ.norm function (Teschendorff et al., Reference Teschendorff, Marabita, Lechner, Bartlett, Tegner, Gomez-Cabrero and Beck2013) of the ChAMP package. The ComBat function of the Surrogate Variable Analysis Bioconductor package was then used to correct for batch effects by slide and then array (Leek, Reference Leek2020). Finally, cell-type proportions were estimated for the most common cell types in saliva using the Epigenetic Dissection of Intra-Sample-Heterogeneity (EpiDISH) Bioconductor package (Zheng et al., Reference Zheng, Breeze, Beck and Teschendorff2018). These procedures yielded DNA methylation values (log2 methylated/unmethylated DNA at a specific probe, i.e., M-values) across 728,396 CpG sites for 276 participants.
Methylome-wide association study (MWAS)
The MWAS was performed on the full sample (N = 276) using regression to identify DNA methylation sites that were associated with resilience (i.e., social, academic, psychological, and across domains), so-called DMPs. Specifically, we fit logistic and ordinary least squares regression models in R for our dichotomous (i.e., resilience across domains) and continuous/discrete (i.e., social, academic, and psychological resilience) outcomes, respectively. To account for the nonindependence of twins within pairs, we corrected the standard errors by fitting our models within a heteroskedasticity-consistent covariance matrix estimator using the sandwich package in R (Zeileis, Reference Zeileis2006). To control for potential confounders, we included sex, age, zygosity, ethnicity, and estimated cell-type proportions as covariates in our models. A p-value threshold of P < 9 x 10-8 was used to declare a DMP methylome-wide significant (Mansell et al., Reference Mansell, Gorrie-Stone, Bao, Kumari, Schalkwyk, Mill and Hannon2019) and P < 1 x 10−5 for suggestive DMPs (Lander & Kruglyak, Reference Lander and Kruglyak1995).
Pathway analysis
To gain insight into the biological pathways affected by resilience, we used ConsensusPathDB (Kamburov et al., Reference Jutte, Miller and Erickson2009, Reference Kamburov, Pentchev, Galicka, Wierling, Lehrach and Herwig2011) to test for overrepresentation of top suggestive MWAS findings located within genes in the biological pathways in the Reactome (Croft et al., Reference Croft, Mundo, Haw, Milacic, Weiser, Wu, Caudy, Garapati, Gillespie, Kamdar, Jassal, Jupe, Matthews, May, Palatnik, Rothfels, Shamovsky, Song, Williams, Birney, Hermjakob, Stein and D’Eustachio2014) database. For a pathway to be considered enriched, a cut point of P < 0.01 was utilized, and at least two genes among the top MWAS findings had to be present in the pathway.
Monozygotic twin difference analyses
Finally, we performed twin difference analyses in R in which we compared MZ co-twins to strengthen causal inferences. Because MZ co-twins cannot differ in their epigenome as a consequence of genetic differences (as they are genetically identical), any differences in the methylome of co-twins point toward environmental mediation. We computed differences in DNA methylation scores for the significant and suggestive DMPs from the MWASs as well as for the four resilience phenotypes. For our twin difference analyses, we regressed DNA methylation difference scores for the DMPs and covariates (i.e., sex, age, and ethnicity, each on the twin-pair level) on resilience (i.e., academic, social, psychological, and across domains) difference scores. DMPs were then compared to a 95% statistical significance threshold (p ≤ 0.05).
Results
Descriptive statistics
Descriptive statistics for resilience across domains, psychological resilience, academic resilience, and social resilience are available in Table 1. While scores for psychological, academic, and social resilience were continuous, the score for general resilience across domains was dichotomous. Approximately half of the participants were considered to be resilient across domains. The majority of participants exhibited high scores for psychological and academic resilience; however, social resilience scores were more variable. Finally, the means and standard deviations of the four resilience phenotypes in MZ twins and DZ twins were equivalent. Pearson correlations between all cell type proportion estimates and the four resilience phenotypes were first examined; none were significant (Table S1).
Table 1. Descriptive statistics
Note. On the left are the descriptive statistics across individuals who are in a monozygotic twin pair, and on the right are the descriptive statistics across individuals who are in a dizygotic twin pair. Means, standard deviations (SD), minimums (Min), maximums (Max), and sample size (N) are presented for each of the four resilience phenotypes and DNA methylation estimates (DNAm). The proportion of discordant co-twins (disc.; i.e., nonidentical scores), co-twin mean difference scores (Tw diff. mean), and the standard deviation for co-twin mean difference scores (Tw diff. SD) are also presented for monozygotic twins.
A large proportion of co-twins differed in their degree of resilience, as indexed dimensionally (Table 1). Most co-twins (71%) had different levels of social resilience, with a mean co-twin difference of 43% of the typical phenotypic standard deviation across the full sample. Roughly half of the co-twins (45%) had different levels of academic resilience, with a mean co-twin difference of 38% of the typical phenotypic standard deviation across the full sample. Finally, although only a third (36%) of co-twins evidenced different levels of psychological resilience, those that differed did so quite a bit, with a mean co-twin difference that was 59% of the typical phenotypic standard deviation across the sample. For our dichotomous phenotype of resilience across domains, 46% of co-twins were discordant.
As a final point, we note that no co-twins across the entire sample had identical DNAm scores for any of the 728,396 CpG sites. The observed differences were quite large. The mean co-twin difference was 300% of the typical DNAm standard deviation across the full sample. What’s more, even when concordance was evaluated somewhat liberally (i.e., a co-twin difference of .001 or less; DNAm range: 0–1), twin pairs remained discordant on 94-98% of the CpG sites. When discordance was evaluated quite liberally (i.e., a co-twin difference of .01 or less), twin pairs were still discordant on 52-80% of the CpG sites.
Methylome-wide association study (MWAS)
The quantile-quantile plots for each of the resilience outcomes are shown in Figure 1. The number of points above the 95% confidence interval, deviating from the line of expected points according to the null hypothesis, indicates a considerable number of statistically significant or suggestive findings for resilience across domains, academic resilience, and social resilience. However, the plot for psychological resilience does not depict points above the 95% confidence interval, suggesting limited significant results for this outcome.
Figure 1. Quantile-quantile plots for MWAS of each resilience domain. Note. The observed p-values (black open circles), on a -log10 scale, are plotted against their expected values (red main diagonal line) under the null hypothesis assuming none of the CpGs have an effect. Shaded grey bands indicate the 95% confidence bands (CI).
The top 10 significant (P < 9 x 10−8; Mansell et al., Reference Mansell, Gorrie-Stone, Bao, Kumari, Schalkwyk, Mill and Hannon2019) and/or suggestive (P < 1 x 10−5; Lander & Kruglyak, Reference Lander and Kruglyak1995) MWAS DMPs and test statistics for each outcome are provided in Table 2 with covariate results for these DMPs provided in Table S2 and full results available in Table S3. Results indicated that, although there were no methylome-wide significant DMPs associated with resilience across domains, there were 90 suggestive DMPs. One of the top suggestive DMPs was located in an intron of SOX30, which is a member of the SOX family of transcription factors involved in determining cell fate and regulating embryonic development (Osaki et al., Reference Osaki, Nishina, Inazawa, Copeland, Gilbert, Jenkins, Ohsugi, Tezuka, Yoshida and Semba1999).
Table 2. Top 10 significant and/or suggestive differentially methylated probes
Note. ‘Probe’ is the name of the CpG probe in the human reference genome hg19/GRCh37, ‘Chr’ is chromosome, ‘start’ is the base pair location of the probe, ‘gene’ is the gene the probe is located in, and ‘genomic feature’ indicates if the probe is located in an intron, exon, or CpG island. Also shown are the signed test statistic values for regression: ‘Z-value’ for the dichotomous outcome of resilience across domains, ‘T-value’ for the continuous outcomes, ‘P-values,’ and ‘beta’ or regression coefficient. The top 10 methylome-wide significant (P - value ≤ 9 × 10−8) and/or suggestive (P - value ≤ 1 × 10−5) MWAS DMPs are displayed for each outcome.
The psychological resilience MWAS yielded no methylome-wide DMPs, but two suggestive ones. The top suggestive DMP was located in an intron of HOXC13, which has been implicated in cancer prognosis and belongs to the homeobox family of genes that encode transcription factors involved in morphogenesis (Panagopoulos et al., Reference Panagopoulos, Isaksson, Billström, Strömbeck, Mitelman and Johansson2003).
There were two methylome-wide significant and 20 suggestive DMPs associated with academic resilience. The top methylome-wide significant DMP was located in an intron of MYO10, which encodes a member of the myosin superfamily proteins and is associated with an increased risk for childhood apraxia of speech (Peter et al., Reference Peter, Wijsman, Nato, Matsushita, Chapman, Stanaway, Wolff, Oda, Gabo and Raskind2016). A top suggestive DMP was located in an intron of BRF1, which encodes a subunit of the RNA polymerase III transcription initiation factor and has been associated with neurodevelopmental abnormalities (Borck et al., Reference Borck, Hög, Dentici, Tan, Sowada, Medeira, Gueneau, Thiele, Kousi, Lepri, Wenzeck, Blumenthal, Radicioni, Schwarzenberg, Mandriani, Fischetto, Morris-Rosendahl, Altmüller, Reymond, Nürnberg, Merla, Dallapiccola, Katsanis, Cramer and Kubisch2015).
Finally, there were six methylome-wide significant and 54 suggestive DMPs associated with social resilience. The top methylome-wide significant DMP was located in an intron and CpG island of AC006372.5, also known as LOC101927914, an uncharacterized RNA gene. The second top methylome-wide significant DMP, as well as a suggestive DMP, was located in an intron of HLA-DRB1. In addition, another suggestive DMP was located in an intron of HLA-DQB2. HLA-DRB1 and HLA-DQB2 are located in the HLA region on chromosome 6, a large region of linkage disequilibrium indicating that these may not be independent signals (Simmonds & Gough, Reference Simmonds and Gough2007).
Sensitivity analyses were also performed on the normalized data with slide and array entered as covariates to verify that the ComBat correction we employed was not inflating results (Zindler et al., Reference Zindler, Frieling, Neyazi, Bleich and Friedel2020); this alternative approach yielded a greater number of significant and suggestive probes, suggesting that our approach was more conservative.
Enriched pathways
The majority of significant or suggestive DMPs were located in unique genes; 76 of 90 for resilience across domains, 2 of 2 for psychological resilience, 16 of 22 for academic resilience, and 47 of 60 for social resilience. All of the significantly enriched pathways are provided in Table 3. Resilience across domains yielded four significantly enriched pathways. The top significant pathway was the “Listeria Monocytogenes Entry into Host Cells,” which is involved in regulating the entry of bacterium that cause the majority of foodborne outbreaks. No prominent theme emerged among these results. There were no significant enriched pathways for psychological resilience, likely due to the small number of significant or suggestive DMPs for this outcome.
Table 3. Enriched pathways
Note. ‘Pathway’ is the name of the significantly enriched pathway from the Reactome database, ‘effective size’ is the number of genes involved in the corresponding pathway, and ‘gene overlap’ provides the names of genes from the MWAS that are present in the pathway. Also shown are the signed test statistic values for the pathway analyses: ‘p-value’ and ‘q-value.’
For academic resilience, we observed eight significantly enriched pathways. The POLR2L and BRF1 genes were found in five pathways implicated in the transcription or initiation of RNA polymerase III. RNA polymerase III serves as a catalyst for the synthesis of small RNAs (e.g., tRNAs, 5S rRNA, snRNA) considered to be essential for various cellular functions (Abascal-Palacios et al., Reference Abascal-Palacios, Ramsay, Beuron, Morris and Vannini2018). The POLR2L gene encodes a subunit of RNA polymerase I, II, and III and is therefore heavily involved in synthesizing messenger RNAs (Acker et al., Reference Acker, Murroni, Mattei, Kedinger and Vigneron1996). In addition, the POLR2L and LIG3 genes were found in three pathways involved in gap-filling and nucleotide excision DNA repairs. As a member of the DNA ligase family, the LIG3 gene is involved in excision repairs and has been linked to increased risk for cancer (Li et al., Reference Li, Suzuki, Liu, Morris, Liu, Okazaki, Li, Chang and Abbruzzese2009; Li et al., Reference Li, Wang, Wang, Guan, Guo, Wang, Zhang, Niu, Zhang, Wang and Yang2018), neural tube defects (Li et al., Reference Li, Wang, Wang, Guan, Guo, Wang, Zhang, Niu, Zhang, Wang and Yang2018), Alzheimer disease (Kwiatkowski et al., Reference Kwiatkowski, Czarny, Toma, Korycinska, Sowinska, Galecki, Bachurska, Bielecka-Kowalska, Szemraj, Maes and Sliwinski2016), and recurrent depression (Czarny et al., Reference Czarny, Kwiatkowski, Toma, Kubiak, Sliwinska, Talarowska, Szemraj, Maes, Galecki and Sliwinski2017).
Social resilience evidenced nine significantly enriched pathways. The HLA-DRB1 and HLA-DQB2 genes appeared in eight of these pathways, most of which are involved in T-cell receptor signaling, indicating that these results were driven by the HLA region on chromosome 6. The HLA region includes several genes – such as the HLA-DRB1 and HLA-DQB2 genes – that play a central role in immune system functioning (Simmonds & Gough, Reference Simmonds and Gough2007). The HLA region is associated with longevity (Joshi et al., Reference Joshi, Pirastu, Kentistou, Fischer, Hofer, Schraut, Clark, Nutile, Barnes, Timmers, Shen, Gandin, McDaid, Hansen, Gordon, Giulianini, Boutin, Abdellaoui, Zhao and Wilson2017), cognitive ability (Payton et al., Reference Payton, Van Den Boogerd, Davidson, Gibbons, Ollier, Rabbitt, Worthington, Horan and Pendleton2006), and mental health disorders (e.g., schizophrenia, autism) (Bennabi et al., Reference Bennabi, Gaman, Delorme, Boukouaci, Manier, Scheid, Si Mohammed, Bengoufa, Charron, Krishnamoorthy, Leboyer and Tamouza2018; Halley et al., Reference Halley, Doherty, Megson, McNamara, Gadja and Wei2013).
Monozygotic twin differences
For our final analyses, we sought to evaluate the extent to which the significant and suggestive DMPs from each of the MWAS models above were environmental in origin via MZ twin differences analyses. Results are provided in Table 4. Two DMPs for resilience across domains differed significantly across MZ pairs. The top DMP was located in Y_RNA, a class of small non-encoding RNAs involved in the repression of Ro60 (i.e., a protein that has been implicated in responses to environmental stress) as well as the initiation of chromosomal DNA replication (Christov et al., Reference Christov, Gardiner, Szüts and Krude2006). The second top DMP was located in an intron of TMEM67, a gene needed to facilitate ciliary structure and function (Yinsheng et al., Reference Yinsheng, Miyoshi, Qin, Fujiwara, Yoshimura and Katayama2022); defects can cause Joubert syndrome (characterized by abnormal brain development) and Meckel syndrome (most commonly characterized by enlarged kidneys).
Table 4. Significant monozygotic twin difference differentially methylated probes
Note. ‘Probe’ is the name of the probe in the human reference genome hg19/GRCh37, ‘Chr’ is chromosome, ‘start’ is the base pair location of the probe, ‘gene’ is the gene the probe is located in, and ‘genomic feature’ indicates if the probe overlaps with introns, exons, or CpG islands. Also shown are the signed test statistic values for regression: ‘Z-value’ for the dichotomous outcome of resilience across domains, ‘T-value’ for the continuous outcomes, ‘P-values,’ and ‘eta’ or regression coefficient. Finally, beta values for all covariates (i.e., ‘gender beta,’ ‘age beta,’ and ‘race beta’) are provided; no covariates’ beta values were significant at p < .05. All significant (P ≤ .05) DMPs are provided for each of the outcomes.
Four DMPs for social resilience also differed significantly across MZ pairs. The top DMP was located in an intron of LINC01250, a long intergenic non-protein coding RNA gene (Dungan et al., Reference Dungan, Qin, Hurdle, Haynes, Hauser and Kraus2021). The second top DMP was located in an intron of PIGG, a protein-coding gene involved in glycosylphosphatidylinositol-anchor biosynthesis; allelic variants of PIGG have been linked to intellectual disability with hypotonia and seizures (Makrythanasis et al., Reference Makrythanasis, Kato, Zaki, Saitsu, Nakamura, Santoni, Miyatake, Nakashima, Issa, Guipponi, Letourneau, Logan, Roberts, Parry, Johnson, Matsumoto, Hamamy, Sheridan, Kinoshita, Antonarakis and Murakami2016).
For academic resilience, two DMPs differed significantly across MZ pairs. The top DMP was located in an intron of RABEP2, a protein-coding gene that enables GTPase activator activity and growth factor activity (Kofler et al., Reference Kofler, Corti, Rivera-Molina, Deng, Toomre and Simons2018). The second top DMP was located in an intron and CpG Island of the aforementioned BRF1. DMPs for psychological resilience did not differ across MZ pairs.
Discussion
The goal of this study was to identify epigenetic correlates of resilience to neighborhood disadvantage in a sample of living humans. MWAS analyses conducted in 276 twins within 141 families revealed a handful of methylome-wide significant DMPs associated with academic as well as social resilience, and suggestive DMPs associated with each of the four resilience phenotypes examined (i.e., psychological, academic, social, and across domains). Pathway analyses revealed significantly enriched pathways for academic and social resilience, as well as resilience across domains. Results for academic resilience to neighborhood disadvantage pointed to DNA methylation in pathways related to DNA repair as well as the transcription and initiation of RNA polymerase III. DNA damage typically triggers a response that includes DNA repair. Dysregulation of DNA damage responses can result in developmental and neurological defects (Lee et al., Reference Lee, Choi, Kim and Kim2016). As mentioned previously, RNA polymerase III is involved in transcribing small RNAs. Dysregulation of small RNAs is thought to be implicated in abnormal brain development (Chang et al., Reference Chang, Wen, Chen and Jin2009). Taken together, these findings suggest that DNA methylation in these two pathways may alter or inhibit the regulation of DNA damage responses and small RNAs.
These enriched pathways also highlight the role of DNA methylation of the BRF1 gene in academic resilience. Mutations in BRF1 have been shown to cause central nervous system and neurodevelopmental anomalies due to a reduction in protein activity. It has been suggested that RNA polymerase III transcription initiated by BRF1 is necessary for typical cognitive development (Borck et al., Reference Borck, Hög, Dentici, Tan, Sowada, Medeira, Gueneau, Thiele, Kousi, Lepri, Wenzeck, Blumenthal, Radicioni, Schwarzenberg, Mandriani, Fischetto, Morris-Rosendahl, Altmüller, Reymond, Nürnberg, Merla, Dallapiccola, Katsanis, Cramer and Kubisch2015), a process that may be affected by DNA methylation of BRF1. The current study extends this line of work by demonstrating that an increase in DNA methylation of BRF1 is associated with academic resilience, a construct that is thought to be correlated with cognitive ability (Mayes et al., Reference Mayes, Calhoun, Bixler and Zimmerman2009; Tiet et al., Reference Tiet, Bird, Davies, Hoven, Cohen, Jensen and Goodman1998).
Results also suggest that DNA methylation in genes located in the HLA region involved in T-cell receptor (TCR) signaling may play a role in social resilience to neighborhood disadvantage. TCR signaling refers to cellular signaling cascades involved in determining cell fate, including cell survival, differentiation, and proliferation. TCRs typically bind to proteins involved in the immune response. Recent studies have demonstrated that proteins involved in the immune response are expressed in the central nervous system and play critical roles in synaptic transmission and plasticity, as well as refinement of connections during brain development (Garay & McAllister, Reference Garay and McAllister2010). Thus, DNA methylation of genes involved in TCR signaling may have downstream effects on brain development. Research on social cognition has demonstrated that the portions of the temporal lobe, the amygdala, and the cingulate cortex are implicated in social behavior via their involvement in the perception of social stimuli and the ability to link these stimuli to emotion, motivation, and cognition (Adolphs, Reference Adolphs2001). Therefore, while additional research is needed to confirm that TCR signaling impacts these brain regions in particular, this may explain its relationship with interpersonal functioning and social resilience (Cook et al., Reference Cook, Greenberg and Kusche1994).
MZ twin difference analyses revealed two suggestive DMPs for resilience across domains, two for academic resilience, and four for social resilience. While none of the significant methylome-wide DMPs differed across MZ twins, the suggestive DMPs that did differ across MZ pairs were located in genes implicated in responses to environmental stress and neurodevelopmental abnormalities.
Since MZ twins are genetically identical, significant findings point toward environmentally engendered DNA methylation in those cases. Alternatively, the absence of significant MZ differences in our methylome-wide significant DMPs suggests that those DMPs may not reflect causal environmental processes per se. This suggests that while some DMPs appear to be environmental in origin, others point to the possibility of genetic or developmental mediation of those methylomic effects. However, null results may also reflect family-wide influences or MZ differences that were too small to capture environmental mediation.
Limitations
The unique twin design of this study coupled with the relatively high degree of disadvantage experienced by participants uniquely positioned us to detect DMPs for resilience that are environmental in origin. However, there are limitations of the current study that are important to consider. First, because DNA methylation can be tissue-specific, etiological interpretations of saliva-based DNA methylation must be made with caution, the minimum interpretation being that DMPs could potentially be biomarkers of resilience. Given that our study did not contain a replication sample, our results are provisional and warrant further investigation with an independent sample. That being said, we are not aware of a second child twin sample experiencing sufficiently high rates of adversity to study resilience to that adversity at this time. Thus, diverse twin samples with higher rates of adverse exposures are needed to facilitate future replication. Next, although our sample is representative of racial and ethnic demographics throughout the state of Michigan, the racial breakdown of the sample is still primarily White, thereby limiting the generalizability of our findings to racially minoritized communities. It would be critical for future methylomic studies of resilience to recruit racially diverse samples. Lastly, while this study focuses specifically on resilience to neighborhood disadvantage, other forms of resilience may have distinct methylomic markers (e.g., resilience to maltreatment). Additional research on other forms of resilience would facilitate a comparison of methylomic markers across distinct forms of resilience.
Implications
Overall, this is one of the first studies to uncover potential DNA methylomic biomarkers of resilience in a sample of living humans. Our findings preliminarily highlight DNA methylation as a potential biological mechanism implicated in resilient outcomes, in that we identified a handful of methylome-wide significant and suggestive DNA methylation sites that predict resilience to neighborhood disadvantage. The etiologic inferences we can make about these DMPs and genes are more limited, however, since the significant DMPs from the MWAS did not differ across MZ twins. Such results support the possibility of genetic or developmental mediation for those DMPs. That said, we did identify a handful of suggestive methylomic correlates of resilience that differed across MZ twins. These environmental changes in the methylome are also at least nominally consistent with the biopsychosocial model’s theory in that they point to the importance of environmental effects, as well as reciprocal feedback between biology and the environment. Although beyond the scope of the current study, we intend to expand on our findings by examining the effect of specific environmental promotive factors (e.g., parental warmth) on DNAm sites implicated in youth resilience in the near future.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0954579424001330.
Funding statement
This study was supported by the National Institute of Mental Health (R01-MH081813 awarded to S.A.B. and K.L.K.), the Eunice Kennedy Shriver National Institute of Child Health and Human Development (R01-HD066040 awarded to S.A.B. and K.L.K.; R01-HD104297 awarded to S.A.B. and S.L.C.; F31HD111273 awarded to A.Y.V.), the National Institute of Health Office of the Director (UH3MH114249 awarded to S.A.B. and L.W.H.), and the National Science Foundation (GRFP awarded to A.Y.V.). The funding sources had no role in the study design, writing of the report, or decision to submit the article for publication.
Competing interests
None.
Open access
