Introduction

Mental illness results from the interplay between genetic susceptibility and environmental risk factors1,2. Previous studies have shown that the effects of environmental factors on traits may be partially heritable3 and moderated by genetics4,5. Major depressive disorder (MDD) is the most common psychiatric disorder with a lifetime prevalence of approximately 14% globally6 and with a heritability of approximately 37%7. There is strong evidence for the role of stressful life events (SLEs) as risk factor and trigger for depression8,9,10,11,12. Genetic control of sensitivity to stress may vary between individuals, resulting in individual differences in the depressogenic effects of SLE, i.e., genotype-by-environment interaction (GxE)4,13,14,15,16. Significant evidence of GxE has been reported for common respiratory diseases and some forms of cancer17,18,19,20,21,22, and GxE studies have identified genetic risk variants not found by genome-wide association studies (GWAS)23,24,25,26,27.

Interaction between polygenic risk of MDD and recent SLE are reported to increase liability to depressive symptoms4,16; validating the implementation of genome-wide approaches to study GxE in depression. Most GxE studies for MDD have been conducted on candidate genes, or using polygenic approaches to a wide range of environmental risk factors, with some contradictory findings28,29,30,31,32. Incorporating knowledge about recent SLE into GWAS may improve our ability to detect risk variants in depression otherwise missed in GWAS33. To date, three studies have performed genome-wide by environment interaction studies (GWEIS) of MDD and SLE34,35,36, but this is the first study to perform GWEIS of depressive symptoms using adult SLE in cohorts of relatively homogeneous European ancestry.

Interpretation of GxE effects may be hindered by gene–environment correlation. Gene–environment correlation denotes a genetic mediation of associations through genetic influences on exposure to, or reporting of, environments2,37. Genetic factors predisposing to MDD may contribute to exposure and/or reporting of SLE38. To tackle this limitation, measures of SLE can be broken down into SLE likely to be independent of a respondent’s own behaviour and symptoms, or into dependent SLE, in which participants may play an active role exposure to SLE39,40. Different genetic influences, including a higher heritability, are reported for dependent SLE compared to independent SLE38,41,42,43,44, suggesting that whereas GxE driven by independent SLE is likely to reflect a genetic moderation of associations between SLE and depression, GxE driven by dependent SLE may result from a genetic mediation of the association through genetically driven personality or behavioural traits. To test this, we analysed dependent and independent SLE scores separately in Generation Scotland (GS).

Stress contributes to many human conditions, with evidence of genetic vulnerability to the effect of SLE45. Therefore, genetic stress-response factors in MDD may also underlie the aetiology of other stress-linked disorders with which MDD is often comorbid46,47 (e.g., cardiovascular diseases48, diabetes49, chronic pain50 and inflammation51). We tested the hypothesis that pleiotropy and shared aetiology between mental and physical health conditions may be due in part to genetic variants underlying SLE effects in depression.

In this study, we conduct GWEIS of depressive symptoms incorporating data on SLE in two independent UK-based cohorts. We aimed to: (i) identify loci associated with depressive symptoms through genetic response to SLE; (ii) study dependent and independent SLE to support a contribution of genetically mediated exposure to stress; (iii) assess whether GxE effects improve the proportion of phenotypic variance in depressive symptoms explained by genetic additive main effects alone; and (iv) test for a significant overlap in the genetic aetiology of the response to SLE and mental and physical stress-related phenotypes.

Materials and methods

The core workflow of this study is summarised in Fig. 1.

Fig. 1: Study flowchart.
figure 1

Overview of the analyses conducted in this study: (i) identify loci associated with depressive symptoms through genetic response to SLE; (ii) test whether results of studying dependent and independent SLE support a contribution of genetically mediated exposure to stress; (iii) assess whether GxE effects improve the proportion of phenotypic variance in depressive symptoms explained by genetic additive main effects alone and (iv) test whether there is significant overlap in the genetic aetiology of the response to SLE and mental and physical stress-related phenotypes. Two core cohorts are used, Generation Scotland (GS) and UK Biobank (UKB). Summary statistics from genome-wide association studies (GWAS) and genome-wide by environment interaction studies (GWEIS) are used to generate polygenic risk scores (PRSs). Summary statistics from Psychiatric Genetic Consortium (PGC) Major Depressive Disorder (MDD) GWAS are also used to generate PRS (PRSMDD). PRS weighted by: additive effects (PRSD and PRSMDD), GxE effects (PRSGxE) and joint effects (the combined additive and GxE effect; PRSJoint), are used for phenotypic prediction. TSLE stands for total number of SLE reported. DSLE stands for SLE dependent on an individual’s own behaviour. Conversely, ISLE stands for independent SLE. N stands for sample size. NnoGS stands for sample size with GS individuals removed. NnoUKB stands for sample size with UKB individuals removed

