Stress and Spirituality in Relation to HPA Axis Gene Methylation Among US Black women: Results from the Black Women’s Health Study and the Study on Stress, Spirituality and Health

Abstract

Background: Few epigenetics studies have been conducted within the Black community to examine the impact of diverse psychosocial stressors and resources for resiliency on the stress pathway (hypothalamus–pituitary–adrenal axis). Methods: Among 1000 participants from the Black Women’s Health Study, associations between ten psychosocial stressors and DNA methylation (DNAm) of four stress-related genes (NR3C1, HSDB1, HSD11B2 and FKBP5) were tested. Whether religiosity or spirituality (R/S) significantly modified these stress-DNAm associations was also assessed. Results: Associations were found for several stressors with DNAm of individual CpG loci and average DNAm levels across each gene, but no associations remained significant after false discovery rate (FDR) correction. Several R/S variables appeared to modify the relationship between two stressors and DNAm, but no identified interaction remained significant after FDR correction. Conclusion: There is limited evidence for a strong signal between stress and DNAm of hypothalamus–pituitary–adrenal axis genes in this general population cohort of US Black women.

Keywords::

• cohort study
• DNA methylation
• psychosocial stress
• religion
• resilience

Dramatic, unacceptable differences in the burden of chronic disease across ethnic communities continue to persist [1] despite decades of national efforts focused on reducing disparities [2]. In the US context, health disparities are often most dramatic with respect to the Black community [3]. Black populations experience a mortality rate from heart disease that is 23% higher than non-Hispanic whites, and also have higher rates of hypertension, obesity and diabetes compared with non-Hispanic whites [4]. Such disparities in the burden of illness will never be eliminated without addressing the psychosocial, environmental and economic conditions in which such disparities are rooted. Compared to their white counterparts, Black populations experience higher levels of inadequate housing [5], violent victimization [6], discrimination [7], childhood physical and sexual abuse [8], shift work [9] and poverty [10]. Numerous studies demonstrated the associations between psychosocial stress and increased risk of chronic disease, including type 2 diabetes and obesity [11,12], hypertension and cardiovascular disease [13] and cancer [14]. Stress can also influence disease risk indirectly, by increasing health-risk behaviors; stressed individuals are more likely to smoke, to sleep less, have lower-quality sleep and exercise less frequently [15–17] compared with those with lower levels of stress.

Black communities also engage in higher levels of community support activities, particularly in the realm of religion and spirituality (R/S). According to the most recent national data surveying over 35,000 Americans, more Blacks rate their religion as very important to them, attend religious services at least once a week, pray daily, attend weekly prayer groups and experience a feeling of spiritual well-being relative to whites [18]. Black communities also have more well-developed social and community relationships with their religious communities than whites [19], which has been shown to be inversely associated with mortality [20] and positively associated with health-promoting behaviors (fruit and vegetable consumption, moderate physical activity and decreased alcohol use) [21]. R/S have been increasingly recognized as important resources of resilience that can support health [22,23]. Measures of R/S have been prospectively associated with several mental health outcomes, including reduced risk of depression [24,25], anxiety or emotional distress [26] and reduced risk of suicide attempts [27,28]. Prospective analyses of chronic disease risk have reported associations between various measures of R/S and lower blood pressure and reduced risk of hypertension [29,30], cardiovascular events [31], obesity [32], mortality [33–35] and better self-rated health [36–39]. R/S practices, beliefs and experiences may be important for improving the health of Black communities and reducing health disparities due to psychosocial and socioeconomic adversity but have not been fully explored in epidemiological research to date. However, R/S can also have deleterious effects on health. For example, negative religious coping (also known as religious struggles) – in which individuals interpret stressful life situations as proof that they have been abandoned by their God, are being punished for their transgressions or come to question their faith altogether – has been associated with greater risk of mortality among medically ill elderly patients [40], and declines in immune status as measured by CD4 [41] and IL-6 levels [42].

Epigenomics is a potentially powerful tool for examining the ways in which psychosocial stress and social disadvantage at the individual and community levels become embodied [43] to increase the risk of many chronic illnesses. One possible biological pathway through which psychosocial stress affects disease etiology is dysregulation of the hypothalamus–pituitary–adrenal (HPA) axis, which is a major neuroendocrine system that controls reactions to stress. Numerous studies have identified DNA methylation (DNAm) of HPA axis genes as a potentially important epigenetic mechanism through which exposure to stressful physical and psychosocial environments across the lifespan may alter glucocorticoid regulation [44–51], which then increases the risk of multiple chronic diseases [52–63].

To date, with the exception of a few studies investigating the impact of psychosocial stress or resilience on DNAm of the NR3C1 [47, 64–68] and the FKBP5 gene [48,51,69–73], no other DNAm analyses focused on HPA axis genes have included significant numbers of Black participants. Furthermore, these extant studies have only assessed a single stress exposure (usually childhood maltreatment) or several conceptually similar stress exposures together. No study to date has systematically compared the relative influence of a broad array of psychosocial stressors on DNAm of multiple HPA axis genes in a Black population. As emerging platforms continue to require the pooling of data across multiple cohort studies, a key unanswered question is which psychosocial measures are most important for prospective cohorts to collect – both in terms of psychosocial stressors and resources for resilience – in building out the ‘social environment’ in environmental epigenomics research. Answering this question will require systematic and comparative research on the relative associations of different stress measures on the same outcome.

Despite researchers having used Black churches as recruitment sites for studies of conditions important to the Black community for many decades [74–76], few studies have assessed how dimensions of R/S within the Black community affect biological mechanisms that influence disease risk. With respect to epigenetic mechanisms, to date only five studies have assessed the associations of R/S with DNAm, all of which assessed meditation/mindfulness [77–79] and yoga [80], and none of which included Black participants. Furthermore, no other types of R/S experience beyond meditation and yoga have been studied in relation to DNAm in any ethnic community.

In this paper, these gaps in the literature are addressed prospectively and retrospectively by assessing the influence of ten measures of psychosocial stress in childhood and adulthood that disproportionately affect Black women on DNAm of four prominent HPA axis genes; and further by assessing the modifying influence of various dimensions of R/S on these relationships. The authors hypothesized that each psychosocial stressor would affect levels of DNAm and that R/S practices, beliefs and experiences would buffer the adverse physiological impact of stress (via DNAm of HPA axis genes), with the exception of measures of religious struggles, which were hypothesized to exacerbate the adverse impact of psychosocial stressors on DNAm.

Methods

Ethics statement

This study was approved by the Partners Human Research Committee and the Boston University Medical Center Institutional Review Board.

Black Women’s Health Study

The Black Women’s Health Study (BWHS) is a prospective cohort study established in 1995, in which 59,000 Black women aged 21–69 years from across the USA first completed health questionnaires. The baseline BWHS questionnaire included data from participants on demographic and lifestyle factors, reproductive history, dietary intake and medical conditions. The BWHS cohort has been followed every two years through mailed and online questionnaires. Follow-up has been successful for 85% of potential persons/year through the last completed follow-up in 2017/2018. Approximately 5% of the BWHS cohort is estimated to be foreign-born based on a nativity question asked in 1997. Of this number, three quarters (73%) were from the Caribbean and Central and South America [81]. In the current sample of 1000 women, 3% (n = 31) provided no information on the country of birth and 3% (n = 31) indicated that they were born outside of the USA, with over half reporting having been born in the Caribbean (59%), UK/western Europe: 24% and Asia: 10% and Africa: 7%.

In 2016, a random subset of BWHS participants who had completed the most recent BWHS survey and had previously provided a blood sample was invited via email to complete the Study on Stress, Spirituality and Health (SSSH) Baseline Spirituality Survey (SS-1). Within 2 weeks, more than 2600 BWHS participants had completed the SS-1 and data collection was stopped. Random sampling was then used to select 1000 participants whose stored blood samples were used to conduct the DNAm assays. The current analysis utilizes the measures of R/S and DNAm from the SSSH, as well other measures of psychosocial stress obtained in previous BWHS questionnaires.

Blood specimen collection

