There is little doubt that stress plays a role in the etiology of depression, but a relative minority of people who face negative life events, even severe events, develop depression (Hammen, Reference Hammen2005). As such, considerable research has focused on better understanding individual differences in sensitivity to depression following stress exposure. Central to this work is the stress sensitization model, which proposes that exposure to life stress during a developmentally sensitive period can increase sensitivity to future stress, thereby increasing depression risk (e.g., Post, Reference Post2016). In essence, the stress sensitization model suggests that early stress exposure alters the effects of later stress exposure on depression, in what could be described as an “environment-by-environment” (E × E) interaction.
Consistent with the stress sensitization model, research suggests that childhood adversity (CA; e.g., maltreatment, parental loss, parental separation) increases individuals’ sensitivity to later proximal stress, thereby increasing risk for depressive symptoms and major depressive disorder (for a review, see Stroud, Reference Stroud, Harkness and Hayden2019). For example, prior work suggests that the link between proximal major stress and depression is stronger among those with a history of CA, as compared to those without such history (McLaughlin, Conron, Koenen, & Gilman, Reference McLaughlin, Conron, Koenen and Gilman2010). Further, adolescents who have experienced CA are more likely to develop depression in the face of lower levels of recent stress (vs. those without CA exposure; e.g., Harkness, Bruce, & Lumley, Reference Harkness, Bruce and Lumley2006; La Rocque, Harkness, & Bagby, Reference La Rocque, Harkness and Bagby2014; Rudolph & Flynn, Reference Rudolph and Flynn2007), suggesting that CA increases stress sensitivity. Notably, among adolescents, the sensitizing effect of CA has been supported across different types of CA, including childhood maltreatment (e.g., Harkness et al., Reference Harkness, Bruce and Lumley2006; La Rocque et al., Reference La Rocque, Harkness and Bagby2014), parental separation and divorce, marital conflict, and caregiver chronic illness (Espejo et al., Reference Espejo, Hammen, Connolly, Brennan, Najman and Bor2007; Rudolph & Flynn, Reference Rudolph and Flynn2007; Starr et al., Reference Starr, Dienes, Stroud, Shaw, Li, Mlawer and Huang2017), suggesting that even relatively less severe forms of CA have the capacity to shape sensitivity to future stress. In contrast, some evidence suggests that exposure to mild CAs may produce a protective “stress inoculation” or “steeling” effect (Seery, Holman, & Silver, Reference Seery, Holman and Silver2010), reducing the association between proximal stressors and depression (although evidence for this effect in youth has been somewhat mixed; Rudolph & Flynn, Reference Rudolph and Flynn2007; Shapero et al., Reference Shapero, Black, Liu, Klugman, Bender, Abramson and Alloy2014). These contrasting findings suggest there may be individual differences in how CA influences the stress response.
Genetic variation in stress sensitization
Some individuals who experience CA may be more susceptible to stress sensitization than others as a function of genetic variation. Substantial work documents gene-by-environment interactions (G × Es), in which genetic variants increase risk of depression following exposure to CA or to more recent stressors, respectively (Assary, Vincent, Keers, & Pluess, Reference Assary, Vincent, Keers and Pluess2018). However, it is also possible that some G × E models can be further qualified, given the interactive effects of stress exposure across different points in development. As such, researchers have proposed a more complex model whereby genetic variation moderates the stress sensitization effects of CA—a “gene-by-environment-by-environment” interaction, or “G × E × E” (Daskalakis, Bagot, Parker, Vinkers, & de Kloet, Reference Daskalakis, Bagot, Parker, Vinkers and de Kloet2013; Homberg & van den Hove, Reference Homberg and van den Hove2012; Keers & Pluess, Reference Keers and Pluess2017; Starr, Hammen, Conway, Raposa, & Brennan, Reference Starr, Hammen, Conway, Raposa and Brennan2014), where the first “E” refers to CA and the second “E” refers to proximal stressors that trigger depressive onset (e.g., recent episodic stress).
There are several plausible biological mechanisms that may explain why genetic risk moderates stress sensitization. For example, CA exposure alters the development of stress regulation-related systems and structures (e.g., hypothalamic–pituitary–adrenal [HPA] axis sensitization, reduced hippocampal volume; Cicchetti & Rogosch, Reference Cicchetti and Rogosch2012; Heim, Newport, Mletzko, Miller, & Nemeroff, Reference Heim, Newport, Mletzko, Miller and Nemeroff2008). Genes may influence variation in plasticity of these systems, making some youth more vulnerable to lasting alteration in stress-related neural circuitry following CA exposure (Heim et al., Reference Heim, Newport, Mletzko, Miller and Nemeroff2008). Genetic variation likely also influences the effects of CA on HPA axis functioning and downstream effects on related neurological structures (Gillespie, Phifer, Bradley, & Ressler, Reference Gillespie, Phifer, Bradley and Ressler2009). These in turn may have lasting effects on stress response (e.g., Gunnar & Fisher, Reference Gunnar and Fisher2006), leading youth to be continually at elevated risk of depression following stress exposure. Moreover, differential susceptibility and vantage sensitivity models (Belsky, Bakermans-Kranenburg, & van IJzendoorn, Reference Belsky, Bakermans-Kranenburg and van IJzendoorn2007; Pluess, Reference Pluess2017) suggest that the same genetic variants that elevate risk under adverse conditions may predict positive outcomes under positive environmental conditions; as such, for youth raised in nurturing, low CA environments, genetic “risk” may predict positive coping and other resiliency factors that in turn buffer against depression following stress exposure.
Supporting these ideas, existing research indicates that genetic variation moderates the effect of early adversity on sensitivity to proximal stress, in a G × E × E effect, mostly using serotonergic polymorphism serotonin transporter linked polymorphic region (5-HTTLPR) (e.g., Grabe et al., Reference Grabe, Schwahn, Mahler, Schulz, Spitzer, Fenske and Freyberger2012; Keers & Pluess, Reference Keers and Pluess2017; Kumsta et al., Reference Kumsta, Stevens, Brookes, Schlotz, Castle, Beckett and Sonuga-Barke2010; Starr et al., Reference Starr, Hammen, Conway, Raposa and Brennan2014; for a nonsupportive example, see Power et al., Reference Power, Lecky-Thompson, Fisher, Cohen-Woods, Hosang, Uher and McGuffin2013) and, in one study, using a variant related to HPA axis functioning (corticotropin releasing hormone receptor 1 (CRHR1) rs110402; Starr et al., Reference Starr, Hammen, Conway, Raposa and Brennan2014).
Beyond single candidate genes: Multilocus approaches
Importantly, nearly all of existing G × E × E studies have used single-variant designs (cf., Keers & Pluess, Reference Keers and Pluess2017). Within the G × E literature, the single-variant candidate gene approach has grown increasingly controversial following a mix of positive and negative meta-analytic findings (Bleys, Luyten, Soenens, & Claes, Reference Bleys, Luyten, Soenens and Claes2018; Culverhouse et al., Reference Culverhouse, Saccone, Horton, Ma, Anstey, Banaschewski and Bierut2017; Karg, Burmeister, Shedden, & Sen, Reference Karg, Burmeister, Shedden and Sen2011; Munafò et al., Reference Munafò, Freimer, Ng, Ophoff, Veijola, Miettunen and Flint2009; Risch et al., Reference Risch, Herrell, Lehner, Liang, Eaves, Hoh and Merikangas2009; Sharpley, Palanisamy, Glyde, Dillingham, & Agnew, Reference Sharpley, Palanisamy, Glyde, Dillingham and Agnew2014). The literature has also weathered prominent criticisms, including for example, use of potentially underpowered samples, selective focus on positive findings, inadequate control of covariates, and insufficient attention to statistical issues related to use of dichotomous genetic variables (de Vries, Roest, Franzen, Munafò, & Bastiaansen, Reference de Vries, Roest, Franzen, Munafò and Bastiaansen2016; Dick et al., Reference Dick, Agrawal, Keller, Adkins, Aliev, Monroe and Sher2015; Duncan & Keller, Reference Duncan and Keller2011). These criticisms, however, have prompted important counterarguments (Caspi, Hariri, Holmes, Uher, & Moffitt, Reference Caspi, Hariri, Holmes, Uher and Moffitt2010; Karg et al., Reference Karg, Burmeister, Shedden and Sen2011; Moffitt & Caspi, Reference Moffitt and Caspi2014; Vrshek-Schallhorn, Corneau, & Starr, Reference Vrshek-Schallhorn, Corneau and Starr2019; Vrshek-Schallhorn, Sapuram, & Avery, Reference Vrshek-Schallhorn, Sapuram and Avery2017). Moreover, fundamental to dominant genetic theoretical models is that the genetic architecture of complex diseases like depression is polygenic, with many genes each contributing small effects in a primarily additive manner (e.g., Sullivan, Daly, & O'Donovan, Reference Sullivan, Daly and O'Donovan2012). Consequently, it is likely difficult to model effect sizes associated with single variants without using very large samples, a design often requiring reliance on low-cost measurements of stress exposure, introducing threats to validity and weakening G × E effects (Karg et al., Reference Karg, Burmeister, Shedden and Sen2011; Monroe, Reference Monroe2008). This problem is likely compounded when examining G × E × E effects, as three-way interactions are underpowered compared to two-way interactions and main effects (Heo & Leon, Reference Heo and Leon2010).
To address these challenges, researchers have devised multilocus genetic profile scores (MGPSs; Bogdan, Hyde, & Hariri, Reference Bogdan, Hyde and Hariri2013). MGPSs are unweighted, additive scores of risk alleles drawn from genetic loci associated with biological systems of interest. MGPSs capture predominantly functional (i.e., biologically meaningful) genetic variation compatible with theory-driven research. They also incorporate multiple variants that yield continuous data, enhancing statistical power (Aliev, Latendresse, Bacanu, Neale, & Dick, Reference Aliev, Latendresse, Bacanu, Neale and Dick2014; Caspi et al., Reference Caspi, Hariri, Holmes, Uher and Moffitt2010; Dick et al., Reference Dick, Agrawal, Keller, Adkins, Aliev, Monroe and Sher2015) and improving effect sizes (vs. those found in traditional single-variant designs). For example, one study found that a dopaminergic MGPS accounted for 11% of reward-related ventral striatal activation (Nikolova, Ferrell, Manuck, & Hariri, Reference Nikolova, Ferrell, Manuck and Hariri2011). Another found that an HPA axis MGPS G × E accounted for 8% of depression variance among adolescents (Starr & Huang, Reference Starr and Huang2019), which is 80 times the presumed G × E effect size for single-variant designs of 0.1% (Duncan & Keller, Reference Duncan and Keller2011).
To date, only one study has tested G × E × E using a polygenic approach. Keers and Pluess (Reference Keers and Pluess2017) created an “environmental sensitivity” genetic risk index, which included 5-HTTLPR and nine other variants located on genes in multiple systems (e.g., brain-derived neurotrophic factor (BDNF), CRHR1, tryptophan hydroxylase 2 (TPH2)). In a birth cohort sample of over 7,000 participants followed over 50 years, the authors found that greater risk score moderated the effect of the quality of childhood material environment (an index of socioeconomic status and familial financial hardship) on sensitivity to material environment in adulthood. Results specifically supported a “for better or for worse” pattern (i.e., differential susceptibility; Belsky et al., Reference Belsky, Bakermans-Kranenburg and van IJzendoorn2007), where genetic “risk” conferred vulnerability at negative environmental conditions, but predicted thriving under positive conditions. This innovative study showed longitudinal support for a polygenic G × E × E model, but the environmental measures were limited. Further, the authors used a risk index capturing genes in multiple systems, but MGPSs capturing genetic variation from a single biological system may reveal insights into biological pathways that function as intermediate phenotypes.
In an attempt to replicate and extend previous single-variant G × E × E findings using the multilocus approach (Starr et al., Reference Starr, Hammen, Conway, Raposa and Brennan2014), we focus on MGPSs capturing genetic variance from two specific biological systems: the HPA axis and serotonin systems.
HPA axis genetic variation
Research suggests that CA exposure predicts alterations in HPA axis regulation, as indexed by basal cortisol, cortisol reactivity, diurnal cortisol rhythms, and trait cortisol (Fries, Shirtcliff, & Pollak, Reference Fries, Shirtcliff and Pollak2008; Gonzalez, Jenkins, Steiner, & Fleming, Reference Gonzalez, Jenkins, Steiner and Fleming2009; Gunnar & Quevedo, Reference Gunnar and Quevedo2007; Stroud, Chen, Doane, & Granger, Reference Stroud, Chen, Doane and Granger2016). Such alterations, in turn, are linked to depression (e.g., Chida & Steptoe, Reference Chida and Steptoe2009; Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Doane, Mineka, Zinbarg, Craske and Adam2013), as well as to elevated rates of depression following recent stress exposure (Schuler et al., Reference Schuler, Ruggero, Goldstein, Perlman, Klein and Kotov2017). Individual single nucleotide polymorphisms (SNPs) linked to HPA axis-related genes (e.g., CRHR1, FK506 binding protein 5 (FKBP5)) appear to moderate the association between stress exposure (both childhood and proximal) and depression (e.g., Bradley et al., Reference Bradley, Binder, Epstein, Tang, Nair, Liu and Newport2008; Vinkers et al., Reference Vinkers, Joëls, Milaneschi, Gerritsen, Kahn, Penninx and Boks2015; Zimmermann et al., Reference Zimmermann, Brückl, Nocon, Pfister, Binder, Uhr and Ising2011), and Starr et al. (Reference Starr, Hammen, Conway, Raposa and Brennan2014) found support for a G × E × E effect predicting depression with CRHR1 SNP rs110402.
More recent efforts using the MGPS approach have shown that a 10-SNP index, constructed with SNPs selected for their association with depression and related phenotypes or HPA axis dysregulation and located in CRHR1, nuclear receptor subfamily 3, group C, member 1 (NR3C1), nuclear receptor subfamily 3 group C member 2 (NR3C2), and FKBP5 (Pagliaccio et al., Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden and Barch2014), interacts with life stress among youth to predict hippocampal amygdala volume and connectivity, higher cortisol reactivity, and flatter diurnal slope following recent proximal chronic stress exposure (Pagliaccio et al., Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden and Barch2014, Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden and Barch2015; Starr, Dienes, Li, & Shaw, Reference Starr, Dienes, Li and Shaw2019a). Findings also extend to affective outcomes, with this MGPS interacting with both CA and recent proximal stress to predict depression (Starr & Huang, Reference Starr and Huang2019). Other overlapping HPA axis MGPSs also moderate the effects of early adversity and proximal stress to predict affective outcomes (Di Iorio et al., Reference Di Iorio, Carey, Michalski, Corral-Frias, Conley, Hariri and Bogdan2017; Feurer et al., Reference Feurer, McGeary, Knopik, Brick, Palmer and Gibb2017). However, no prior work has tested a G × E × E model using multilocus HPA axis genetic variation.
Serotonergic genetic variation
Interest in serotonergic genetic risk spiked following the breakthrough finding of Caspi et al. (Reference Caspi, Sugden, Moffitt, Taylor, Craig, Harrington and Poulton2003) that a polymorphism in the serotonin transporter gene (5-HTTLPR) moderates the association between both proximal and childhood stress exposure, and depression. As noted, multiple studies support 5-HTTLPR's role as a moderator of stress sensitization (Grabe et al., Reference Grabe, Schwahn, Mahler, Schulz, Spitzer, Fenske and Freyberger2012; Kumsta et al., Reference Kumsta, Stevens, Brookes, Schlotz, Castle, Beckett and Sonuga-Barke2010; Starr et al., Reference Starr, Hammen, Conway, Raposa and Brennan2014). However, G × E research on 5-HTTLPR has been famously unreliable (Culverhouse et al., Reference Culverhouse, Saccone, Horton, Ma, Anstey, Banaschewski and Bierut2017; Munafò, Durrant, Lewis, & Flint, Reference Munafò, Durrant, Lewis and Flint2009) though some meta-analytic evidence suggests significant G × E effects when high quality stress assessment methods were used and CA was examined (Karg et al., Reference Karg, Burmeister, Shedden and Sen2011). Using the MPGS approach, Vrshek-Schallhorn et al. (Reference Vrshek-Schallhorn, Stroud, Mineka, Zinbarg, Adam, Redei and Craske2015b) showed that a five-SNP multilocus profile capturing serotonergic genetic variation interacted with recent stressful life events to predict subsequent depression in two samples, a finding recently replicated in a third sample (Starr, Vrshek-Schallhorn, & Stroud, Reference Starr, Vrshek-Schallhorn and Stroud2019b). However, no studies have examined this MGPS (or any serotonergic MGPS) in the context of CA or stress sensitization.
Additional considerations of the environment and development
Some researchers have observed that careful definition of the “candidate environment” in serotonergic G × E models is an often-overlooked step (e.g., Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Stroud, Mineka, Hammen, Zinbarg, Wolitzky-Taylor and Craske2015). Environmental stress is an extraordinarily heterogeneous construct, and certain biological systems—and corresponding genetic variation—may respond selectively to certain kinds of environmental threat. For example, major interpersonal stressful life events have been supported as a candidate environment for serotonergic genetic risk for depression, when risk is indexed with either 5-HTTPLR or a serotonergic MGPS (Starr, Vrshek-Schallhorn, & Stroud, Reference Starr, Vrshek-Schallhorn and Stroud2019b; Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Mineka, Zinbarg, Craske, Griffith, Sutton and Adam2014; Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Stroud, Mineka, Zinbarg, Adam, Redei and Craske2015). Interpersonal stress is an especially potent risk factor for depression, especially for more severe events (Brown, Bifulco, & Harris, Reference Brown, Bifulco and Harris1987; Hammen, Reference Hammen2005; Sheets & Craighead, Reference Sheets and Craighead2014; Stroud, Davila, Hammen, & Vrshek-Schallhorn, Reference Stroud, Davila, Hammen and Vrshek-Schallhorn2011; Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Stroud, Mineka, Hammen, Zinbarg, Wolitzky-Taylor and Craske2015), and some research supports a role for the serotonin system in regulating sensitivity to social events (e.g., Way & Taylor, Reference Way and Taylor2010). Indeed, in the primary sample used here, we previously demonstrated that adolescents with high serotonergic MGPS showed intensified associations between major interpersonal stressors, but not major noninterpersonal or minor stressors, and depression (Starr, Vrshek-Schallhorn, & Stroud, Reference Starr, Vrshek-Schallhorn and Stroud2019b).
If youth at high serotonergic genetic risk are more sensitive to major interpersonal stressors, this effect may be amplified by CA exposure, as stress sensitization may increase their reactivity to these salient stressors (i.e., selective sensitization; Slavich, Monroe, & Gotlib, Reference Slavich, Monroe and Gotlib2011). In contrast, adolescents with high serotonergic MGPS may be less attuned to proximal stressors that are noninterpersonal in nature, and therefore stress sensitization may be less robust in response to these experiences. Thus, as an a priori choice, we specifically tested major interpersonal events when testing G × E × E effects for the serotonergic MGPS (while also testing overall episodic stress to probe generalizability to stress broadly). In contrast, based on prior work showing that the HPA axis MGPS moderates a broad range of proximal stressors, including chronic and acute, and interpersonal and noninterpersonal (Starr & Huang, Reference Starr and Huang2019), for the HPA axis MGPS we examined overall episodic stress as an index of stress exposure.
Adolescence has been conceptualized as a second “sensitive period” wherein developmental changes in neuroendocrine physiology increase sensitivity to the environment (e.g., Del Giudice, Ellis, & Shirtcliff, Reference Del Giudice, Ellis and Shirtcliff2011; Viner et al., Reference Viner, Ozer, Denny, Marmot, Resnick, Fatusi and Currie2012). Adolescence is characterized by increased exposure to and reactivity to environmental stress, higher stress sensitization effects among CA-exposed youth (La Rocque et al., Reference La Rocque, Harkness and Bagby2014; Rudolph & Hammen, Reference Rudolph and Hammen1999), and increased risk for first onsets of depression (Avenevoli, Swendsen, He, Burstein, & Merikangas, Reference Avenevoli, Swendsen, He, Burstein and Merikangas2015). As such, it is important to investigate G × E × E effects during this developmental period.
The present study
We hypothesized two three-way interactions predicting depression: (a) HPA axis MGPS, CA, and recent episodic stress, and (b) serotonergic MGPS, CA, and major recent interpersonal episodic stress. For both MPGSs, we expected that as genetic risk scores increased, stress sensitization (as captured by the two-way interaction between CA and proximal stress) would increase, such that those high on both MGPS and CA would show the strongest associations between episodic stress and depression. Both hypotheses were tested in a community sample of mid-adolescents. The serotonergic MGPS hypothesis was also tested in a replication sample of early adolescent girls. Notably, clinically significant first onsets of depression are just beginning to occur in this developmental period, with relatively low rates of diagnosable depression in community samples (Lewinsohn, Hops, Roberts, & Seeley, Reference Lewinsohn, Hops, Roberts and Seeley1993). Subthreshold symptoms are more prevalent in this age group and are predictive of later onset of full-syndrome depression (Judd et al., Reference Judd, Akiskal, Maser, Zeller, Endicott, Coryell and Keller1998; Shankman et al., Reference Shankman, Lewinsohn, Klein, Small, Seeley and Altman2009). As such, and because of increased power conferred by continuous variables, depression was coded on a dimensional scale that captured both diagnosable and subthreshold symptoms.
Methods: Primary sample
Participants
The full sample included 241 adolescents aged 14–17 years (M age = 15.90, SD = 1.09) who participated with their primary caregiver, and were recruited from a mid-sized metropolitan area using multiple recruitment methods, including advertisements, a commercial mailing list, and ResearchMatch, a national research health volunteer registry (for more recruitment details, see Starr et al., Reference Starr, Dienes, Stroud, Shaw, Li, Mlawer and Huang2017). Exclusion criteria included prior diagnosis of bipolar or psychotic disorder; any major physical, neurological, or pervasive developmental disorder; English reading or language difficulties; and prior participation of another household member. The median parent-reported annual family income was $80,000–$89,999 (approximately 60th percentile for the state), with 24.1% of youth receiving free or reduced-cost lunch at school (an index of economic hardship). Most participating parents were biological mothers (87.6%). To account for population stratification, most analyses were conducted in adolescents reporting European ancestry (adolescent race was parent-reported; participant reports of racial identification correspond closely with ancestry informative markers; Sinha, Larkin, Elston, & Redline, Reference Sinha, Larkin, Elston and Redline2006). The sample included 192 White adolescents (M age = 15.89 years, SD = 1.08, 53.1% femaleFootnote 1). Power was estimated in G*Power (Faul, Erdfelder, Lang, & Buchner, Reference Faul, Erdfelder, Lang and Buchner2009) using parameters derived from single-SNP G × E × E research (Starr et al., Reference Starr, Hammen, Conway, Raposa and Brennan2014), with interaction R 2 estimates of .02 for HPA axis genetic risk and .01 for serotonergic genetic risk. Notably, this is a conservative estimate as the MGPS approach very likely confers additional power due to additional variants increasing effect size and dimensional approaches boosting power relative to dichotomous ones. Power estimates ranged from .50 (White sample) to .60 (full sample) for HPA axis genetic risk and .28 to .34 for serotonergic risk. Although these estimates are likely conservative, they suggest findings should be interpreted with caution, and in particular, the serotonergic MGPS findings should be considered in conjunction with the replication study results.
Measures
Genotyping and MGPS calculation
DNA was analyzed at the University of Wisconsin-Madison Biotechnology Center. Following extraction with standard salting-out procedures, concentration was measured using a PicoGreen® kit (Life Technologies, Grand Island, NY), and genotyping used KBiosciences competitive allele-specific polymerase chain reaction (PCR) SNP genotyping assay based on dual FRET (KASPar) with endpoint fluorescence detection (additional details available on request).
Following Pagliaccio et al.'s (Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden and Barch2014), the HPA axis MGPS included risk alleles of 10 SNPs from four HPA axis-related genes: CRHR1 (rs4792887 T-allele, rs110402 G-allele, rs242941 T-allele, rs242939 G-allele, rs1876828 G-allele), NR3C1 (rs41423247 G-allele, rs10482605 T-allele, rs10052957 A-allele), NR3C2 (rs5522 G-allele), and FKB5 (rs1360780 T-allele). These SNPs were chosen for associations with depression and related phenotypes and/or dysregulated cortisol and were pruned from a list of 15 to reduce linkage disequilibrium (Pagliaccio et al., Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden and Barch2014). Distributions of genotypes are available upon request; all were in Hardy–Weinberg equilibrium (HWE) with the exception of rs1876828, χ2 (1) = 4.12, p = .041, although this SNP was in HWE when tests were stratified by racial group, suggesting the HWE failure was due to racial admixture. Re-running analyses excluding this SNP did not impact results. HPA axis genetic profile scores were computed by coding SNPs for presence (1) or absence (0) of at-risk genotypes (.5 codes were assigned to hetero-zygotes in cases of allelic rather than genotypic effects). Individual SNP codes were summed, with higher MGPS reflecting higher genetic risk (possible range 0–10, sample range 2–9). We permitted up to two missing genotypes (20%) per person, rescaling using available SNPs when necessary (no participants were missing >2 SNPs). Based on the prior literature, higher scores can be interpreted as increased liability to HPA system dysregulation.
In accordance with Vrshek-Schallhorn et al. (Reference Vrshek-Schallhorn, Stroud, Mineka, Zinbarg, Adam, Redei and Craske2015b), genotype data for the serotonergic MGPS were obtained from four SNPs selected for their functional effects (causal biological impact) (5-hydroxytryptamine receptor (HTR)1A rs6295 G-allele in forward coding direction, HTR2A rs6314 C-allele, HTR2C rs6318 C-allele, TPH2 rs11178997 T-allele) and one SNP implicated in depression risk via meta-analysis (rs4570625 G-allele), located on or near the four genes linked to serotonergic functioning. Risk direction was coded based on evidence of association with depression and related outcomes. All SNPs were in HWE. Serotonergic genetic profile scores were computed by coding SNPs based on the number of risk alleles (0–2), except for rs6318 which is X-linked (coded as C-allele presence [1] or absence [0]). Scores were summed to create an additive index reflecting number of risk alleles (possible range 0–9, sample range 2–9). We permitted missing data for one SNP (20%), rescaling the MGPS using available data (no participants were missing > 1 SNP).
Proximal stress
To assess past year episodic stress (i.e., discrete events with a brief onset and relatively short duration), adolescents were administered the UCLA Life Stress Interview (LSI) (Rudolph & Hammen, Reference Rudolph and Hammen1999; Rudolph et al., Reference Rudolph, Hammen, Burge, Lindberg, Herzberg and Daley2000). For each event reported, adolescents provided information about the surrounding context (e.g., circumstances and resources to cope with it, predictability, and prior experience with similar events), its duration, and the consequences to obtain the degree of objective impact. Interviewers prepared narrative accounts of each event (including the surrounding circumstances and consequences, but excluding participants’ subjective reactions), which were presented to an independent rating team blind to all other study data. Consistent with prior work (Rudolph et al., Reference Rudolph, Hammen, Burge, Lindberg, Herzberg and Daley2000), for each event, the team rated: (a) the objective impact (rated 1 [no negative impact] to 5 [extremely severe negative impact], with half-points permitted); and (b) interpersonal status (coded 1/0: rated interpersonal when the primary context involved relations with others [e.g., conflicts] or affected the participants’ relations [e.g., friend moving away]). Events rated 2.5+ were considered major events (similar results were obtained when using a more conservative 3.0 threshold). In the full sample, participants reported an average of 2.95 events (SD = 2.08), with 92% of the sample reporting at least one event and 45% reporting at least one major interpersonal event. A second team, blind to the original ratings, rerated a set of events with excellent reliability (ICC = .86). Composite scores were formed by summing the severity ratings for (a) all events, and (b) major interpersonal events.
Supplemental analyses examined chronic stress as an alternative phenotype. The Chronic Stress Interview (CSI), a portion of the LSI, is a semi-structured interview adapted for use with late adolescents (Hammen and Brennan, Reference Hammen and Brennan2001). The CSI evaluates and objectively codes the nature and quality of ongoing conditions in the prior six months adolescents’ lives across multiple domains: close friendships, peer relations, family life, romantic life, academics, and behavioral issues. The CSI isolates objective assessments of ongoing stressful circumstances from the adolescents’ subjective perceptions of stress, with interviewers rating chronic stress based on behaviorally specific objective features of individual domains on a scale from 1 (exceptionally good conditions) to 5 (extreme adversity; half-points permitted). Scores across all domains were averaged. Second raters re-coded 20% of interviews; ICC = .85.
Childhood adversity
The lifetime adversity section of the Youth Life Stress Interview (YLSI) (Rudolph & Flynn, Reference Rudolph and Flynn2007) was used to assess adolescents’ exposure to negative family events and circumstances during their lifetime (up until the year prior to the interview, ensuring no overlap with the time period assessed for CA as compared to recent episodic events and chronic stress). This interview was administered to participating caregivers due to time constraints and better recollection of early childhood events. A general probe was used first to assess exposure to particularly stressful events and circumstances. Next, a series of probes was used assessing specific types of adversity (death of a loved one, parental separations, parental separation/ divorce, exposure to serious marital conflict, physical or mental illness of a loved one, multiple family transitions, chaotic family living circumstances [e.g., neglect], family legal problems, financial difficulties, or other very difficult adversity) was used.
For each adversity endorsed, participants gave information about the context (i.e., circumstances) and the consequences, and interviewers prepared narrative accounts, which excluded participants’ subjective reactions. Narratives were presented to an independent rating team who were blind to participants’ subjective reactions and all other data. The team rated the objective impact (i.e., severity) on a 1 to 5 scale (half-points permitted) for each adversity, and then rated the overall severity of CA (taking into account all individual events, weighted by relative importance) on a scale from 1 (no adversity) to 10 (extreme, severe adversity). The overall severity rating was used as the primary index of CA, consistent with prior research (Chen, Stroud, Vrshek-Schallhorn, Doane, & Granger, Reference Chen, Stroud, Vrshek-Schallhorn, Doane and Granger2017). However, as some prior publications (e.g., Starr & Huang, Reference Starr and Huang2019) have summed the severity scores of individual events, this index was used in supplemental analyses. Across the full sample, participants reported an average of 4.56 events (SD = 2.85), with the vast majority (97%) reporting at least one event, and most participants reporting at least one more severe event (78% with an event with 2.5 + severity and 60% with an event with 3.0 + severity). An independent coding team rerated a set of participants (overall severity rating: ICC = .86; individual adversities: ICC = .97).
Depression
The Schedule for Affective Disorders and Schizophrenia for school-aged children, present and lifetime version (K-SADS-PL; Kaufman et al., Reference Kaufman, Birmaher, Brent, Rau, Flynn, Moreci and Ryan1997) assessed current depressive symptoms. Symptoms were rated: 0 = no symptoms; 1 = mild symptoms; 2 = moderate, subthreshold symptoms; 3 = meets fourth edition of Diagnostic and Statistical Manual of Mental Disorders (DSM-IV) criteria; 4 = DSM-IV criteria met with high severity. The maximum score of either current major depressive disorder or dysthymia was used. Consistent with point prevalence rates in other adolescent community samples (Lewinsohn et al., Reference Lewinsohn, Hops, Roberts and Seeley1993), 15% of the sample reported current clinically significant symptoms (5% meeting diagnostic criteria). Interrater reliability (assessed by having second coders, drawn from the original pool of interviewers, rate audiotaped interviews) was perfect (ICC = 1.0).
Pubertal status for covariation
Pubertal status was assessed via adolescent self-report using the Pubertal Development Scale (PDS; Petersen, Crockett, Richards, & Boxer, Reference Petersen, Crockett, Richards and Boxer1988). Five sex-specific items, assessing growth spurt in height, skin and body hair changes, and facial hair and vocal chord changes in boys and breast development and menarche in girls, are rated on a 4-point scale, from no development (1) to development seems completed (4); except for menarche which is rated 1 or 4. The PDS has shown good reliability and validity (Petersen et al., Reference Petersen, Crockett, Richards and Boxer1988).
Procedure
Participants and a primary caregiver attended a lab session; informed consent/assent was obtained, after which adolescents completed interviews and questionnaires and provided saliva samples using Oragene™ (DNA Genotek, Ontario, Canada) collection kits. Families were paid $160 for participation and entered into raffles. The University of Rochester Research Subjects Review Board approved all procedures.
Primary sample results
All primary sample analyses were conducted using SPSS 25 (IBM Corp.), with interactions tested using the PROCESS macro (Hayes, Reference Hayes2013). Data were complete for all major variables. Descriptive data for all variables are available in Table 1. Figure 1 displays a scatterplot of CA by total episodic severity at different levels of each MGPS.
Figure 1. Scatterplot of childhood adversity by total recent episodic stress in the portion of the primary sample identifying as White, with values for serotonergic multilocus genetic profile score (MGPS) represented by shape (by quartile) and values for HPA axis MGPS represented by shading [in print publication]/ color [in online publication].
Table 1. Descriptive statistics and bivariate correlations among major study variables in primary and replication samples
Note: All data are from full samples unless otherwise noted. Correlations from White subsamples available upon request. *p < .05, **p < .01, ***p < .001.
Preliminary analyses
Because of the potentially confounding effect of population stratification by race in genetic analyses, we tested for differences between self-reported racial groups. For the HPA axis MGPS, non-White youth showed significantly higher MGPS than White youth, t (239) = 2.10, p = .036, and also scored marginally higher on depressive symptoms (t (58.03) = 2.38, p = .073) and total episodic stress (t (52.19) = 1.94, p = .057). Race also moderated the association between HPA axis MGPS and depression (interaction b = .19, SE = .08, p = .024, stronger associations for White youth; see also Starr & Huang, Reference Starr and Huang2019). In contrast for the serotonergic MGPS, self-identified White race is unrelated to MGPS and does not moderate its association with depression (ps > .05; see also Starr et al., Reference Starr, Vrshek-Schallhorn and Stroud2019b). Thus, for the HPA axis MGPS, analyses were restricted to the subsample reporting European ancestry for the primary models. For consistency, we did the same for the serotonergic MGPS within the primary sample. For both MGPSs, we replicated all analyses using the full sample.
Bivariate correlations are presented in Table 1. Results of two-way interactions (e.g., CA × Proximal Stress interactions, MGPS × environment interactions) will not be presented because they are not the present focus, and they have been tested and described elsewhere (Starr et al., Reference Starr, Dienes, Stroud, Shaw, Li, Mlawer and Huang2017; Starr & Huang, Reference Starr and Huang2019; Starr, Vrshek-Schallhorn, & Stroud, Reference Starr, Vrshek-Schallhorn and Stroud2019b).
HPA axis MGPS × proximal episodic stress × CA
We next examined whether HPA axis genetic variation moderated stress sensitization to depression following CA, by testing a G (HPA axis MGPS) × E (total recent episodic severity) × E (CA) three-way interaction, with all predictor variables centered and depression entered as the outcome. The three-way interaction was entered along with all constituent main effects and two-way interactions. The three-way interaction was significant, p = .006, with the G × E × E effect accounting for 3% of the variance in depression (see Table 2). Simple slope tests revealed that the CA × Proximal Stress interaction (capturing stress sensitization) was significant at high MGPS (M + SD), b = .019, SE = .007, p = .008. For these high MGPS adolescents, proximal stress was associated with depression at high but not low CA exposure. For low MGPS adolescents (M-SD), there was no CA × Proximal Stress interaction (b = −.009, SE = .009, p = .332), suggesting no stress sensitization. This interaction is illustrated in Figure 2. A Johnson–Neyman test (JNT) revealed a region of significance beginning at MGPS = 5.37 (74th percentile), meaning that the CA × Proximal Stress interaction predicted depression significantly for those in the top quartile on MGPS of the sample, but not significantly for those below.
Figure 2. Prediction of depression ratings by recent episodic stress (centered) as a function of childhood adversity history at (a) high hypothalamic–pituitary–adrenal (HPA) axis multilocus genetic profile score (MGPS) (M + 1 SD), and (b) low HPA axis MGPS (M − 1 SD).
Table 2. Primary sample: Results of linear regression model testing three-way interaction between total proximal episodic stress, childhood adversity, and HPA axis multilocus genetic profile score (MGPS), predicting depression
Note: Model conducted in the portion of the sample reporting European ancestry.
As an alternative interpretation to this effect, we examined whether the genetic moderation of recent stress exposure was intensified by CA exposure (i.e., an alternative decomposition of the same G × E × E interaction). Supporting this, the G × E was significant at high levels of CA (b = .03, SE = .01, p = .003), with those at high MGPS showing stronger proximal stress–depression associations. However, the G × E was not significant among those with low levels of CA (p = .622), with a JNT region of significance beginning at 3.61 (64th percentile).
As tests of robustness, we ran subsequent models with age, pubertal maturation, and gender as covariates. In these models, the three-way interaction remained significant, when: (a) the main effects of these covariates were included (p = .006); and (b) the two-way G × Covariate, CA × Covariate, and Proximal Stress × Covariate interactions (resulting in nine additional covariates terms) were included (p = .003). Finally, the three-way interaction remained significant when using the full sample, with race (White) and the three race interactive terms (i.e., G × White, CA × White, Proximal Stress × White) as covariates (and the full set of covariates, for a total of 16 covariates), p = .004.
HPA axis MGPS sensitivity analyses
To determine whether multilocus genetic variation (vs. a single dominant SNP) accounted for the three-way interaction, we conducted “n − 1” sensitivity analyses (Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Stroud, Mineka, Zinbarg, Adam, Redei and Craske2015) by constructing 10 modified MGPSs, each excluding one SNP, and then tested the G × E × E model with each modified MGPS. All three-way interaction terms remained significant with the exclusion of individual SNPs (ps < .05) (see Supplementary Table 1).
Serotonergic MGPS × Proximal Episodic Stress × CA
First, following the same steps as outlined above, we tested a three-way interaction between serotonergic MGPS, total episodic stress, and CA in predicting depression, which was not significant, b = .001, SE = .116, p = .486. Next, as planned, we restricted the proximal episodic stress variable to major interpersonal stress. The three-way interaction was significant (p = .000006) and explained 9% of the variance in depression (see Table 3). Simple slope tests indicated that at high levels of MGPS (M + SD), the E × E stress-sensitization effect was significant (b = .05, SE = .01, p = .0004). Specifically, as illustrated in Figure 3, among those at high MGPS, interpersonal major stress was predictive of depression at high CA (p = .0000002) but not low CA (p = .4741), consistent with stress sensitization. In contrast, for adolescents at low MGPS, the E × E effect was also significant, but in the opposite direction, b = −.05, SE = .02, p = .0055 (consistent with a “steeling” or “stress inoculation” effect; Seery et al., Reference Seery, Holman and Silver2010), with reduced associations between proximal stress and depression for those with high CA. JNT analyses revealed significant E × E effects when MGPS ≥ 6.95 (76th percentile) or, in the opposite direction, when MGPS ≤ 5.74 (24th percentile). As an alternative interpretation, we decomposed the three-way interaction by CA level, and found that there was a significant G × E effect at high (b = .10, SE = .02, p = .0000005), but not low (b = −.04, SE = .02, p = .1162), levels of CA. JNT analyses suggested that above CA ratings of 3.13 (64th percentile), MGPS significantly predicted increased associations between major interpersonal stressors and depression, but also that at very low ratings of CA ( ≤ 1.13, 14th percentile), MGPS conferred a protective effect, predicting significantly weaker associations between major interpersonal stress and depression.
Figure 3. Prediction of depression ratings by recent interpersonal major stress (centered) as a function of childhood adversity history, at (a) high serotonergic multilocus genetic profile score (MGPS) (M + 1 SD), and (b) low serotonergic MGPS (M − 1 SD).
Table 3. Results of linear regression models testing three-way interaction between major proximal major interpersonal stress, childhood adversity, and serotonin multilocus genetic profile score, predicting depression in two samples
Note: In both samples, models were conducted in the portion of the sample reporting European ancestry.
In tests of robustness, the interaction term remained significant when accounting for main effects of age, gender, and pubertal status (p = .000008) as well as when including interactive effects of these covariates and the three genetic and environmental variables (p = .000003). Finally, the interaction remained significant in the full sample, with race and race interaction variables as covariates (p = .0173), and when accounting for the full set of covariates (p = .0036; R 2 for the three-way interaction = 3% in the models with covariates).
Serotonergic MGPS sensitivity analyses
In n − 1 analyses (described above), all interaction models analyses remained significant, suggesting that additive genetic effects account for effects rather than any single SNP (Supplementary Table 1).
Alternate environmental operationalizations
We conducted post hoc exploratory analyses using alternative environmental phenotypes. First, using the sum of the individual adversities as an index of CA, the G × E × E interaction terms remained significant for the HPA axis MGPS (using total episodic as the proximal variable, b = .002, SE = .001, p = .008) and serotonergic MGPS (using major interpersonal stress, b = .006, SE = .002, p = .001), with patterns of decomposition supporting increased stress sensitization at higher levels of both MGPSs. Second, using a more conservative threshold of “major” interpersonal events (3.0 versus 2.5), the G × E × E remained significant (p = .0001) in the serotonergic MGPS model. Finally, to replicate prior G × E × E work (Starr et al., Reference Starr, Hammen, Conway, Raposa and Brennan2014) which examined chronic (versus acute) stress, we re-ran models with total chronic stress as the proximal variable, using the original CA variable. The G × E × E effect was significant for the HPA axis MGPS (b = .11, SE = .04, p = .005), but not for the serotonergic MGPS (b = .03, SE = .04, p = .425). For completeness, we also tested these models using the summed severity variable for CA, and found significant G × E × E effects for both HPA axis MGPS (b = .007, SE = .003, p = .008) and serotonergic MGPS (b = .02, SE = .01, p = .032).
Replication sample method
Participants and overview
We sought to replicate the significant G × E × E effects for the serotonergic MGPS in an independent sample.Footnote 2 Adolescent girls (N = 132) and their primary female caregivers (herein called mothers) were recruited using flyers, word-of-mouth, and local schools. At baseline (T1), mothers and daughters each completed separate diagnostic and objective stress interviews, and adolescents provided a DNA sample using an Oragene saliva collection kit. Three girls were excluded who did not provide DNA. One-year later, 85.3% (n = 110) participated in a follow-up (T2) that included the same interviews. Participants (T1 M = 12.36, SD = .76; 86% White) who provided DNA were included in analyses, regardless of their participation at T2 (missing data were estimated; see below). Participants were mostly from middle- to upper-class families (family income ≥ $61,000: n = 79 [61.24%]). The Williams College Institutional Review Board approved study procedures. Using parameters from the Primary Sample's serotonergic MGPS results (R 2 = .09), power was estimated at .94, suggesting this sample is adequately powered to detect a significant G × E × E effect.
Measures and procedure
DNA collection and genotyping
DNA extractions were performed using the Oragene PrepIT L2P DNA Purification Kit and genotyping was carried out using allele-specific PCR with primers designed using PrimerPicker Software (LGC Genomics, Hertfordshire, UK). Genotypes were ascertained via fluorescent detection with a Synergy II Plate reader (Biotek Instruments, Winooski, VT). Full genotyping details are available on request. The serotonergic MGPS was calculated as in the primary sample.
Childhood adversity
CA was assessed at T1 in the same way as in the primary sample, except that both mothers and daughters each separately completed the YLSI given the early adolescent age of the sample, consistent with prior work (Rudolph & Flynn, Reference Rudolph and Flynn2007). As in prior reports (e.g., Rudolph & Flynn, Reference Rudolph and Flynn2007), information provided by mothers and daughters was combined into a single narrative. If mothers and daughters endorsed the same adversity, the narratives reflected both of their reports. If only the mother or only the daughter endorsed the adversity, the narrative was based upon only one person's report. A second team, blind to original ratings, recoded narratives for reliability (total severity: ICC = .99).
Proximal episodic stress
At T2, proximal episodic stress occurring between T1 and T2 was assessed using the same assessment and coding methods as the primary sample except that adolescents and mothers each completed separate interviews to assess adolescent's exposure to episodic stress, consistent with prior work in early adolescent samples (Rudolph et al., Reference Rudolph, Hammen, Burge, Lindberg, Herzberg and Daley2000). When mothers and daughters reported the same event, information from mothers and adolescents was combined into a single narrative. If only one reported the event, the narrative reflected only her report (e.g., Rudolph et al., Reference Rudolph, Hammen, Burge, Lindberg, Herzberg and Daley2000). As in the primary sample, the sum of the severity ratings for major interpersonal events was used. A second team, blind to the original ratings, rerated a set of events (n = 132; ICCs = .92 (objective impact); .98 (interpersonal status)).
Depressive symptoms
T2 depressive symptoms were assessed using the K-SADS-PL (Kaufman et al., Reference Kaufman, Birmaher, Brent, Rau, Flynn, Moreci and Ryan1997) current major depression module. A similar 4-point dimensional scale was employed as in the primary sample, with slightly different increments (0 = no symptoms; 1 = mild symptoms; 2 = moderate, subthreshold symptoms; 3 = meets DSM-IV criteria). Interrater reliability: ICC = .95. Consistent with other early adolescent community samples (Davila et al., Reference Davila, Stroud, Starr, Miller, Yoneda and Hershenberg2009; Rohde, Beevers, Stice, & O'Neil, Reference Rohde, Beevers, Stice and O'Neil2009), depression rates were relatively low, with 9% reporting at least subthreshold symptoms and 1% reporting diagnosable depression.
Pubertal status (covariate)
At T1, pubertal status was assessed via adolescent self-report using the PDS (Petersen et al., Reference Petersen, Crockett, Richards and Boxer1988).
Replication sample results
Preliminary analyses and racial effects
Bivariate correlations are presented in Table 1. Tests of racial effects revealed some racial differences related to the serotonergic MGPS and the environmental variables. White adolescents showed significantly higher MGPS than non-White adolescents, t (127) = −2.764, p = .007 and lower levels of CA, r = −.189, p = .043, though there was no evidence that White ancestry moderated the relation between the serotonergic MGPS and T2 depressive symptoms. Thus, to ensure that population stratification was not driving the results, we restricted analyses to the subsample of White adolescents (n = 111), but repeated analyses in the full sample (N = 129).
Serotonergic MGPS × major interpersonal stress × CA
Path analysis was conducted in Mplus 8 (Muthen & Muthen, Reference Muthen and Muthen1998–2018) because of missing data due to attrition. Analyses were conducted with maximum likelihood estimation with robust standard errors (MLR). Little's MCAR test indicated that data were missing completely at random (χ2 (2) = 1.66; p = .44). Full information maximum likelihood was used to handle missing data (Savalei & Rhemtulla, Reference Savalei and Rhemtulla2012). The model included serotonergic MGPS, major interpersonal stress, and CA, and their two- and three-way interaction effects predicting T2 depressive symptoms. Predictor variables were standardized. Interaction variables were computed by multiplying the standardized predictor variables with each other. Covariances were included between the disturbances of: (a) predictor variables; and (b) each of the predictor variables and the interaction variables.
Consistent with the results of the primary sample, the three-way interaction between serotonergic MGPS, major interpersonal stress, and CA was significant (b = .081 [.005, .158], SE = .039, p = .037; see Table 3). To decompose this interaction, we conducted a series of nine simple slope tests at high (M + SD), mean, and low (M-SD) levels of MGPS and high, mean, and low levels of CA. Consistent with the hypothesis that serotonergic MGPS increases stress sensitization, the association between major interpersonal stress and T2 depressive symptoms was significant only among those with high levels of both serotonergic MGPS and CA (b = .185 [.060, .311], SE = .064, p = .004). None of the other simple slopes were significant (all ps > .253). In a robustness test, T1 pubertal status and T1 age were added to the model, and a covariance was added between the covariates. The three-way interaction was near significant (p = .051). Finally, we tested the model in the full sample (N = 129), and the three-way interaction remained significant (p = .033), even when adding T1 age and T1 pubertal status as covariates (p = .044). Full results of the follow-up models and simple slope tests are available upon request.
Discussion
The current study showed that both HPA axis and serotonergic multilocus genetic variation predicted increased stress sensitization—depressive reactivity to recent events as a function of CA—among adolescents. These findings extend previous single-variant G × E × E findings (e.g., Starr et al., Reference Starr, Hammen, Conway, Raposa and Brennan2014) by using a more statistically powerful MGPS approach. An independent sample fully replicated the serotonergic G × E × E results (genotyping data were unavailable to replicate the HPA axis MGPS findings), boosting confidence in the results.
Results illustrate the interactive effects of stress exposures on depression during different developmental periods (early adolescence, mid-adolescence) and show that these effects are in turn potentiated by genetic variation, capturing the dynamic complexity of depression risk processes at multiple levels of analysis. Further, findings show that previously demonstrated G × E models, where both HPA axis and serotonergic multilocus genetic risk moderate the association between recent stress exposure and depression (Feurer et al., Reference Feurer, McGeary, Knopik, Brick, Palmer and Gibb2017; Starr & Huang, Reference Starr and Huang2019; Starr et al., Reference Starr, Dienes, Li and Shaw2019a; Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Stroud, Mineka, Zinbarg, Adam, Redei and Craske2015), appear to be further modified by CA exposure, with stronger G × E effects for adolescents who have experienced higher levels of adversity earlier in childhood. This could in part account for heterogeneity in G × E findings across studies, as studies that have recruited from populations with relatively high levels of CA exposure (e.g., selected for high maltreatment rates) would be expected to show more robust two-way G × E effects, whereas more normative or high-functioning samples may yield weaker G × Es.
An important next step is to identify specific mechanisms that place adolescents high on serotonergic and HPA axis MGPS at higher risk for stress sensitization. For example, both MGPSs are likely to influence HPA axis activity, which has long been argued as an important factor in stress sensitization (e.g., Heim et al., Reference Heim, Newport, Mletzko, Miller and Nemeroff2008). The SNPs on the HPA axis MGPS were largely selected for their association with cortisol dysregulation, and initial research has linked the MGPS with cortisol reactivity (Pagliaccio et al., Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden and Barch2014). Serotonin also has downstream effects on the HPA axis, and serotonergic genetic risk has been meta-analytically linked to HPA axis functioning (Miller, Wankerl, Stalder, Kirschbaum, & Alexander, Reference Miller, Wankerl, Stalder, Kirschbaum and Alexander2013). Thus, it is very possible that both MGPSs predict aberrant HPA axis responses to CA exposures, such as a lower threshold for activation or greater cortisol output. In turn, increased HPA axis activity during sensitive windows may have important implications for developing neural circuits and neuroendocrine systems. These may include breakdown in the ability of glucocorticoid receptors to downregulate HPA activation in response to elevated cortisol levels (i.e., a negative feedback loop), through desensitization or altered density of this inhibitory receptor, and alterations in brain structures related to stress regulation (Heim & Binder, Reference Heim and Binder2012; Heim et al., Reference Heim, Newport, Mletzko, Miller and Nemeroff2008). Consequently, youth at high genetic risk may be less able to respond adaptively when proximal stressors occur later in development (e.g., Starr et al., Reference Starr, Dienes, Li and Shaw2019).
In addition to purely biological mechanisms, social factors should be considered. CA exposure disrupts secure attachment formation in early childhood and subsequent interpersonal development (e.g., Baer & Martinez, Reference Baer and Martinez2006; Huh, Kim, Yu, & Chae, Reference Huh, Kim, Yu and Chae2014; Johnson et al., Reference Johnson, Cohen, Gould, Kasen, Brown and Brook2002; Shaw & Vondra, Reference Shaw and Vondra1993), which may predict poor quality interpersonal support systems later in life (Hankin, Kassel, & Abela, Reference Hankin, Kassel and Abela2005). Those who are at high genetic risk may then be especially sensitive to this lack of interpersonal support (Kaufman et al., Reference Kaufman, Yang, Douglas-Palumberi, Houshyar, Lipschitz, Krystal and Gelernter2004), especially following recent stress exposure (Kilpatrick et al., Reference Kilpatrick, Koenen, Ruggiero, Acierno, Galea, Resnick and Gelernter2007). Evidence also suggests that those at high genetic risk may show more interpersonal problems following CA (Harkness et al., Reference Harkness, Bagby, Stewart, Larocque, Mazurka, Strauss and Kennedy2015; Huang & Starr, Reference Huang and Starr2019), which may further erode social support systems, leaving youth more vulnerable to stress reactivity.
Critically, the results indicate that the serotonergic MGPS's influence over stress sensitization on depression emerged only for major interpersonal recent stress, and not for total recent stress which included noninterpersonal stress. This extends prior evidence for a “candidate environment” for serotonergic genetic risk on depression, namely evidence that the serotonin transporter polymorphism interacted with interpersonal but not noninterpersonal stress (Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Doane, Mineka, Zinbarg, Craske and Adam2013; Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Mineka, Zinbarg, Craske, Griffith, Sutton and Adam2014) and that this serotonergic MGPS produced a significantly greater G × E interaction with interpersonal recent major stress than with noninterpersonal recent major stress (Starr et al., Reference Starr, Dienes, Li and Shaw2019a). Because it is unclear whether other neurobiological systems’ genetic variation might also interact with interpersonal forms of stress in particular, and because interpersonal stress is significantly more potent for depression (Hammen, Reference Hammen2005), future G × E and G × E × E research should continue to test interpersonal stress as a candidate environment for depression, potentially in addition to other stress variables.
Our results are consistent with several studies documenting single-variant stress sensitization G × E × E effects for a serotonin-linked polymorphism, 5-HTTLPR (Grabe et al., Reference Grabe, Schwahn, Mahler, Schulz, Spitzer, Fenske and Freyberger2012; Kumsta et al., Reference Kumsta, Stevens, Brookes, Schlotz, Castle, Beckett and Sonuga-Barke2010; Starr et al., Reference Starr, Hammen, Conway, Raposa and Brennan2014) and one study supporting an HPA axis-related genotype G × E × E (CRHR1 rs110402; Starr et al., Reference Starr, Hammen, Conway, Raposa and Brennan2014). Although these single-variant studies produced significant effects, sensitivity analyses suggest that our G × E × E effects are best accounted for by additive genetic variance, rather than individual SNPs, implying that polygenic approaches are more effective than single SNPs in capturing this kind of genetic risk. This is the first study testing a G × E × E model using neurobiological system-specific MGPSs (also see Keers & Pluess, Reference Keers and Pluess2017).
By capturing polygenic risk, the multilocus approach appears to offers superior power over single-variant approaches. Although increasing statistical power has universal advantages in research, it is perhaps especially important for G × E × E research. The candidate gene literature has been plagued by concerns about nonreliability of findings, prompting some to argue for sample size requirements that are near-impossible to implement without sacrificing the quality of environmental assessments, which, ironically, are related to the strength of G × E effects (Karg et al., Reference Karg, Burmeister, Shedden and Sen2011) Use of low-quality environmental assessment may be even more problematic in G × E × E research (given that two key environmental contexts are incorporated into the model). Further, three-way interactions generally require more power to detect, and use of MGPSs may help counter this inherent limitation of testing G × E × E models.
In both samples, CA was defined cumulatively, or experienced over the course of childhood until one year prior to baseline, rather than occurring during a discrete developmental window. Cumulative adversity (compared to specific events) contributes more strongly to allostatic load (Evans, Reference Evans2003; Evans, Kim, Ting, Tesher, & Shannis, Reference Evans, Kim, Ting, Tesher and Shannis2007; Lupien et al., Reference Lupien, Ouellet-Morin, Hupbach, Tu, Buss, Walker, McEwen, Cicchetti and Cohen2006), which may reflect HPA axis dysregulation, suggesting relevance to stress sensitization. Moreover, research suggests that exposure to early adversity has a cumulative (but nonadditive) effect on first onsets of psychological disorders during adolescence (McLaughlin et al., Reference McLaughlin, Conron, Koenen and Gilman2010). Nonetheless, this cumulative definition precludes us from examining the important issue of developmental timing. Numerous researchers have speculated that plasticity of neural system peaks during specific developmental windows, and that adversities occurring during these sensitive periods may be more likely to lead to stress sensitization effects, perhaps especially among genetically susceptible children (Bosch et al., Reference Bosch, Riese, Reijneveld, Bakker, Verhulst, Ormel and Oldehinkel2012; Heim & Binder, Reference Heim and Binder2012). However, little is known about when, precisely, these windows may occur; adversities tend to cluster together and rarely occur in isolation, making this question difficult to test even in the best designed studies. Despite these challenges, it is an imperative question that demands further attention, as pinpointing periods of maximal vulnerability would both enrich understanding of risk trajectories and allow for more precise prevention efforts.
Important caveats and study limitations should be noted. First, sample size was small by genetics standards and analyses may have been underpowered. Although our independent replication offers evidence for robustness, additional replication of both MGPS findings is needed, particularly in larger samples. Second, in both samples, CA encompassed relatively normative experiences (e.g., parental divorce, hospitalizations, grandparent deaths) as well as more severe (rarer) contexts, such as maltreatment or trauma (with normative experiences more heavily represented, given the community nature of the sample). Although this is not unlike other G × E × E research (Keers & Pluess, Reference Keers and Pluess2017; Starr et al., Reference Starr, Hammen, Conway, Raposa and Brennan2014), a large contingent of the stress sensitization literature has focused on childhood maltreatment, which may have different effects. Moreover, CA was quantified as a unidimensional construct, measured across all types of events, belying the crucial heterogeneity inherent in adverse experiences. Indeed, numerous studies have suggested that severity, chronicity, and the qualitative nature of CAs independently predict negative outcomes (Gabrielli, Jackson, Tunno, & Hambrick, Reference Gabrielli, Jackson, Tunno and Hambrick2017; Jackson, Gabrielli, Fleming, Tunno, & Makanui, Reference Jackson, Gabrielli, Fleming, Tunno and Makanui2014; Warmingham, Handley, Rogosch, Manly, & Cicchetti, Reference Warmingham, Handley, Rogosch, Manly and Cicchetti2019). Third, the primary study utilized cross-sectional data (with current depression and recent stressors assessed at the same time point) and both samples utilized retrospective assessments of CA, which may have been subject to biased recall. The primary sample (but not the replication sample) specifically relied on parental report of CA, which may have excluded events that occurred outside of the parent's awareness or that the parent was disinclined to report (e.g., their own perpetration of maltreatment). Fourth, although comparable to other community samples (Lewinsohn et al., Reference Lewinsohn, Hops, Roberts and Seeley1993), both samples were community recruited and rates of clinical depression were limited; results should be replicated in samples with more severe depression. Fifth, given genotyping data were unavailable to replicate the HPA axis MGPS results, replication is needed. Sixth, pubertal status (a covariate) was assessed via self-report rather than physician examination; although this practice very common, it is not ideal (Rasmussen et al., Reference Rasmussen, Wohlfahrt-Veje, Tefre de Renzy-Martin, Hagen, Tinggaard, Mouritsen and Main2015). Finally, studies have conducted analyses comparing MGPS scores to randomly generated profile composed of SNPs with similar profiles (see Di Iorio et al., Reference Di Iorio, Carey, Michalski, Corral-Frias, Conley, Hariri and Bogdan2017); this was not possible here because neither sample used genome-wide chips. Limitations were balanced by important study strengths, including the use of gold-standard assessments for both CA and proximal stress, replication in an independent sample with consistent methods, use of established MPGSs, and continuously defined, interview-assessed depression as the outcome variable.
Conclusion
We provide the first evidence that additive genetic variation in the serotonin system and the HPA axis respectively moderates stress sensitization in adolescents. For both biological systems, more “risky” multilocus genetic scores predicted greater associations between recent stressors and depression as a function of CA. These results highlight the utility of multilocus genetic profile scores and offer neurobiological insights into developmental and environmental etiological pathways to adolescent depression.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0954579420000474
Acknowledgements
The authors would like to thank participating families in both studies, as well as Fanny Mlawer, Y. Irina Li, and Meghan Huang, and the staffs of the Rochester Internalizing Disorder and Emotional Adjustment Lab and the Williams College Youth Emotion Center for their assistance with data collection, coding, and management, and the University of Wisconsin Biotechnology Center DNA Sequencing Facility (Joshua Hyman, Director) for providing genotyping services. Institutional funds from University of Rochester supported the primary study (L.R.S., Principal Investigator). Institutional funds from Williams College supported the replication study (C.B.S., Principal Investigator). The authors declare no conflicts of interest.