Full size image

Cohort descriptions

GS

GS is a family-based population cohort representative of the Scottish population52. At baseline, blood and salivary DNA samples were collected, stored and genotyped at the Wellcome Trust Clinical Research Facility, Edinburgh. Genome-wide genotype data were generated using the Illumina HumanOmniExpressExome-8 v1.0 DNA Analysis BeadChip (San Diego, CA, USA) and Infinium chemistry53. The procedures and details for DNA extraction and genotyping have been extensively described elsewhere54,55. In total, 21,525 participants were re-contacted to participate in a follow-up mental health study (Stratifying Resilience and Depression Longitudinally, STRADL), of which 8541 participants responded providing updated measures in psychiatric symptoms and SLE through self-reported mental health questionnaires56. Samples were excluded if: they were duplicate samples, had diagnoses of bipolar disorder, no SLE data (non-respondents), were population outliers (mainly non-Caucasians and Italian ancestry subgroup), had sex mismatches or were missing >2% of genotypes. Single nucleotide polymorphisms (SNPs) were excluded if: missing >2% of genotypes, Hardy–Weinberg equilibrium test p < 1 × 10−6, or minor allele frequency <1%. Further details of the GS and STRADL cohort are available elsewhere52,56,57,58. All components of GS and STRADL obtained ethical approval from the Tayside Committee on Medical Research Ethics on behalf of the NHS (reference 05/s1401/89). After quality control, individuals were filtered by degree of relatedness (pi-hat < 0.05), maximising retention of those individuals reporting a higher number of SLE. The final dataset comprised data on 4919 unrelated individuals (1929 men; 2990 women) and 560,351 SNPs.

Independent GS datasets

Additional datasets for a range of stress-linked medical conditions and personality traits were created from GS (N = 21,525) excluding respondents and their relatives (N = 5724). Following the same quality control criteria detailed above, we maximised unrelated non-respondents for retention of cases, or proxy cases (see below), to maximise the information available for each phenotype. This resulted in independent datasets with unrelated individuals for each trait. Differences between respondents and non-respondents are noted in the figure legend of Table 1.

Table 1 GS samples with stress-related phenotypes
Full size table

UK Biobank (UKB)

This study used data from 99,057 unrelated individuals (47,558 men; 51,499 women) from the initial release of UKB genotyped data (released 2015; under UKB project 4844). Briefly, participants were removed based on UKB genomic analysis exclusion, non-white British ancestry, high missingness, genetic relatedness (kinship coefficient > 0.0442), QC failure in UK BiLEVE study and gender mismatch. GS participants and their relatives were excluded and GS SNPs imputed to a reference set combining the UK10K haplotype and 1000 Genomes Phase 3 reference panels59. After quality control, 1,009,208 SNPs remained. UKB received ethical approval from the NHS National Research Ethics Service North West (reference: 11/NW/0382). Further details on UKB cohort description, genotyping, imputation and quality control are available elsewhere60,61,62.

All participants provided informed consent.

Phenotype assessment

SLEs

GS participants reported SLE experienced over the preceding 6 months through a self-reported brief life events questionnaire based on the 12-item list of threatening experiences39,63,64 (Supplementary Table 1a). The total number of SLE reported (TSLE) consisted of the number of ‘yes’ responses. TSLE were subdivided into SLE potentially dependent or secondary to an individual’s own behaviour (DSLE, questions 6–11 in Supplementary Table 1a), and independent SLE (ISLE, questions 1–5 in Supplementary Table 1a; pregnancy item removed) following Brugha et al.39,40. Thus, three SLE measures (TSLE, DSLE and ISLE) were constructed for GS. UKB participants were screened for ‘illness, injury, bereavement and stress’’ (Supplementary Table 1b) over the previous 2 years using six items included in the UKB Touchscreen questionnaire. A score reflecting SLE reported in UKB (TSLEUKB) was constructed by summing the number of ‘yes’ responses.

Psychological assessment

GS participants reported whether their current mental state over the preceding 2 weeks differed from their typical state using a self-administered 28-item scaled version of the General Health Questionnaire (GHQ)65,66,67. Participants rated the degree and severity of their current symptoms with a four-point Likert scale (following Goldberg et al.67). A final log-transformed GHQ was used to detect altered psychopathology and thus, assess depressive symptoms as results of SLE. In UKB participants, current depressive symptoms over the preceding 2 weeks were evaluated using four psychometric screening items (Supplementary Table 2), including two validated and reliable questions for screening depression68, from the Patient Health Questionnaire (PHQ) validated to screen mental illness69,70. Each question was rated in a four-point Likert scale to assess impairment/severity of symptoms. Due to its skewed distribution, a four-point PHQ score was formed from PHQ (0 = 0; 1 = 1–2; 2 = 3–5; 3 = 6 or more) to create a more normal distribution.