From 2012 to 2016, all BWHS participants were invited to provide a blood sample by going to a Quest Diagnostic Laboratory Patient Service Center at a location convenient to the participant [82]. During this period, each woman was sent a packet with an introductory letter, brochure, consent form, instructions for how to locate a blood collection site and a preprinted laboratory requisition form. Blood samples were collected and processed by Quest Diagnostics (NJ, USA; www.QuestDiagnostics.com), which is an accredited national clinical laboratory (CLIA). Samples, with signed informed consent, were provided by 13,300 BWHS participants. Women who did not provide blood samples were similar to participants who provided blood samples with regard to age, BMI, education, income, alcohol consumption, vigorous exercise, menopausal hormone use and prevalence of diabetes, hypertension and high cholesterol. Stored blood samples from this 2012–2016 period were used to perform the DNAm assays for analysis.

After the collection of blood samples, EDTA tubes were spun at approximately 1500–2500 revolutions per minute for 15 to 20 min. The buffy coat (white blood cell) layer was removed by pipet and aliquoted into cryotubes, frozen at -20°C and shipped on dry ice to the Boston University Core Genetics facility for long-term storage at -80°C. The QiAmp (Qiagen, CA, USA) 96-spin blood protocol was used to extract genomic DNA from peripheral blood leukocytes. Samples used for the methylation assays had been stored 1–5 years prior to DNA extraction; methylation assay was performed 1 year later.

Selection of genes & CpG sites

Based on the previous literature, the authors focused on four genes known to play prominent roles in HPA axis activity, and for some of which altered DNAm has been associated with stress and the development of various chronic conditions [53]. The HPA axis is the major neuroendocrine system that controls stress reactivity and regulates stress hormone (e.g., glucocorticoid) levels in the body. NR3C1, the GR gene, is one of the most studied genes in epigenetic analyses investigating the impact of psychosocial stress on DNAm [44,45,83]. The GR binds to glucocorticoids such as cortisol and regulates the production of many immune, inflammatory and metabolic proteins inside the cell nucleus. The authors’ lab was unable to successfully run the assay developed for a highly cited paper by Oberlander et al. [83] that measures DNAm in the NR3C1 1F promoter. Therefore, we instead included four CpGs located in a CpG island shore downstream of the CpG island located in the NR3C1 proximal promoter. CpGs were identified via bioinformatic analysis using the Genomatix software (Genomatix Software, Inc., MI, USA). This region was selected in light of the growing evidence that CpG island shores are enriched with functional methylation sites that control gene expression [84–86], and the authors’ recent research among BWHS participants showing that childhood abuse victimization was associated with methylation of one of these CpG (CpG7) island shores in a dose-response manner [65].

Three CpG sites located in the HSD11B2 gene, which cover three transcription factor binding sites, were selected. One CpG overlaps with one region where DNAm has been previously shown to control HSD11B2 in human cell lines and in rats [87]. With respect to the HSD11B1 gene, two CpGs were included that were previously analyzed in the only extant human study of HSD11B1 methylation using pyrosequencing, which demonstrated an association between HSD11B1 methylation and being born large for gestational age [88]. The HSD11B2 enzyme regulates cortisol levels by degrading cortisol to cortisone, whereas HSD11B1 catalyzes the conversion of the inactive cortisone to active cortisol [89]. Lastly, seven loci in the FKBP5 gene, all located in intron 7 close to glucocorticoid response elements, were analyzed. This assay interrogates the same bins of CpGs used in previous studies assessing the associations of early childhood exposures and FKBP5 methylation [72,90]. FKBP5 is a protein cochaperone that is acutely induced by stress and can influence GR signaling and glucocorticoid resistance/sensitivity [91,92].

DNA methylation sequencing & creation of dependent variables

Buffy coat samples were separated from the whole blood samples at Quest Diagnostics and immediately stored at -80°C. Genomic DNA was isolated according to standard protocols and stored at -80°C for future use [93–95]. A total of 500 ng of genomic DNA was treated using the EZ DNA Methylation Kit (Zymo Research, CA, USA) according to the manufacturer’s protocol. The final elution volume was 46 μl. Conversion efficiency was confirmed by verifying the methylation status of non-CpG cytosines. Bisulfite-treated DNA was aliquoted and stored at -80°C until ready for use. For each pyrosequencing assay, a PCR reaction was performed using HotStar Taq Polymerase (Qiagen) in combination with standard reagents and 0.6 ul of a 10 μM forward and reverse primer, one of which was biotinylated. Cycling conditions and target sequences for each gene are listed in Supplementary Table 1. The PCR product was then bound to Streptavidin Sepharose HP beads (GE Healthcare Bio-Sciences AB, Uppsala, Sweden). The Sepharose beads containing the immobilized PCR product were purified, washed and denatured using 0.2 M NaOH solution, and then rewashed using the pyrosequencing Vacuum Prep Tool (Qiagen) according to the manufacturer’s protocol. A total of 10 μl of the purified single-stranded PCR product and 0.5 μM of the pyrosequencing primer were then pyrosequenced on the PSQ96 HS System (Qiagen), according to the manufacturer’s instructions.

The methylation status of each CpG site was analyzed individually using QCpG software (Qiagen). For each sample, the methylation level at each CpG site was calculated as the percentage of the methylated cells over the sum of methylated and unmethylated cells. Samples were sequenced in duplicate, with a correlation of more than 99.8%. The average methylation level across all CpGs for each gene was calculated and used in further analyses. A pyrogram peak pattern from every sample was visually inspected to confirm the quality of the reaction, and four controls wells (no template, high, medium and low methylation controls) were included on each plate to ensure data quality. The four target CpG locations for NR3C1 were Chr5: 143401272, 143401140, 143401129 and 143401102; the two target CpG locations for HSD11B1 were Chr1: 209704596 and 209704625; the three target CpG locations for HSD11B2 were Chr16: 67430021, 67430027 and 67430031; and the seven target CpG locations for FKBP5 were Chr6: 35590945, 35590934, 35590737, 35590720, 35590712, 35590662 and 35590610 (Genome Reference Consortium Human Reference 38 Genome [GRCh38]). PCR and pyrosequencing were carried out by EpigenDx (MA, USA). The dependent variables included in this study represent the methylation level at each CpG site, as well as one variable for each gene reflecting the average methylation level across all CpGs within each individual gene. See Supplementary Figure 1 for visual representations of the genomic locations of the investigated CpGs.

Genotyping

Genotyping was provided by the Dana-Farber/Harvard Cancer Center Genotyping and Genetics for Population Sciences core facility, a unit of the Partners HealthCare Center for Personalized Genetic Medicine. Genotyping of aN SNP, rs1360780, was completed on DNA extracted from whole blood normalized to 10 ng/uL concentrations and was completed using the TaqMan^® OpenArray^® platform (Life Technologies, CA, USA). Samples were plated randomly with respect to time of recruitment and 5% of samples were duplicated for quality control purposes.

Independent variables

Childhood psychosocial stressors

Childhood loss of a parent or guardian

On the 2017 questionnaire, participants were asked: ‘Before the age of 18, did you lose (either from death or prolonged separation) your: mother, father, guardian, not applicable’ [96]. The variable was categorized as any loss (loss of mother, father or guardian) versus none.

Family member in prison

Also in 2017, a question adapted from White et al. [97] asked, ‘During the first 18 years of your life, did anyone in your household serve time in prison?’ Yes/No.

Childhood financial hardship

On the 2011 follow-up questionnaire, financial hardship as a child was assessed using the following two questions about life experience up to age 11: ‘Was there at least one time when your household did not have enough money for food or housing?’ and ‘Was there at least one time when your household received public assistance or welfare?’ Women who responded positively to either question were considered to have experienced financial hardship as a child.

Childhood abuse victimization