Stress-related traits

Targeted GS stress-related phenotypes and sample sizes are shown in Table 1 and detailed elsewhere52. These conditions were selected from literature review based on previous evidence of a link with stress45 (see also Supplementary Material: third section). Furthermore, we created additional independent samples using mapping by proxy, where individuals with a self-reported first-degree relative with a selected phenotype were included as proxy cases. This approach provides greater power to detect susceptibility variants in traits with low prevalence71.

Statistical analyses

SNP-heritability and genetic correlation

A restricted maximum likelihood approach was applied to estimate SNP-heritability (h2SNP) of depressive symptoms and self-reported SLE measures, and within samples bivariate genetic correlation between depressive symptoms and SLE measures using GCTA72.

GWAS analyses

GWAS were conducted in PLINK73. In GS, age, sex and 20 principal components (PCs) were fitted as covariates. In UKB, age, sex and 15 PCs recommended by UKB were fitted as covariates. The genome-wide significance threshold was p = 5 × 10–8.

GWEIS analyses

GWEIS were conducted on GHQ (the dependent variable) for TSLE, DSLE and ISLE in GS and on PHQ for TSLEUKB in UKB fitting the same covariates detailed above to reduce error variance. GWEIS were conducted using an R plugin for PLINK73 developed by Almli et al.74 (https://epstein-software.github.io/robust-joint-interaction). This method implements a robust test that jointly considers SNP and SNP–environment interaction effects from a full model (Y ~ β0 + βSNP + βSLE + βSNPxSLE + βCovariates) against a null model where both the SNP and SNP×SLE effects equal 0, to assess the joint effect (the combined additive main and GxE genetic effect at a SNP) using a nonlinear statistical approach that applies Huber–White estimates of variance to correct possible inflation due to heteroscedasticity (unequal variances across exposure levels). This robust test should reduce confounding due to differences in variance induced by covariate interaction effects if present75. Additional code was added (courtesy of Prof. Michael Epstein;74 Supplementary Material) to generate beta-coefficients and the p-value of the GxE term alone. In UKB, correcting for 1,009,208 SNPs and one exposure, we established a Bonferroni-adjusted threshold for significance at p = 2.47 × 10–8 for both joint and GxE effects. In GS, correcting for 560,351 SNPs and three measures of SLE we established a genome-wide significance threshold of p = 2.97 × 10–8.

Post-GWAS/GWEIS analyses

GWAS and GWEIS summary statistics were analysed using FUMA76 including: gene-based tests, functional annotation, gene prioritisation and pathway enrichment (Supplementary Material).

Polygenic profiling and prediction

Polygenic risk scores (PRSs) weighting by GxE effects (PRSGxE) were generated using PRSice-277 (Supplementary Material) in GS using GxE effects from UKB-GWEIS. In UKB, PRSGxE were constructed using GxE effects from all three GS-GWEIS (TSLE, DSLE and ISLE as exposures) independently. PRS were also weighted in both samples using either UKB-GWAS or GS-GWAS statistics (PRSD), and summary statistics from Psychiatric Genetic Consortium (PGC) MDD-GWAS (released 2016; PRSMDD) that excluded GS and UKB individuals when required (NnoGS = 155,866; NnoUKB = 138,884). Furthermore, we calculated PRS weighted by the joint effects (the combined additive main and GxE genetic effects; PRSJoint) from either the UKB-GWEIS or GS-GWEIS. PRS predictions of depressive symptoms were permuted 10,000 times. Multiple regression models fitting PRSGxE and PRSMDD, and both PRSGxE and PRSD were tested. All models were adjusted by same covariates used in GWAS/GWEIS. Null models were estimated from the direct effects of covariates alone. The predictive improvement of combining PRSGxE and PRSMDD/PRSD effects over PRSMDD/PRSD effects alone was tested for significance using the likelihood ratio test (LRT).

Prediction of PRSD, PRSGxE and PRSJoint of stress-linked traits were adjusted by age, sex and 20 PCs; and permuted 10,000 times. Empirical-p-values after permutations were further adjusted by false discovery rate (FDR, conservative threshold at Empirical-p = 6.16 × 10–3). The predictive improvement of fitting PRSGxE combined with PRSD and covariates over prediction of a phenotype using the PRSD effect alone with covariates was assessed using LRT, and LRT-p-values adjusted by FDR (conservative threshold at LRT-p = 8.35 × 10–4).

Results

Phenotypic and genetic correlations

Depressive symptom scores and SLE measures were positively correlated in both UKB (r2 = 0.22, p < 2.2 × 10–16) and GS (TSLE-r2 = 0.21, p = 1.69 × 10−52; DSLE-r2 = 0.21, p = 8.59 × 10−51; ISLE-r2 = 0.17, p = 2.33 × 10−33). Significant bivariate genetic correlation between depression and SLE scores was identified in UKB (rG = 0.72; p < 1 × 10−5, N = 50,000), but not in GS (rG = 1, p = 0.056, N = 4919; Supplementary Table 3a).

SNP-heritability (h 2 SNP)

In UKB, a significant h2SNP of PHQ was identified (h2SNP = 0.090; p < 0.001; N = 99,057). This estimate remained significant after adjusting by TSLEUKB effect (h2SNP = 0.079; p < 0.001), suggesting a genetic contribution unique to depressive symptoms. The h2SNP of TSLEUKB was also significant (h2SNP = 0.040, p < 0.001; Supplementary Table 3b). In GS, h2SNP was not significant for GHQ (h2SNP = 0.071, p = 0.165; N = 4919). However, in an ad hoc estimation from the baseline sample of 6751 unrelated GS participants (details in Supplementary Table 3b) we detected a significant h2SNP for GHQ (h2SNP = 0.135; p < 5.15 × 10−3), suggesting that the power to estimate h2SNP in GS may be limited by sample size. Estimates were not significant for either TSLE (h2SNP = 0.061, p = 0.189; Supplementary Table 3b) or ISLE (h2SNP = 0.000, p = 0.5), but h2SNP was significant for DSLE (h2SNP = 0.131, p = 0.029), supporting a potential genetic mediation and gene–environment correlation.

GWAS of depressive symptoms

No genome-wide significant SNPs were detected by GWAS in either cohort. Top findings (p < 1 × 10−5) are summarised in Supplementary Table 4. Manhattan and QQ plots are shown in Supplementary Figures 1-4. There was no evidence of genomic inflation (all λ1000 < 1.01).

Post-GWAS analyses

Gene-based test identified six genes associated with PHQ using the UKB-GWAS statistics at genome-wide significance (Bonferroni-corrected p = 2.77 × 10−6; DCC, p = 7.53 × 10−8; ACSS3, p = 6.51 × 10−7; DRD2, p = 6.55 × 10−7; STAG1, p = 1.63 × 10−6; FOXP2, p = 2.09 × 10−6; KYNU, p = 2.24 × 10−6; Supplementary Figure 8). Prioritised genes based on position, expression quantitative trait loci (eQTL) and chromatin interaction mapping are detailed in Supplementary Table 5. No genes were detected in GS-GWAS gene-based test (Supplementary Figures 9). No tissue-specific enrichment was detected from GWAS in either cohort. Significant gene-sets and GWAS catalogue associations for UKB-GWAS are reported in Supplementary Table 6. These included the biological process: positive regulation of long-term synaptic potentiation, and GWAS catalogue associations: brain structure, schizophrenia, response to amphetamines, age-related cataracts (age at onset), fibrinogen, acne (severe), fibrinogen levels and educational attainment; all adjusted-p < 0.01. There was no significant gene-set enrichment from GS-GWAS.

GWEIS of depressive symptoms

Manhattan and QQ plots are shown in Supplementary Figures 1-4. There was no evidence of GWEIS inflation for either UKB or GS (all λ1000 < 1.01). No genome-wide significant GWEIS associations were detected for SLE in UKB. GS-GWEIS using TSLE identified a significant GxE effect (p < 2.97 × 10−8) at an intragenic SNP on chromosome 11 (rs12789145, p = 4.95 × 10−9, β = 0.06, closest gene: PIWIL4; Supplementary Figure 5), and using DSLE at an intronic SNP in ZCCHC2 on chromosome 18 (rs17070072, p = 1.46 × 10−8, β = −0.08; Supplementary Figure 6). In their corresponding joint effect tests, both rs12789145 (p = 2.77 × 10−8) and rs17070072 p = 1.96 × 10−8) were significant. GWEIS for joint effect using DSLE identified two further significant SNPs on chromosome 9 (rs12000047, p = 2.00 × 10−8, β = −0.23; rs12005200, p = 2.09 × 10−8, β = −0.23, LD r2 > 0.8, closest gene: CYLC2; Supplementary Figure 7). None of these associations replicated in UKB (p > 0.05), although the effect direction was consistent between cohorts for the SNP close to PIWL1 and SNPs at CYLC2. No SNP achieved genome-wide significant association in the GS-GWEIS using ISLE as exposure. Top GWEIS results (p < 1 × 10−5) are summarised in Supplementary Tables 7-10.