Procedures for collection of childhood abuse data in the BWHS have been previously published [98]. On the 2005 BWHS follow-up questionnaire, participants were asked questions about physical and sexual abuse (‘abuse victimization’) across the lifespan, including exposure as a child (up to age 11), teenager (aged 12–18) and adult (aged 19 and older). The nine-item abuse instrument was adapted from the Conflict Tactics Scale [99] and the Pregnancy Abuse Assessment Screen [100]. Response categories were ‘never’, ‘1–3 times’, or ‘≥4 times’. Any childhood physical abuse was defined as any report of a perpetrator having ‘pushed, grabbed, or shoved me’, ‘threw something at me that could hurt me’, ‘kicked, bit, or punched me’, ‘hit me with something including hand or fist’ or ‘physically attacked me in some other way’ at a frequency of ≥4 times; or either ‘choked or burned me’ or ‘seriously harmed someone I loved’ at any frequency. Any childhood sexual abuse was defined as any report of a perpetrator having ‘exposed genitals against my will’ ≥4 times; or ‘been sexual with me against my will’ at any frequency.

As in previous research [101], the authors created physical and sexual abuse severity scores by assigning 1 point for each report of a physical abuse item occurring ≥4 times, with the exception of ‘choked or burned’ or ‘seriously harmed someone I loved’, where 1 point was assigned for reports that these occurred 1–3 times and 2 points for reports that these occurred ≥4 times because these events were considered more severe. The resulting abuse severity scores, which ranged from 0 to 9, were categorized as mild (score = 1), moderate (score = 2) and severe (score ≥3). Severe abuse was defined as 3 or more types of physical abuse and sexual abuse occurring ≥4 times. Further details on these abuse definitions have been described in previous publications from the BWHS [101–105] and other studies [106,107].

Lack of childhood emotional support

The 2011 BWHS questionnaire queried women about receipt of nurturing and emotional support in childhood using the following two questions: ‘When you were growing up, did people in your family show confidence in you and encourage you to achieve?’ and ‘When you were growing up, did you feel that there was someone to take care of you and protect you?’ Response categories were ‘never’, ‘almost never’, ‘sometimes’, ‘fairly often’ and ‘very often’ [108]. Women who answered ‘never’ or ‘almost never’ to both questions were considered to have experienced a lack of emotional support in childhood [65,109].

Adulthood psychosocial stressors

Sleep duration

The 2009 BWHS questionnaire asked when the participant usually fell asleep and when she awoke over the past two years, from which sleep duration was calculated. The total hours of sleep were categorized as <7, 7–8 and ≥9. Only seven of the 1000 participants in these analyses reported ≥9 h of sleep and therefore estimates for that category are not presented.

Depression diagnosis

A self-reported physician diagnosis of depression defined as depression treated with antidepressant use was ascertained in 2013.

Depressive symptoms

Experience of depressive symptoms was measured using the Center for Epidemiologic Studies Depression Scale (CES-D) [110] in 2005. For participants with missing CES-D scores in 2005, their score from the 1999 questionnaire was used, if available. The CES-D score is routinely categorized as ‘absence of depressive symptoms’ (score <16), moderate depressive symptoms (score 16–23) and severe depressive symptoms (score ≥24) [111].

Interpersonal racism

The 2009 follow-up questionnaire included five items on experiences of racism adapted from an instrument developed by Williams et al. [112]. These items assessed the frequency of interpersonal racism in daily life according to the following experiences: ‘you receive poorer service than other people in restaurants or stores’, ‘people act as if they think you are not intelligent’, ‘people act as if they are afraid of you’, ‘people act as if they think you are dishonest’ and ‘people act as if they are better than you’. Response options were ‘never’, ‘a few times a year’, ‘once a month’, ‘once a week’ and ‘almost every day’. An interpersonal racism score was created by averaging subjects’ responses to the five questions; the score was then divided into four quantiles ranging from low to high.

Lifetime institutional racism

In 2009, three questions ascertained lifetime exposure to institutional racism by asking each participant whether she had ever been ‘treated unfairly due to your race’ in three contexts: on the job, in housing or by the police. Response categories were ‘yes’ and ‘no’. [81] A lifetime institutional racism score summed the positive responses to obtain an overall score of 0, 1, 2 or 3, with a score of 3 indicating exposure to all three types of institutional racism.

Religion & spirituality

The SS-1 includes a broad range of measures capturing religious and spiritual practices, beliefs and experiences, and was developed with input from SSSH investigators, cohort principal investigators and national R/S experts. The questionnaire was reviewed by the BWHS Community Advisory Board, with some modifications to language subsequently made to reflect cultural preferences while maintaining item validity. A subset of R/S items across conceptual themes was selected for this analysis. A recent psychometric analysis of R/S items and scales from the SS-1 found this selection of items from established scales and de novo measures to have reliable psychometric properties [113].

Religious service attendance

One item asked how frequently respondents attended religious services. Response categories were ‘never’, ‘rarely’, ‘about once a month’, ‘2–3 times per month’, ‘once a week’ and ‘several times per week’. For stratification, the variable was dichotomized to less than once a week versus once a week or more.

Quality of congregational relationships

Three items assessed the quality of individuals’ relationships with their congregational community: ‘how often do you show someone in your congregation or religious community that you love or care for them?’ [114], ‘how often are people in your congregation or religious community critical of you or your lifestyle?’ [115] and ‘how often do you feel ignored or neglected by people in your congregation or religious community?’ (de novo). Response categories for all three were: ‘never’, ‘once in a while’, ‘fairly often’ and ‘very often.’ In stratified models, the ‘love and care’ variable was dichotomized to ’very often’ versus all other response options, while the ‘critical of you’ and ’ignored or neglected’ variables were categorized as ‘never’ versus all other response options.

Positive religious coping

A single item assessed positive religious coping: ‘to what extent is your religion or spirituality involved in understanding or dealing with stressful situations?’ [115]. Response categories were ‘not at all’, ’not very much’, ’somewhat’ and ‘very much so’. In stratified analyses, a binary variable comparing responses of ‘very much so’ with all other responses was used.

A positive religious coping scale was also created, based on eight subitems reflecting four domains of the longer R-COPE [116]: ‘I saw my situation as part of God’s plan’, ‘I tried to see how God might be trying to strengthen me in these situations’, ‘I tried to make sense of the situation with God’, ‘I worked together with God to relieve my worries’, ‘I did what I could and put the rest in God’s hands’, ‘I took control over what I could, and gave the rest up to God’, ‘I sought God’s love and care’ and ’I trusted that God would be by my side’. Two de novo items developed as a result of focus groups the authors’ team conducted with Blacks and Hispanics/Latinos exploring the relationship between their spirituality and their health: ’I felt hopeful that God would help me get through one day at a time’ and ‘I looked to my faith in God for hope about the future’. Response categories were ‘not at all’, ‘somewhat’, ‘quite a bit’ and ‘a great deal’. An average positive coping score was calculated by summing responses for each item and dividing by eight. Women were separated into two groups according to a score in the top quartile versus those in the bottom three quartiles.

Spiritual struggles

The potential negative impact of one’s religion or spirituality on her health was assessed using six of the 35 ‘negative religious coping’ items from Pargament’s R-COPE [116], selected in collaboration with the scale’s creator, which capture negative religious beliefs/spiritual struggles one has when trying to cope with stressful life situations: ‘I wondered what I did for God to punish me’; ‘I wondered if God allowed this event to happen to me because of my wrongdoings’; ‘I believed the devil or evil spirits were responsible for my situation’; ‘I felt as though the devil or an evil spirit was trying to turn me away from God’; ‘I wondered whether God had abandoned me’; and ‘I questioned God’s love or care for me’. Along with two items from the Religious and Spiritual Struggles Scale [117]: ‘I felt confused about my religious or spiritual beliefs’ and ‘I felt troubled by doubts or questions about my religion or spirituality’. Response categories were ‘not at all’, ‘somewhat’, ‘quite a bit’ and ‘a great deal.’ As with the measure of positive religious coping, women’s responses were summed and divided by eight to generate an average negative religious coping/spiritual struggles score, hereafter referred to as ‘spiritual struggles’. Women were separated into two groups according to those with a score in the top quartile versus those in the bottom three quartiles and this variable was used in interaction analyses.

Nontheistic daily spiritual experiences

The Daily Spiritual Experiences Scale (DSES) [118] captures both God-focused and nontheistic everyday experiences of transcendence. The four nontheistic items of the DSES were used to create a nontheistic DSES subscale (α = 0.75): ‘I experience a connection to all of life’, ‘I feel deep inner peace or harmony’, ‘I am touched by the beauty of creation’ and ‘I feel a selfless caring for others.’ Response categories were ‘never’, ‘once in a while’, ‘some days’, ‘every day’ and ‘many times a day’. A score was created by summing responses to each item and dividing by four. Women were separated into two groups according to those with a score in the top quartile versus those in the bottom three quartiles.

Statistical analyses

All analyses were carried out using R software (version 4.0.2). To account for the random effect (intercept) of plate, linear mixed models were used to estimate beta values and 95% CIs for the association between each psychosocial stressor and percent of cells methylated in each of the select CpG sites across the four HPA axis genes under study (i.e., HSD11B1/2, NR3C1 and FKBP5). Four CpGs were assessed in NR3C1, seven CpGs were assessed in FKBP5, two CpGs were assessed in HSD11B1 and three CpGs were assessed in HSD11B2. Associations between psychosocial stressors and the mean methylation levels were first assessed for all CpG loci within each gene. Sensitivity analyses assessing the relationship between stressors and individual CpG sites was then performed. Minimally adjusted models included terms for age at the time of blood draw, white blood cell count and proportion of lymphocytes, monocytes, basophils, eosinophils and neutrophils. Multivariable models additionally controlled for alcohol consumption (never, past, current 1–6 drinks/week, current ≥7 drinks/week), body mass index (continuous kg/m²), years of education (≤12, 13–15, ≥16), marital status (married or living as married, single, widowed, divorced or separated), physical activity (none, <1, 1, 2, 3–4, 5–6, 7–9, ≥10 h per week) and smoking status (never, former or current). With respect to the FKBP5 gene, previous research has shown that a single SNP, rs1360780, modulates the impact of stress-related exposures and early life trauma on FKBP5 DNAm [72]. Analyses of FKBP5 were therefore stratified into two FKBP5 genotypes: CT/TT (‘risk’ genotype) and CC (‘protective’ genotype).

To explore whether measures of R/S practices and beliefs modified associations between the various psychosocial exposures and DNAm levels of individual CpGs, interaction analyses were performed by adding an interaction term (stressor * binary R/S variable) for each of the main effect associations that reached a minimal association threshold of p < 0.05, prior to false discovery rate (FDR) correction. A final set of comprehensive interaction analyses was then completed to test whether certain stress exposures were only significantly associated with mean DNAm in the presence of an interaction with the R/S variables. Specifically, an interaction for every combination of stress and R/S variables in models with mean DNAm of the four genes as the dependent variable was performed, resulting in 400 interaction models tested.

Results

Table 1 displays the distribution of psychosocial stressors in childhood and adulthood used in the analysis and participant characteristics at the time of blood draw among the 1000 Black Women’s Health Study participants who completed the SSSH SS-1 survey. Table 2 displays the R/S measures used in the analysis among the 1000 Black Women’s Health Study participants who completed the SS-1. Tables 3–5 display the average DNAm distributions for each gene and FKBP5 genotype.

Table 1. Characteristics of Black Women’s Health Study participants with methylation data (n = 1000).

Variable	n	%
Childhood stressors
Death of a parent or caretaker
No	786	79.3%
Yes, any	205	20.7%
Parent or sibling imprisoned in childhood
No	924	93.1%
Yes	69	6.9%
Financial hardship in childhood
No	627	66.4%
Yes	317	33.6%
Lack of emotional support in childhood
No	649	65.7%
Yes	339	34.3%
Severity of sexual abuse in childhood
None	366	57.6%
Moderate	146	23.0%
Severe	123	19.4%
Severity of physical abuse in childhood
None	366	40.4%
Mild	179	19.8%
Moderate	96	10.6%
Severe	265	29.3%
Adulthood stressors
CESD score
<16	774	78.1%
16–23	118	11.9%
24+	99	10.0%
Interpersonal racism score
Quartile 1	259	25.9%
Quartile 2	180	18.0%
Quartile 3	331	33.1%
Quartile 4	230	23.0%
Lifetime institutional racism
No to all	300	30.0%
Yes to 1	321	32.1%
Yes to 2	245	24.5%
Yes to 3	134	13.4%
Sleep per day
<7 h	536	53.7%
7 h	293	29.3%
8 h	148	14.8%
9 h	16	1.6%
10+ h	6	0.6%
Missing	1	0.0%
Marital status
Widowed, divorced, separated	291	30.3%
Single	233	24.2%
Married or living as married	437	45.5%
Other characteristics (at time of blood draw)
Alcohol consumption
Never	696	69.6%
Past	59	5.9%
Current 1–6 drinks/week	206	20.6%
Current 7+ drinks/week	39	3.9%
Cigarette smoking
Never	727	72.7%
Current	75	7.5%
Past	198	19.8%
Vigorous physical activity
None or very little	656	65.7%
1–2 h per week	152	15.2%
3–4 h per week	113	11.3%
5+ h per week	78	7.8%
Length of education (years)
≤12	71	7.1%
12–15	241	24.1%
>16	687	68.8%
Age (mean, SD)	56.3 (7.5)
BMI (mean, SD)	30.7 (6.9)

CESD: Center for Epidemiologic Studies Depression Scale; SD: Standard deviation.

Table 2. Characteristics of religion and spirituality variables among Black Women’s Health Study participants with methylation data (n = 1000).

Variable	n	%
Religious service attendance
Missing	0	–
Less than once a week	393	39.3
Once a week or more	607	60.7
Showing love to religious community
Missing	330	–
Never, once in a while or fairly often	304	45.4
Very often	366	54.6
Criticism from religious community
Missing	330	–
Never	458	68.4
Once in a while, fairly often or very often	212	31.6
Ignored or neglected by members of one’s religious community
Missing	330	–
Never	431	64.3%
Once in a while, fairly often or very often	239	35.7
Importance of one’s religion/spirituality in coping with stressful life situations
Missing	17	–
Not at all	21	2.1
Not very much	37	3.8
Somewhat	179	18.2
Very much so	746	75.9
Nontheistic daily spiritual experiences
Missing	1	–
Bottom three quartiles	767	76.8
Top quartile	232	23.2
Positive religious coping score
Missing	18	–
Bottom three quartiles	819	83.4
Top quartile	163	16.6
Negative religious coping score
Missing	18	–
Bottom three quartiles	761	77.5
Top quartile	221	22.5

Table 3. Average percent methylation levels for all CpGs and mean DNA methylation in NR3C1, HSD11B1 and HSD11B2 among 1000 Black Women’s Health Study participants.

DNA methylation	Mean	SD	Minimum	Q1	Median	Q3	Maximum	Valid (n)
HSD11B1 CpG 5	85.61	4.17	72.40	82.50	85.50	88.55	98.50	1000
HSD11B1 CpG 6	89.81	2.10	82.70	88.30	90.40	91.30	100	1000
HSD11B2 CpG 83	4.30	2.15	0	3.10	4.10	5.40	13.40	997
HSD11B2 CpG 84	11.04	2.44	0	9.40	10.80	12.50	19.60	997
HSD11B2 CpG 85	2.66	1.79	0	1.70	2.90	3.70	9.60	997
HSD11B1 mean DNAm	87.71	2.00	80.40	86.50	87.80	89.05	96.50	1000
HSD11B2 mean DNAm	6.00	1.76	0	4.80	5.87	7.03	12.93	997
NR3C1 mean DNAm	75.22	3.06	59.43	73.80	75.58	77.10	82.45	992
NR3C1 CpG 4	89.91	4.61	72.80	87.10	90.05	92.90	100	994
NR3C1 CpG 5	83.41	5.36	42.80	80.50	83.90	87.10	95.90	998
NR3C1 CpG 6	75.60	5.98	37.40	73.80	76.50	79.10	86.60	998
NR3C1 CpG 7	51.95	4.98	33.60	48.70	51.90	55.00	70.70	997

DNAm: DNA methylation; Q: Quartile; SD: Standard deviation.

Table 4. Average percent methylation levels for all CpGs and mean DNA methylation in FKBP5 (CC genotype) among 1000 Black Women’s Health Study participants.