Post-GWEIS analyses: gene-based tests

All results are shown in Supplementary Figures 10-17. Two genes were associated with PHQ using the joint effect from the UKB-GWEIS (ACSS3 p = 1.61 × 10−6; PHF2, p = 2.28 × 10−6; Supplementary Figure 11). ACSS3 was previously identified using the additive main effects, whereas PHF2 was only significantly associated using the joint effects. Gene-based tests identified MTNR1B as significantly associated with GHQ on the GS-GWEIS using DSLE in both GxE (p = 1.53 × 10−6) and joint effects (p = 2.38 × 10−6; Supplementary Figures 14-15).

Post-GWEIS analyses: tissue enrichment

We prioritised genes based on position, eQTL and chromatin interaction mapping in brain tissues and regions. In UKB, prioritised genes using GxE effects were enriched for upregulated differentially expressed genes from adrenal gland (adjusted-p = 3.58 × 10−2). Using joint effects, prioritised genes were enriched on upregulated differentially expressed genes from artery tibial (adjusted-p = 4.34 × 10−2). In GS, prioritised genes were enriched: in upregulated differentially expressed genes from artery coronary (adjusted-p = 4.55 × 10−2) using GxE effects with DSLE; in downregulated differentially expressed genes from artery aorta tissue (adjusted-p = 4.71 × 10−2) using GxE effects with ISLE; in upregulated differentially expressed genes from artery coronary (adjusted-p = 5.97 × 10−3, adjusted-p = 9.57 × 10−3) and artery tibial (adjusted-p = 1.05 × 10−2, adjusted-p = 1.55 × 10−2) tissues using joint effects with both TSLE and DSLE; and in downregulated differentially expressed genes from lung tissue (adjusted-p = 3.98 × 10−2) and in up- and downregulated differentially expressed genes from the spleen (adjusted-p = 4.71 × 10−2) using joint effects with ISLE. There was no enrichment using GxE effect with TSLE.