DNA methylation locus	Mean	SD	Minimum	Q1	Median	Q3	Maximum	Valid (n)
FKBP5 CpG 539	64.42	5.03	46.80	61.10	64.60	67.90	80.80	329
FKBP5 CpG 540	61.67	6.46	16.80	57.80	61.60	65.10	80.70	329
FKBP5 CpG 541	76.31	3.73	52.10	74.10	76.60	78.90	85.60	329
FKBP5 CpG 542	96.03	2.51	68.20	94.80	95.80	96.80	100	329
FKBP5 CpG 543	92.66	2.18	86.40	91.70	92.60	93.30	100	329
FKBP5 CpG 544	79.84	6.82	57.10	75.50	79.50	84.30	100	329
FKBP5 CpG 545	81.45	4.32	37.20	79.25	81.90	83.90	90.20	328
FKBP5 mean DNAm	78.92	2.82	65.50	77.03	78.75	80.89	86.81	328

DNAm: DNA methylation; Q: Quartile; SD: Standard deviation.

Table 5. Average percent methylation levels for all CpGs and mean DNA methylation in FKBP5 (CT/TT genotype) among 1000 Black Women’s Health Study participants.

DNA methylation locus	Mean	SD	Minimum	Q1	Median	Q3	Maximum	Valid (n)
FKBP5 CpG 539	61.42	6.16	41.30	57.30	62.00	65.70	77.70	658
FKBP5 CpG 540	60.49	5.81	35.30	56.80	60.25	64.10	78.90	658
FKBP5 CpG 541	75.78	4.90	54.40	73.60	76.80	78.80	100	659
FKBP5 CpG 542	95.84	2.76	53.60	94.80	95.60	96.50	100	659
FKBP5 CpG 543	92.41	3.37	46.20	91.50	92.50	93.10	100	659
FKBP5 CpG 544	79.01	7.63	50.20	73.40	78.50	84.50	100	659
FKBP5 CpG 545	81.07	6.05	34.10	79.10	81.90	84.20	90.60	657
FKBP5 mean DNAm	78.00	3.29	65.51	76.06	78.24	79.99	87.37	656

DNAm: DNA methylation; Q: Quartile; SD: Standard deviation.

Table 6 shows differences in mean methylation levels for each gene according to the levels of the psychosocial stress variables; CpG-specific results are shown in Supplementary Tables 2–4. Figure 1 shows a heatmap of mean and CpG-specific associations found with stress exposures, although no associations shown in Figure 1 remained significant after FDR correction.

Table 6. Multivariable-adjusted differences in mean NR3C1, HSD11B1/2 and FKBP5 percent methylation levels according to experience of psychosocial stress among 1000 Black Women’s Health Study participants.

Exposure	HSD11B1 mean DNAm	HSD11B2 mean DNAm	NR3C1 mean DNAm	FKBP5 mean DNAm (CC genotype)	FKBP5 mean DNAm (CC/TT genotype)
	Beta (95% CI)	Beta (95% CI)	Beta (95% CI)	Beta (95% CI)	Beta (95% CI)
Childhood stressors
Death of a parent or caretaker
No	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)
Yes, any	-0.085 (-0.360, 0.190)	0.111 (-0.130, 0.353)	-0.398 (-0.831, 0.035)	-0.246 (-0.830, 0.337)	0.041 (-0.467, 0.548)
p-value	0.543	0.367	0.072	0.409	0.875
FDR p-value	0.900	0.881	0.878	0.818	0.933
Incarceration of a family member
No	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)
Yes	0.263 (-0.173, 0.699)	0.235 (-0.148, 0.617)	-0.371 (-1.059, 0.317)	-0.056 (-1.031, 0.918)	0.351 (-0.439, 1.142)
p-value	0.238	0.230	0.291	0.91	0.384
FDR p-value	0.881	0.881	0.881	0.967	0.768
Financial hardship
No	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)
Yes	-0.002 (-0.235, 0.232)	0.087 (-0.123, 0.297)	0.057 (-0.311, 0.425)	-0.114 (-0.625, 0.397)	0.271 (-0.171, 0.713)
p-value	0.988	0.418	0.762	0.661	0.230
FDR p-value	0.988	0.881	0.904	0.967	0.575
Lack of emotional support
No	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)
Yes	-0.195 (-0.436, 0.045)	0.060 (-0.151, 0.270)	-0.178 (-0.556, 0.200)	-0.634 (-1.147, -0.122)	-0.024 (-0.457, 0.409)
p-value	0.111	0.580	0.357	0.016	0.913
FDR p-value	0.878	0.900	0.881	0.160	0.933
Sexual abuse severity
None	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)
Moderate	-0.148 (-0.500, 0.203)	-0.091 (-0.397, 0.215)	-0.308 (-0.843, 0.227)	-0.814 (-1.587, -0.042)	0.190 (-0.432, 0.813)
Severe	0.132 (-0.239, 0.504)	-0.028 (-0.352, 0.295)	0.045 (-0.522, 0.612)	-0.161 (-1.008, 0.687)	0.622 (-0.034, 1.279)
p-value	0.444	0.844	0.470	0.121	0.179
FDR p-value	0.881	0.904	0.881	0.523	0.575
Physical abuse severity
None	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)
Mild	0.075 (-0.243, 0.394)	-0.100 (-0.378, 0.178)	0.266 (-0.229, 0.76)	0.037 (-0.583, 0.657)	0.338 (-0.243, 0.918)
Moderate	0.002 (-0.397, 0.4)	-0.121 (-0.468, 0.227)	-0.040 (-0.659, 0.579)	-0.791 (-1.561, -0.02)	0.780 (0.049, 1.511)
Severe	-0.084 (-0.365, 0.197)	0.072 (-0.174, 0.318)	0.076 (-0.36, 0.511)	-0.204 (-0.76, 0.353)	0.185 (-0.336, 0.706)
p-value	0.824	0.600	0.732	0.209	0.194
FDR p-value	0.904	0.900	0.904	0.523	0.575
Adulthood stressors
CESD score
<16	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)
16–23	-0.097 (-0.440, 0.247)	0.195 (-0.108, 0.497)	-0.391 (-0.939, 0.156)	-0.079 (-0.831, 0.672)	0.336 (-0.289, 0.962)
≥24	-0.201 (-0.577, 0.175)	-0.020 (-0.353, 0.313)	-0.044 (-0.641, 0.553)	-0.070 (-0.800, 0.659)	-0.777 (-1.521, -0.033)
p-value	0.529	0.434	0.375	0.967	0.052
FDR p-value	0.900	0.881	0.881	0.967	0.520
Interpersonal racism score
Quartile 1	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)
Quartile 2	0.005 (-0.338, 0.348)	0.198 (-0.102, 0.498)	0.103 (-0.440, 0.645)	-0.078 (-0.849, 0.692)	-0.167 (-0.787, 0.453)
Quartile 3	0.045 (-0.247, 0.337)	0.139 (-0.117, 0.395)	-0.018 (-0.481, 0.445)	0.040 (-0.580, 0.659)	-0.136 (-0.673, 0.401)
Quartile 4	-0.093 (-0.411, 0.226)	0.425 (0.146, 0.704)	-0.114 (-0.619, 0.39)	-0.204 (-0.879, 0.471)	0.261 (-0.324, 0.846)
p-value	0.843	0.026	0.897	0.888	0.472
FDR p-value	0.904	0.780	0.928	0.967	0.787
Lifetime institutional racism
None	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)
Yes to 1	-0.177 (-0.459, 0.105)	0.19 (-0.057, 0.438)	-0.052 (-0.499, 0.395)	-0.147 (-0.747, 0.453)	0.021 (-0.492, 0.534)
Yes to 2	-0.09 (-0.395, 0.214)	0.327 (0.06, 0.594)	-0.036 (-0.518, 0.447)	-0.006 (-0.625, 0.613)	0.091 (-0.478, 0.66)
Yes to 3	-0.019 (-0.389, 0.35)	0.174 (-0.151, 0.498)	-0.322 (-0.909, 0.264)	-0.857 (-1.684, -0.03)	-0.246 (-0.913, 0.42)
p-value	0.633	0.117	0.736	0.192	0.803
FDR p-value	0.904	0.878	0.904	0.523	0.933
Sleep duration
7 or 8 h	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)	Referent (0.00)
<=6 h	-0.021 (-0.245, 0.203)	-0.138 (-0.334, 0.059)	-0.206 (-0.560, 0.148)	0.101 (-0.384, 0.586)	-0.076 (-0.482, 0.33)
<=9 h	0.505 (-0.878, 1.888)	0.603 (-0.609, 1.816)	0.643 (-1.536, 2.822)	0.249 (-2.665, 3.164)	-0.107 (-2.618, 2.404)
p-value	0.752	0.217	0.419	0.913	0.933
FDR p-value	0.904	0.881	0.881	0.967	0.933