Post-GWEIS analyses: gene-sets enrichment

Significant gene-sets and GWAS catalogue hits from GWEIS are detailed in Supplementary Tables 11-14, including for UKB Biocarta: GPCR pathway; Reactome: opioid signalling, neurotransmitter receptor binding and downstream transmission in the postsynaptic cell, transmission across chemical synapses, gastrin CREB signalling pathway via PKC and MAPK; GWAS catalogue: post bronchodilator FEV1/FVC ratio, migraine and body mass index. In GS, enrichment was seen using TSLE and DLSE for GWAS catalogue: age-related macular degeneration, myopia, urate levels and Heschl’s gyrus morphology; and using ISLE for biological process: regulation of transporter activity. All adjusted-p < 0.01.

Cross-cohort prediction

In GS, PRSD weighted by the UKB-GWAS of PHQ significantly explained 0.56% of GHQ variance (Empirical-p < 1.10−4), similar to PRSMDD weighted by PGC MDD-GWAS (R2 = 0.78%, Empirical-p < 1.10−4). PRSGxE weighted by the UKB-GWEIS GxE effects explained 0.15% of GHQ variance (Empirical-p = 0.03, Supplementary Table 15). PRSGxE fitted jointly with PRSMDD significantly improved prediction of GHQ (R2 = 0.93%, model p = 6.12 × 10−11; predictive improvement of 19%, LRT-p = 5.91 × 10−3) compared with PRSMDD alone. Similar to PRSGxE with PRSD (R2 = 0.69%, model p = 2.72 × 10−8; predictive improvement of 23%, LRT-p = 0.01). PRSJoint weighted by the UKB-GWEIS also predicted GHQ (R2 = 0.58%, Empirical-p < 1.10−4), although the variance explained was significantly reduced compared with the model fitting PRSGxE and PRSD together (LRT-p = 4.69 × 10−7), suggesting that additive and GxE effects should be modelled independently for polygenic approaches (Fig. 2a).

Fig. 2: Prediction of depression scores by PRSGxE, PRSD, PRSMDD and PRSJoint.
figure 2

Variance of depression score explained by PRSGxE PRSD, PRSMDD and PRSJoint as single effect; and combining both PRSD and PRSMDD with PRSGxE in single models. Prediction was conducted using a Generation Scotland (GS) and b UK Biobank (UKB) as target sample. PRSGxE were weighted by cross- sample genome-wide by environment interaction studies (GWEIS) using GxE effect. PRSD were weighted by cross-sample genome-wide association studies (GWAS) of depressive symptoms effect. PRSMDD was weighted by Psychiatric Genetic Consortium (PGC) Major Depressive Disorder (MDD)-GWAS summary statistics. PRSJoint were weighted by cross-sample GWEIS using joint effect. A nominally significant gain in variance explained of General Health Questionnaire (GHQ) of about 23% was seen in GS when PRSGxE was incorporated into a multiple regression model along with PRSD; and of about 19% when PRSGxE was incorporated into a multiple regression model along with PRSMDD. Such a gain was not seen in UKB, but it must be noted that both PRSD and PRSMDD also explains much less variance of PHQ in UKB than of GHQ in GS. Also note, a noticeably reduction of variance explained by PRSJoint compared with combined polygenic risk scores (PRS)/effects