Models are adjusted for age, white blood cell count, smoking, body mass index, education, alcohol, marital status and physical activity with a random effect for plate.

DNAm: DNA methylation; CESD: Center for Epidemiologic Studies Depression Scale; FDR: False discovery rate.

Figure 1. Heatmap of CpG-specific and mean DNA methylation changes of HSD11B1/2, NR3C1 and FKBP5 that are nominally associated with childhood and adulthood stress among Black Women’s Health Study participants with methylation data (n = 1000; note: no associations survive false discovery rat).

Blank cells represented here remained significant after false discovery rat correction. Cells shaded red represent an association with a beta coefficient >0 and a nominal p-value <0.05. Cells shaded blue represent an association with a beta coefficient <0 and a nominal p-value <0.05. Models are adjusted for age, white blood cell count, smoking, BMI, education, alcohol, marital status and physical activity with a random effect for plate.

CESD: Center for Epidemiologic Studies Depression Scale.

NR3C1 methylation

In multivariable models assessing associations with the mean DNAm levels across the four CpG sites in NR3C1 (Table 6), no significant associations were observed with any of the psychosocial stressors examined after FDR. Nominally significant results were observed in individual CpG sites for loss of a parent or caretaker (CpG4 β: -0.719 [-1.239, -0.199], p-value = 0.007 and CpG5 β: -0.626 [-1.237, -0.015], p-value = 0.045). See Supplementary Table 2 for CpG-specific results.

HSD11B1/2 methylation

In multivariable models for mean DNAm level across the two examined CpG sites in HSD11B1 and the three examined HSD11B2 CpG sites, nominal associations were observed with one stressor, interpersonal racism, in relationship to mean HSD11B2 DNAm (Table 6). Mean percent DNAm levels were higher in individuals in the highest quartile of the interpersonal racism score (β: 0.425 [0.146, 0.704], p-value = 0.026). This association with interpersonal racism was driven by the higher DNAm levels in two of the three HSD11B2 CpG sites examined: CpG83 and CpG84 (Supplementary Table 2). None of these associations survived FDR correction.

FKBP5 methylation

FKBP5 CC genotype

Among women with the protective CC genotype for FKBP5, nominally significant associations were observed between lack of emotional support in childhood and decreased mean FKBP5 DNAm (β: -0.634 [-1.147, -0.122], p-value = 0.016). Lack of emotional support in childhood and CESD depression score were also associated with decreased DNAm at CpG540 and CpGs541–542, respectively (Supplementary Table 3). None of these associations survived FDR correction.

FKBP5 CT/TT genotype

Among participants with the CT or TT (risk) genotype, none of the psychosocial stressors examined were associated with mean DNAm across any FKBP5 CpG sites (Table 6). Nominally significant associations were observed in relation to increased DNAm of individual CpG sites for severe sexual abuse in childhood (CpG541), physical abuse in childhood (CpG540) and financial hardship in childhood (CpG545). See Supplementary Table 6 for CpG-specific results. None of these associations survived FDR correction.

The influence of R/S on psychosocial stress & DNA methylation

No stress × DNAm associations with a nominally significant p-value <0.05 showed a significant interaction with any R/S variables for mean DNAm. However, when testing systematically for significant interactions between any stress variable and any R/S variable in models with mean DNAm of the four genes as the dependent variable, nine significant interactions were found (Table 7). None of the significant p-values survived FDR correction after correcting for all 400 tests conducted. Stress variables from the nine nominally significant interaction models were then retested in a model stratified by the corresponding R/S variable from the interaction models. Among these, three stratification models were nominally significant (Table 8). Interestingly, financial hardship in childhood was only associated with mean NR3C1 methylation and only among those whose R/S was not involved in coping with stress, whereas among those with high levels of R/S coping there was no association between childhood financial hardship and mean NR3C1 methylation. Financial hardship in childhood was not associated with the mean of CpG-specific DNAm of NR3C1 in the main effects analysis.

Table 7. Significant interactions between stress and measures of religiosity/spirituality in the presence of mean DNA methylation of all HPA axis genes.

Gene	Stressor	R/S variable	p-value for interaction term^†
NR3C1 mean DNAm	Financial hardship	Importance of one’s religion or spirituality in coping with stressful life situations	0.003866671
HSD11B2 mean DNAm	Lifetime institutional racism	Showing members of your religious community that you love or care for them	0.003900336
HSD11B2 mean DNAm	CESD score	Nontheistic daily spiritual experiences	0.004238321
HSD11B2 mean DNAm	Lifetime institutional racism	Nontheistic daily spiritual experiences	0.008571150
HSD11B1 mean DNAm	Physical abuse severity	Nontheistic daily spiritual experiences	0.009625497
HSD11B1 mean DNAm	CESD score	Feeling criticized by members of one’s religious community	0.009891924
FKBP5 mean DNAm	(CC genotype)	Physical abuse severity Feeling ignored or neglected by members of one’s religious community	0.002672422
FKBP5 mean DNAm (CC genotype)	Physical abuse severity	Presence of religious/spiritual struggles	0.003870147
FKBP5 mean DNAm (CC genotype)	Sexual abuse severity	Importance of one’s religion or spirituality in coping with stressful life situations	0.009831949

Models are adjusted for age, white blood cell count, smoking, BMI, education, alcohol, marital status and physical activity with a random effect for plate.

† No p-values survived FDR correction.

DNAm: DNA methylation; CESD: Center for Epidemiologic Studies Depression Scale; FDR: False discovery rate; R/S: Religiosity or spirituality.

Table 8. Stratified results for stressors and mean DNA methylation stratified by religiosity/spirituality variables with a nominally significant interaction term.

Gene	Stressor	Stratified group	Estimate	p-value
NR3C1 mean DNAm	Financial hardship	Importance of one’s religion or spirituality in coping with stressful life situations	–	–
	Yes	Not at all, not very much, somewhat	1.218 (0.493, 1.943)	0.001
	Yes	Very much so	-0.152 (-0.583, 0.28)	0.491
HSD11B2 mean DNAm	Lifetime institutional racism	Showing members of your religious community that you love or care for them	–	–
	Yes to 1	Never, once in a while or fairly often	0.289 (-0.212, 0.79)	0.313
	Yes to 2	Never, once in a while or fairly often	-0.018 (-0.576, 0.54)	–
	Yes to 3	Never, once in a while or fairly often	0.481 (-0.137, 1.098)	–
	Yes to 1	Very often	-0.195 (-0.612, 0.222)	0.001
	Yes to 2	Very often	0.570 (0.134, 1.006)	–
	Yes to 3	Very often	-0.327 (-0.872, 0.218)	–
HSD11B2 mean DNAm	CESD score	Nontheistic daily spiritual experiences	–	–
	16–23	Bottom three quartiles	0.084 (-0.246, 0.414)	0.829
	≥24	Bottom three quartiles	-0.052 (-0.423, 0.32)	–
	16–23	Top quartile	1.051 (0.182, 1.921)	0.060
	≥24	Top quartile	0.226 (-0.568, 1.021)	–
HSD11B2 mean DNAm	Lifetime institutional racism	Nontheistic daily spiritual experiences	–	–
	Yes to 1	Bottom three quartiles	0.187 (-0.102, 0.476)	0.284
	Yes to 2	Bottom three quartiles	0.283 (-0.033, 0.598)	–
	Yes to 3	Bottom three quartiles	0.281 (-0.104, 0.666)	–
	Yes to 1	Top quartile	0.173 (-0.318, 0.665)	0.106
	Yes to 2	Top quartile	0.617 (0.101, 1.133)	–
	Yes to 3	Top quartile	0.040 (-0.560, 0.640)	–
HSD11B1 mean DNAm	Physical abuse severity	Nontheistic daily spiritual experiences	–	–
	Mild	Bottom three quartiles	0.374 (0.011, 0.738)	0.163
	Moderate	Bottom three quartiles	0.043 (-0.419, 0.505)	–
	Severe	Bottom three quartiles	-0.021 (-0.338, 0.296)	–
	Mild	Top quartile	-0.868 (-1.592, -0.145)	0.128
	Moderate	Top quartile	-0.098 (-0.949, 0.752)	–
	Severe	Top quartile	-0.287 (-0.969, 0.394)	–
HSD11B1 mean DNAm	CESD score	Religious community criticizing you	–	–
	16–23	Never	0.031 (-0.554, 0.615)	0.166
	≥24	Never	-0.562 (-1.149, 0.024)	–
	16–23	Once in a while, fairly often or very often	-0.303 (-0.997, 0.391)	0.423
	≥24	Once in a while, fairly often or very often	0.455 (-0.649, 1.558)	–
FKBP5 mean DNAm	Physical abuse severity	Religious community neglecting you	–	–
(CC genotype)	Mild	Never	0.718 (-0.391, 1.827)	0.376577
	Moderate	Never	-0.636 (-1.967, 0.695)	–
	Severe	Never	0.017 (-0.978, 1.013)	–
	Mild	Once in a while, fairly often or very often	Models don’t converge	–
	Moderate	Once in a while, fairly often or very often	–	–
	Severe	Once in a while, fairly often or very often	–	–
FKBP5 mean DNAm	Physical abuse severity	Negative religious coping	–	–
(CC genotype)	Mild	Bottom three quartiles	0.008 (-0.671, 0.688)	0.227255
	Moderate	Bottom three quartiles	-0.692 (-1.491, 0.106)	–
	Severe	Bottom three quartiles	0.233 (-0.409, 0.876)	–
	Mild	Top quartile	0.274 (-1.714, 2.262)	0.019208
	Moderate	Top quartile	-3.589 (-7.671, 0.494)	–
	Severe	Top quartile	-2.166 (-3.849, -0.482)	–
FKBP5 mean DNAm	Sexual abuse severity	Positive religious coping	–	–
(CC genotype)	Moderate	Bottom three quartiles	-0.658 (-1.514, 0.198)	0.223621
	Severe	Bottom three quartiles	0.223 (-0.754, 1.201)	–
	Moderate	Top quartile	Models don’t converge	–
	Severe	Top quartile	–	–

Models are adjusted for age, white blood cell count, smoking, BMI, education, alcohol, marital status and physical activity with a random effect for plate.

DNAm: DNA methylation; CESD: Center for Epidemiologic Studies Depression Scale; FDR: False discovery rate.

Discussion

This study is the first prospective epigenomic analysis investigating the relation of a large and diverse set of psychosocial stressors to DNAm of four HPA axis genes, using data from 1000 BWHS participants who completed the SSSH survey. To the best of our knowledge, this is the first data assessing the potential modifying effects of various measures of R/S on associations between psychosocial stressors in childhood and adulthood and DNAm. In total, we examined 16 unique CpGs across the FOUR HPA axis genes. We assessed the associations of ten different childhood and adulthood psychosocial stressors in relation to DNAm of these genes. We further assessed the extent to which eight different measures of R/S modified the impact of these stressors on DNAm. We believe that study designs such as this, which test diverse types of psychosocial factors simultaneously across many CpG sites in multiple related genes in a large sample size will be crucial in critically assessing findings from previously published studies analyzing a single stressor and/or gene. Importantly, the results show that despite previous findings between stress and DNAm of HPA axis genes, consistent and reliable DNAm signals could not be detected across multiple CpGs within the same regions in all of the studied genes, and no DNAm signals remained significant after FDR correction.

The results also did not replicate our own team’s previous finding of an association of childhood abuse victimization with increased DNAm of CpG7 in NR3C1 in a smaller sample (n = 295) of BWHS participants [65]. In the current analysis of 1000 BWHS women, no association was seen at this CpG site. It is unlikely that differences in covariates included in the multivariable models explain the inconsistencies because age-adjusted models from the two studies also provide differing results. Sampling variation is a likely explanation.

Childhood maltreatment and sexual abuse have been previously associated with decreased DNAm of the same FKBP5 CpGs analyzed in studies that included Black participants [51,72] and several that did not [90,119]. However, a larger and more recent cohort analysis of German participants published by Klinger-König et al. also found no main effect or interaction effect between childhood trauma and DNAm of these CpGs in FKBP5 [119]. The authors mentioned that although the childhood trauma questionnaire (CTQ) subscales analyzed did show some significant results, their findings were inconsistent between the cohorts they studied and were also inconsistent in the direction of the effects. Consequently, they concluded that, although they found some single CpG associations between FKBP5 DNAm levels and childhood maltreatment, as we did, FKBP5 DNAm was not a strong biomarker of maltreatment.

Overall, the results indicate that there is not a reliable signal between retrospective assessment of childhood adversity and DNAm of these HPA axis genes in Black women. Importantly, this contradicts many previous findings that have come largely from smaller studies analyzing single exposures and/or DNAm in single-candidate genes or at-risk study populations. When analyzing multiple childhood and adulthood stressful experiences together across a panel of multiple HPA axis genes, we found no signals of altered DNAm that were stable across more than two nearby CpGs. Furthermore, none of the significant associations at the level of a single individual CpG survived FDR correction.

A recently published systematic review that gives a broad overview of all published literature examining childhood maltreatment in relation to DNAm has also demonstrated that there is often a divergence between existing studies in terms of findings, measurements of adversity and DNAm, and measurement of and adjustment for environmental confounding [120]. Even for genes such as NR3C1, where there is greater consensus among previous studies, DNAm of NR3C1 has not been associated with maltreatment at the genome-wide level of significance in any existing epigenome-wide association study (EWAS) published to date [120]. As the review authors note, one major source of differences between published studies on the topic is the use of high-risk versus population-based samples. For example, signficant associations between childhood maltreatment and DNAm are typically found in studies carried out in high-risk samples (e.g., children with substantiated histories of maltreatment, socially disadvantaged youth and psychiatric inpatients), while studies using population-based samples have difficulty replicating such results. This could either be because population-based studies may underrepresent more severe maltreatment cases due to lower enrollment of at-risk populations, or perhaps because population-based samples may be more robust to false positives due to considerably larger sample sizes and frequent use of replication samples. Previously reported associations between NR3C1 and FKBP5 DNAm and childhood adversity have also tended to be reported among high-risk populations such as trauma cohorts [72] and holocaust survivors [90], whereas population-based studies such as ours and the Klinger-König et al. study [119] found no significant associations. It is clear that greater attention must be paid to this potential source of sampling variation and bias.

Importantly, these results do not negate the idea that these HPA axis genes may play important roles in disease pathogenesis or are significant biosensors of the external environment, but suggest instead that either: DNAm of these genes is not a reliable marker of the specific stress exposures we studied (but may be for other exposures); DNAm of these genes is not associated with these stress exposures in the general population but only among very high-risk groups; or analyzing DNAm in blood from late adulthood in relation to retrospectively assessed childhood adversity is not an adequate study design to capture these effects, as DNAm of these genes may be too dynamic over time. Indeed, with an average age of 56.3 years in this study, DNAm patterns of these genes would have to have been altered and then remained stable for 30–50 years in order to detect an association. Previous literature supports this idea, with recent candidate gene studies conducted among children showing that adversity in childhood is associated with DNAm of NR3C1 when measured close in time to the exposure in childhood, but that NR3C1 DNAm in those participants is then dynamic and subject to subsequent change [47].

The previously published evidence considered alongside our current analyses suggests that if there is an effect between stress and DNAm of HPA axis genes, it may be only significant and reliably detected among high-risk subgroups of the general population or close in time to the actual maltreatment or stress exposure. Detecting such an effect in the general population may be extremely sensitive to the timing of DNAm assessment after adversity exposure, due to the dynamic nature of DNAm in HPA axis genes over time. Perhaps it is only in very high-risk populations that stress or adversity lead to persistent HPA axis epigenetic dysregulation that can be measured late into adulthood.