Full size image

In UKB (Fig. 2b), both PRSD weighted by the GS-GWAS of GHQ and PRSMDD significantly explained 0.04 and 0.45% of PHQ variance, respectively (both Empirical-p < 1.10−4; Supplementary Table 15). PRSGxE derived from the GS-GWEIS GxE effect did not significantly predicted PHQ (TSLE-PRSGxE Empirical-p = 0.382; DSLE-PRSGxE Empirical-p = 0.642; ISLE-PRSGxE Empirical-p = 0.748). Predictive improvements using the PRSGxE effect fitted jointly with PRSMDD or PRSD were not significant (all LRT-p > 0.08). PRSJoint significantly predicted PHQ (TSLE-PRSJoint: R2 = 0.032%, Empirical-p < 1.10−4; DSLE-PRSJoint: R2 = 0.012%, Empirical-p = 4.3 × 10−3; ISLE-PRSJoint: R2 = 0.032%, Empirical-p < 1.10−4), although the variance explained was significantly reduced compared with the models fitting PRSGxE and PRSD together (all LRT-p < 1.48 × 10−3).

Prediction of stress-related traits

Prediction of stress-related traits in independent samples using PRSD, PRSGxE and PRSJoint are summarised in Fig. 3a and Supplementary Table 16. Significant trait prediction after FDR adjustment (Empirical-p < 6.16 × 10−3, FDR-adjusted Empirical-p < 0.05) using both UKB and GS PRSD was seen for: depression status, neuroticism and schizotypal personality. PRSGxE weighted by the GS-GWEIS GxE effect using ISLE significantly predicted depression status mapped by proxy (Empirical-p = 7.00 × 10−4, FDR-adjusted Empirical-p = 9.54 × 10−3).

Fig. 3: Polygenic risk score (PRS) prediction in independent Generation Scotland (GS) datasets.
figure 3

a Heatmap illustrating PRS prediction of a wide range of traits from GS listed in the x axis (Table 1). (R) refers to traits using mapping by proxy approach (i.e., where first-degree relatives of individuals with the disease are considered proxy cases and included into the group of cases). Y axis shows the discovery sample and the effect used to weight PRS. Numbers in cells indicate the % of variance explained, also represented by colour scale. Significance is represented by asterixes according to the following significance codes: **p < 0.01; *p < 0.05; in grey Empirical-p-values after permutation (10,000 times) and in yellow FDR-adjusted Empirical-p-values. b Predictive improvement by GxE effect in independent GS datasets. Heatmap illustrating the predictive improvement as a result of incorporating PRSGxE into a multiple model along with PRSD and covariates (full model), over a model fitting PRSD alone with covariates (null model); predicting a wide range of traits from GS listed in the x axis (Table 1). Covariates: age, sex and 20 PCs. (R) refers to traits using mapping by proxy approach (i.e., where first-degree relatives of individuals with the disease are consider proxy cases and included into the group of cases). PRSGxE are weighted by genome-wide by environment interaction studies (GWEIS) using GxE effects. PRSD were weighted by the genome-wide association studies (GWAS) of depressive symptoms additive main effects. The y axis shows the discovery sample used to weight PRS. Numbers in cells indicate the % of variance explained by the PRSGxE, also represented by colour scale. Notice that those correspond to the PRSGxE predictions in Fig. 3a when PRSGxE are weighted by GxE effects. Significance was tested by likelihood ratio tests (LRT): full model including PRSD + PRSGxE vs. null model with PRSD alone (covariates adjusted). Significance is represented by asterixes according to the following significance codes: ***p < 0.001; **p < 0.01; *p < 0.05; in grey LRT-p-values and in yellow FDR-adjusted LRT-p-values

Full size image

Nominally significant predictive improvements (LRT-p < 0.05) of fitting PRSGxE, over the PRSD effect alone, using summary statistics generated from both UKB and GS were detected for schizotypal personality, heart diseases and chronic obstructive pulmonary disease (COPD) by proxy (Fig. 3b). PRSGxE weighted by GS-GWEIS GxE effect using ISLE significantly improved prediction over PRSD effect alone of depression status mapped by proxy after FDR adjustment (LRT-p = 1.96 × 10−4, FDR-adjusted LRT-p = 2.35 × 10−2).

Discussion

This study performs GWAS and incorporates data on recent adult SLEs into GWEIS of depressive symptoms, identifies new loci and candidate genes for the modulation of genetic response to SLE; and provides insights to help disentangle the underlying aetiological mechanisms increasing the genetic liability through SLE to both depressive symptoms and stress-related traits.

SNP-heritability of depressive symptoms (h2SNP = 9–13%), were slightly higher than previous estimates from African-American populations34, and over a third larger than estimates in MDD from European samples78. h2SNP for PHQ in UKB (9.0%) remained significant after adjusting for SLE (7.9%). Thus, although some genetic contributions may be partially shared between depressive symptoms and reporting of SLE, there is still a relatively large genetic contribution unique to depressive symptoms. Significant h2SNP of DSLE in GS (13%) and TSLEUKB in UKB (4%), which is mainly composed of dependent SLE items, were detected similar to previous studies (8 and 29%)34,42. Conversely, there was no evidence for heritability of independent SLE. A significant bivariate genetic correlation between depressive symptoms and SLE (rG = 0.72) was detected in UKB after adjusting for covariates, suggesting that there are shared common variants underlying self-reported depressive symptoms and SLE. This bivariate genetic correlation was smaller than that estimated from African-American populations (rG = 0.97; p = 0.04; N = 7179)34. Genetic correlations between SLE measures and GHQ were not significant in GS (N = 4919; rG = 1; all p > 0.056), perhaps due to a lack of power in this smaller sample.

Post-GWAS gene-based tests detected six genes significantly associated with PHQ (DCC, ACSS3, DRD2, STAG1, FOXP2 and KYNU). Previous studies have implicated these genes in liability to depression (see Supplementary Table 17), and three of them are genome-wide significant in gene-based tests from the latest meta-analysis of major depression that includes UKB (DCC, p = 2.57 × 10−14; DRD2, p = 5.35 × 10−14; and KYNU, p = 2.38 × 10−6; N = 807,553)79. This supports the implementation of quantitative measures such as PHQ to detect genes underlying lifetime depression status80. For example, significant gene ontology analysis of the UKB-GWAS identified enrichment for positive regulation of long-term synaptic potentiation, and for previous GWAS findings of brain structure81, schizophrenia82 and response to amphetamines83.

The key element of this study was to conduct GWEIS of depressive symptoms and recent SLE. We identified two loci with significant GxE effect in GS. However, none of these associations replicated in UKB (p > 0.05). The strongest association was using TSLE at 53-kb downstream of PIWIL4 (rs12789145). PIWIL4 is brain expressed and involved in chromatin modification84, suggesting it may moderate the effects of stress on depression. It encodes HIWI2, a protein thought to regulate OTX2, that is critical for the development of forebrain and for coordinating critical periods of plasticity disrupting the integration of cortical circuits85,86. Indeed, an intronic SNP in PIWIL4 was identified as the strongest GxE signal in attention deficit hyperactivity disorder using mother’s warmth as environmental exposure87. The other significant GxE identified in our study was in ZCCHC2 using DSLE. This zinc-finger protein is expressed in blood CD4+ T cells and is downregulated in individuals with MDD88 and in those resistant to treatment with citalopram89. No GxE effect was seen using ISLE as exposure.

No significant locus or gene with GxE effect was detected in the UKB-GWEIS. Nevertheless, joint effects (the combined additive main and GxE genetic effects) identified two genes significantly associated with PHQ (ACSS3 and PHF2; see Supplementary Table 17). PHF2 was recently detected as genome-wide significant at the latest meta-analysis of depression79. Notably, PHF2 paralogs have previously been linked with MDD through stress-response in three other studies90,91,92. Joint effects analyses in GS also detected an additional significant association upstream CYLC2, a gene nominally associated (p < 1 × 10−5) with obsessive-compulsive disorder and Tourette’s syndrome93. Gene-based test from the GS-GWEIS identified a significant association with MTNR1B, a melatonin receptor gene, using DSLE (both GxE and joint effect; Supplementary Table 17). Genes prioritised using GxE effects were enriched in differentially expressed genes from several tissues including the adrenal gland, which releases cortisol into the blood-stream in response to stress, thus playing a key role in the stress-response system, reinforcing a potential role of GxE in stress-related conditions.