The initial results prior to FDR correction suggest that positive religious coping (using one’s religion or spirituality to cope with stressful life situations) may modify the biological impact of childhood financial hardship by attenuating NR3C1 DNAm. This is the first study to assess the impact of R/S measures on the relationship between sources of psychosocial stress and DNAm, and suggests that certain sources of adversity and stress may only have a lasting biological impact in those without resources for resiliency to deal with these sources of stress. Further research is needed to better understand the range of psychosocial stressors affected by R/S and how strongly resources for resilience such as R/S truly modify the biological impact of these stressors. To date, there are no epigenetic studies assessing the effect modification role of R/S on DNAm directly through an interactive effect. Further research is needed to better understand the range of influences that R/S have on human health across religious traditions, using a broad range of R/S measures and across different ethnic communities. Most importantly, future research is also needed to understand how these patterns of DNAm associated with psychosocial stress and spirituality ultimately affect risk of disease, and to understand whether there are particular chronic diseases most affected by certain stressors or dimensions of R/S.

This study has several limitations that should be noted. The results are not necessarily generalizable to all Black women in the USA, as women in the BWHS tend to have slightly higher levels of education than Black women nationally. Although the broad range of stressors analyzed concurrently led to a greater multiple testing penalty imposed on the results, we believe that the approach of analyzing all adulthood and childhood stressors (that we had access to) concurrently is an important way to determine what is important for DNAm. This approach will be helpful especially when cohort study investigators are trying to select measures of psychosocial stress to include in their data collection efforts, as comparative studies such as these will help in prioritizing key measures of stress to include in future research. Furthermore, despite analyzing ten different measures of psychosocial stress and eight different measures of R/S, these measures do not represent the full array of psychosocial stressors or religious and spiritual practices, experiences and beliefs important to Black women in the USA. Similarly, we focused on one biological pathway: dysregulation of the HPA axis as captured in DNAm. There are other important biological pathways through which social adversity and resources for resiliency such as R/S influence disease etiology (e.g., the inflammation pathway [121,122] and immune function [123]). Future research should explore these additional pathways. Conducting untargeted scans of the entire human epigenome to discover potentially novel biological systems that are dysregulated by stress in childhood and adulthood will be crucial to advancing knowledge about psychosocial stress, resiliency and DNAm. Finally, in order for epigenomics research to maximally contribute to understanding how health disparities are produced, future studies should examine a broader range of R/S measures and assess them within diverse minority communities that experience a disproportionate burden of chronic illness. Specifically, future work is needed that analyzes the entire stress–DNAm–disease risk pathway, and then robustly assess the potential buffering effects of different dimensions of R/S and other resources for resilience in attenuating the adverse impact of such stressors (or exacerbating the adverse impact of stressors).

Despite these limitations, this study represents the most comprehensive study to date of the impact of varied and diverse psychosocial stressors on DNAm of HPA axis genes in a US sample of Black women. This study also represents the first investigation of any R/S measures (beyond yoga or meditation) on DNAm, and the first study on R/S and DNAm to include such a large sample of Black women, and to assess psychosocial stress, DNAm and R/S in any population.

Conclusion

When testing direct associations between ten different measures of psychosocial stress and DNAm of HPA axis genes, no consistent and reliable DNAm signals could not be detected across multiple CpGs within the same regions in our sample of 1000 BWHS participants, and no associations survived FDR correction. Further interaction analysis between all combinations of stress and R/S variables with mean DNAm of all four genes as the dependent variables demonstrated that positive religious coping (using one’s religion or spirituality to cope with stressful life situations) may modify the biological impact of childhood financial hardship by attenuating NR3C1 DNAm, and suggests that certain sources of adversity and stress may only have a lasting biological impact in those without resources for resiliency to deal with these sources of stress. Future work investigating these dynamics will need to pay careful attention to the choice of study population and timing of DNAm measurement relative to adversity exposure, and should investigate the impact of resources for resiliency and R/S on the relationship between stress and DNAm.

Future perspective

Given that virtually every stressor studied here is present in communities of color and low-income communities at a much higher rate than in white and more affluent communities, and that such stressors are associated with increased risk of developing a myriad illnesses, efforts to better understand the pathways and mechanism through which such social adversity becomes embodied should remain an urgent public health priority. Efforts to eliminate disparities in health will not be successful without attending to the early formative experiences of trauma and ongoing burden of multiple psychosocial stressors that people of color cope with every day, and attending to cultural resources that may effectively mitigate these adverse impacts. Only by addressing these manifestations of structural inequality and supporting the mobilization of resources for resilience within underserved and at-risk ethnic communities will we eliminate health disparities.

Identifying and developing stable biomarkers of stress and adversity will be a crucial step in preventative clinical efforts to combat chronic disease. Future research should explore whether DNAm of other genes/pathways or different biological mechanisms altogether, are stable biomarkers of stress and adversity in the general population that accurately predict risk of different chronic diseases. These results also emphasize the importance of studying not only psychosocial stress but also resources for resiliency. These resources for resiliency differ across different cultural communities, which may hold different R/S beliefs or practice different religious traditions. Better understanding the role of religiosity and spirituality in shaping the contours of human health may prove a valuable, understudied leverage point in the quest to reduce and eliminate health disparities.

Summary points

• Among 1000 participants from the Black Women’s Health Study, associations between ten psychosocial stressors and DNA methylation (DNAm) levels we tested across 15 different loci within four stress-related genes (NR3C1, HSDB1, HSD11B2 and FKBP5).

• Several stressors were significantly associated with DNAm of individual CpG loci in initial analyses, but no associations remained significant after correcting for multiple testing by false discovery rate (FDR).

• Importantly, we failed to replicate previously reported associations between childhood abuse and DNAm at loci within FKBP5.

• The results indicate that DNAm of genes in the stress pathway (i.e., HPA axis genes) may not be a stable biomarker of childhood and adulthood stress in a retrospective study design using population-based cohort study data.

• Several dimensions of religiosity or spirituality (R/S) modified the relationship between some stressors (lifetime discrimination and financial hardship in childhood) and DNAm, although these interaction analyses became nonsignificant after the FDR correction.

• Further research investigating the role of religion and spirituality as important resources for resiliency within the Black community may prove to be necessary parts of the picture when analyzing relationships between stress and methylation of genes in the stress pathway.

• Future work should also interrogate whether other biological mechanisms or DNAm of other genes/pathways may be more stable biomarkers of stress and adversity in general populations.

Ethical conduct of research

The authors state that they have obtained appropriate institutional review board approval or have followed the principles outlined in the Declaration of Helsinki for all human or animal experimental investigations. In addition, for investigations involving human subjects, informed consent has been obtained from the participants involved.

Acknowledgments

The authors wish to thank N Vicas for superb research assistance. The authors also wish to thank the Dana-Farber/Harvard Cancer Center in Boston, MA, for the use of the Genotyping and Genetics for Population Sciences Core, which provided TaqMan OpenArray genotyping service.

Supplementary data

To view the supplementary data that accompany this paper please visit the journal website at: www.tandfonline.com/doi/suppl/10.2217/epi-2021-0275

Financial & competing interests disclosure

This study on stress, spirituality and health is supported by The John Templeton Foundation (grant nos. 48424 and 59607). The Black Women’s Health Study is supported by R01 CA 058420 and U01 CA 164974. Dana-Farber/Harvard Cancer Center is supported in part by an NCI Cancer Center Support (grant no. National Institute of Health 5 P30 CA06516). The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed.

No writing assistance was utilized in the production of this manuscript.

Additional information

Funding

This study on stress, spirituality and health is supported by The John Templeton Foundation (grant nos. 48424 and 59607). The Black Women’s Health Study is supported by R01 CA 058420 and U01 CA 164974. Dana-Farber/Harvard Cancer Center is supported in part by an NCI Cancer Center Support (grant no. National Institute of Health 5 P30 CA06516). The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed. No writing assistance was utilized in the production of this manuscript.