The different instruments and sample sizes available make it hard to compare results between cohorts. Whereas GS contains deeper phenotyping measurements of stress and depressive symptoms than UKB, the sample size is much smaller, which may be reflected in the statistical power required to reliably detect GxE effects. Furthermore, the presence and size of GxE effects are dependent on their parameterisation (i.e., the measurement, scale and distribution of the instruments used to test such interaction)94. Thus, GxE may be incomparable across GWEIS due to differences in both phenotype assessment and stressors tested. Although our results suggest that both depressive symptom measures are correlated with lifetime depression status, different influences on depressive symptoms from the SLE covered across studies may contribute to lack of stronger replication. Instruments in GS cover a wider range of SLE and are more likely to capture changes in depressive symptoms as consequence of their short-term effects. Conversely, UKB could capture more marked long-term effects, as SLE were captured over 2 years compared with the 6 months in GS. New mental health questionnaires covering a wide range of psychiatric symptoms and SLE in the latest release of UKB data provides the opportunity to create similar measures to GS in the near future. Further replication in independent studies with equivalent instruments is required to validate our GWEIS findings. Despite these limitations and a lack of overlap in the individual genes prioritised from the two GWEIS, replication was seen in the predictive improvement of using PRSGxE derived from the GWEIS GxE effects to predict stress-related phenotypes.

The third aim of this study was to test whether modelling GxE effects could improve predictive genetic models, and thus help to explain deviation from additive models and missing heritability for MDD95. Multiple regression models suggested that inclusion of PRSGxE weighted by GxE effects could improve prediction of an individual’s depressive symptoms over use of PRSMDD or PRSD weighted by additive effects alone. In GS, we detected a predictive gain of 19% over PRSMDD weighted by PGC MDD-GWAS, and a gain of 23% over PRSD weighted by UKB-GWAS (Fig. 2a). However, these findings did not surpass stringent Bonferroni correction and could not be validated in UKB. This may reflect in the poor predictive power of the PRS generated from the much smaller GS discovery sample. The results show a noticeably reduced prediction using PRSJoint weighted by joint effects, which suggests that the genetic architecture of stress-response is at least partially independent and differs from genetic additive main effects. Overall, our results from multiple regression models suggest that for polygenic approaches main and GxE effects should be modelled independently.

SLE effects are not limited to mental illness45. Our final aim was to investigate shared aetiology between GxE for depressive symptoms and stress-related traits. Despite the differences between the respondents and non-respondents (Table 1 legend), a significant improvement was seen in predicting depressive status when mapping by proxy cases using GxE effect from GS-GWEIS with independent SLE (FDR-adjusted LRT-p = 0.013), but not with dependent SLE. GxE effects using statistics generated from both discovery samples, despite the differences in measures, nominally improved the phenotypic prediction of schizotypal personality, heart disease and the proxy of COPD (LRT-p < 0.05). Other studies have also found evidence supporting a link between stress and depression in these phenotypes (see Supplementary Material for extended review) and suggest, for instance, potential pleiotropy between schizotypal personality and stress-response. Our findings point to a potential genetic component underlying a stress-response-depression-comorbidities link due, at least in part, to shared stress-response mechanisms. A relationship between SLE, depression and coping strategies such as smoking suggests that genetic stress-response may modulate adaptive behaviours such as smoking, fatty diet intake, alcohol consumption and substance abuse. This is discussed further in the Supplementary Material.

In this study, evidence for SNPs with significant GxE effects came primarily from the analyses of dependent SLE and not from independent SLE. This supports a genetic effect on probability of exposure to, or reporting of SLE, endorsing a gene–environment correlation. Chronic stress may influence cognition, decision making and behaviour eventually leading to higher risk taking96. These conditions may also increase sensitivity to stress among vulnerable individuals, including those with depression, who also have a higher propensity to report SLE, particularly dependent SLE38. A potential reporting bias in dependent SLE may be mediated as well by heritable behavioural, anxiety or psychological traits such as risk taking42,97. Furthermore, individuals vulnerable to MDD may behave in a manner that exposes them more frequently to high risk and stressful environments14. This complex interplay, reflected in the form of a gene–environment correlation effect, would hinder the interpretation of GxE effects from GWEIS as pure interactions. A mediation of associations between SLE and depressive symptoms, through genetically driven sensitivity to stress, personality or behavioural traits would support the possibility of subtle genotype-by-genotype (GxG) interactions, or genotype-by-genotype-by-environment (GxGxE) interactions, contributing to depression98,99. In contrast, PRS prediction of the stress-related traits: schizotypal personality, heart disease and COPD, was primarily from derived weights using independent SLE, suggesting that a common set of variants moderate the effects of SLE across stress-related traits and that larger sample sizes will be required to detect the individual SNPs contributing to this. Thus, our findings support the inclusion of environmental information into GWAS to enhance the detection of relevant genes. The results of studying dependent and independent SLE support a contribution of genetically mediated exposure to and/or reporting of SLE, perhaps through sensitivity to stress as mediator.

This study emphasises the relevance of GxE in depression and human health in general and provides the basis for future lines of research.