Polygenic risk scores for major depressive disorder and neuroticism as predictors of antidepressant response: Meta-analysis of three treatment cohorts.

There are currently no reliable approaches for correctly identifying which patients with major depressive disorder (MDD) will respond well to antidepressant therapy. However, recent genetic advances suggest that Polygenic Risk Scores (PRS) could allow MDD patients to be stratified for antidepressant response. We used PRS for MDD and PRS for neuroticism as putative predictors of antidepressant response within three treatment cohorts: The Genome-based Therapeutic Drugs for Depression (GENDEP) cohort, and 2 sub-cohorts from the Pharmacogenomics Research Network Antidepressant Medication Pharmacogenomics Study PRGN-AMPS (total patient number = 760). Results across cohorts were combined via meta-analysis within a random effects model. Overall, PRS for MDD and neuroticism did not significantly predict antidepressant response but there was a consistent direction of effect, whereby greater genetic loading for both MDD (best MDD result, p < 5*10-5 MDD-PRS at 4 weeks, β = -0.019, S.E = 0.008, p = 0.01) and neuroticism (best neuroticism result, p < 0.1 neuroticism-PRS at 8 weeks, β = -0.017, S.E = 0.008, p = 0.03) were associated with less favourable response. We conclude that the PRS approach may offer some promise for treatment stratification in MDD and should now be assessed within larger clinical cohorts.

Introduction

Major Depressive disorder (MDD) is a leading cause of disability worldwide [1]. Antidepressants such as Selective Serotonin Reuptake Inhibitors (SSRIs) are first line treatments for MDD but up to one third of patients do not respond satisfactorily [2, 3]. There are currently no robust methods for predicting whether an individual patient will respond well to SSRIs and there is often a lag period of several weeks before clinical response, making decisions on switching to a different class of antidepressant difficult. Individual genetic variation may dictate likelihood of response to SSRIs [4] and, as such, stratifying patients into sub-groups based on genetic profiles may allow for more efficient targeting of treatment.

Polygenic risk scoring (PRS) [5] is a method which allows an individual’s genetic loading for a trait to be calculated using genome-wide single nucleotide polymorphism (SNP) data and the output of genome-wide association study (GWAS) summary statistics from another study of the same or related phenotype. As current GWAS results do not capture the full extent of genetic effects on any given trait, typically a series of scores are created at different association p-value cut offs, allowing for the capture of more variance than that explained by only genome-wide significant loci. Additionally, as the underlying genetic architecture of the trait is unknown creating a range of scores can allow for the optimum p value threshold to be determined, should one detect a significant correlation.

It has been shown that a PRS can be of clinical use in predicting traits in independent samples. For example, for coronary heart disease, PRS improved the 10 year risk prediction in those over age 60 [6]. PRS approaches can also predict response to treatment, as demonstrated recently with an association between PRS for schizophrenia and less favourable response to lithium in bipolar disorder [7]. Here we test the hypothesis that PRS for MDD and PRS for neuroticism are associated with less favourable response to SSRIs, specifically citalopram and its active S-enantiomer escitalopram, in patients with MDD. Neuroticism is of particular interest in this regard because it has a known association with both serotonergic neurotransmission [8] and response to antidepressants [9, 10], and those with higher phenotypic neuroticism are less likely to respond as well to antidepressant therapy [11].

The analysis investigated three cohorts, GENDEP, AMPS-1 and AMPS-2 separately and then combine the results via meta-analysis.

Methods

Cohort descriptions, genotyping and imputation

The Pharmacogenomics Research Network Antidepressant Medication Pharmacogenomics Study (PGRN-AMPS) is a study of citalopram/escitalopram for treatment of MDD performed at the Mayo Clinic. An initial batch of 530 subjects (N = 499 subjects of European ancestry that passed quality control) was genotyped for a pharmacogenomics GWAS of SSRIs [12]. An additional 229 patients recruited in the PGRN-AMPS were subsequently genotyped for the International SSRI Pharmacogenomics Consortium (ISPC) GWAS [13]. Depressive symptoms were assessed on the Hamilton Depression Rating Scale (HAMD) with a maximum score of 51, a scale developed to rate both the psychiatric as well as the psychomotor and somatic symptoms of the condition[14]. Full genotyping and imputation of these cohorts (here referred to as AMPS-1 and AMPS-2) have been described previously [12, 13].

Genome Based Therapeutic Drugs for Depression (GENDEP) is a cohort of 868 individuals, recruited from across Europe, treated with two classes of antidepressants: escitalopram (an SSRI) and nortriptyline (a tricyclic antidepressant). For the purposes of this study, only those patients in GENDEP treated with an SSRI were assessed (n = 267). Depressive symptoms were assessed on the 10-item Montgomery-Asberg Depression Rating Scale (MADRS) with a maximum score of 60, with measurements taken weekly for 12 weeks from baseline. MADRS differs from HAMD in that it focuses exclusively on the psychiatric symptoms only and not the accompanying psychomotor and somatic symptoms of MDD [14]. Full genotyping and imputation methodology in GENDEP is described in previous reports [15].

Principal component generation and PRS construction

Principal genetic components were derived using PLINK. For all models the top 4 principal components were used as covariates in the model to account for hidden population structure. To ensure that an ethnically homogeneous sample was used in the AMPS-1 and AMPS-2 cohorts those whose Principal genetic components 1 to 4 were outside two standard deviations from the mean were excluded as outliers.

PRS were constructed via PLINK [16] with SNP weights based on outputs from the Smith et al. (2016) neuroticism GWAS [17] and the “probable MDD” phenotype of Howard et al (2018) MDD GWAS from UK Biobank[18]. SNPs were filtered by MAF < 0.01, HWE p<1*10−6 and imputation score < 0.8 before Linkage Disequilibrium (LD) clumping. SNPs were clumped using LD parameters of r2 >0.05 in a 500kb window. Selection of SNPs for each clump was based on which SNP had the lowest p value. If 2 SNPs in a clump had the same P value the SNP with the largest beta coefficient was selected. The scores generated were average scores with no-mean-imputation flag. Six profile scores were created for each trait using p value cut offs of p < 5*10−8, p < 5*10−5, p < 0.01, p < 0.05, p < 0.1 and p < 0.5. Risk scores were then standardised to mean = 0, SD = 1[19].

Due to low numbers and therefore the potential for noise within outcome data, instead of assessing change in outcomes across the full range of polygenic scores we chose to investigate only the difference between the extreme ends of the PRS scale. To do this, we split the standardised scores into quintiles and looked at the difference between the top and bottom quintile of each PRS p-value cut off within each cohort. For the GENDEP cohort the top and bottom quintile from each centre was selected to account for variation between recruitment centres. It is also important to note that an individual may be in the top quintile for one PRS P-value cut-off but not in another. As such, the two fifths of individuals used in each regression will change depending on the PRS p value cut off used.

Phenotype definition

For all three cohorts the primary outcome of interest was percentage change in depression score from baseline at four weeks. This was calculated by subtracting the score at four weeks from baseline, and dividing this difference by the score at baseline. A secondary outcome at eight weeks was also assessed, calculated using the same method. To be included in the analysis, an individual had to have a score recorded at baseline, four weeks and eight weeks.

Statistical modelling

Modelling was performed in R using the lm function. All models were adjusted for age, sex and the first 4 principal components. The GENDEP models were additionally adjusted for recruitment centre which was treated as a factor variable. The R2 for the PRS term of the model was derived using the methodology described in Selzam et al[20]. Due to the results being largely null we did not perform any correction for multiple testing.

Meta-analysis

A random effects Meta-analysis was performed using the rma.uni function of the metaphor package with method set to “REML”[21].

Results

Demographic and clinical characteristics of the three cohorts (GENDEP, AMPS-1, and AMPS-2) are presented in Table 1. The percentage female and age range of the three cohorts were broadly similar. The scores at baseline, 4 week and 8 week time points in AMPS-1 and AMPS-2 show a similar trend with a similar percentage drop at 4 and 8 week time points. The baseline scores of the GENDEP cohort are higher than in AMPS-1 and AMPS-2 due to the cohort being scored using MADRS and not HAMD as is the case with AMP-1 and AMPS-2. At Both the 4 week and 8 week time point the GENDEP cohort showed a smaller percentage reduction than in the AMPS-1 and AMPS-2. This difference may be explained by the differing depression measures picking up on differing aspects of MDD, differing healthcare settings and levels of severity at baseline. The within cohort drop from baseline at both 4 and 8 weeks was statistically significant for all three cohorts.

Table 1

Demographic and clinical characteristics.

Cohort	Total N	N used per regression	N female of total N (%)	Age of total N, mean(SD)	Baseline score*, mean (SD)	4 week score*, mean (SD)	8 week score*, mean (SD)	% drop in mean score at 4 weeks from mean score at baseline	% drop in mean score at 8 weeks from mean score at baseline
AMPS-1	357	142	229 (64.1)	40.9 (13.5)	22 (4.88)	11.9 (6.7)	8.83 (5.92)	46	60
AMPS-2	138	55	85 (61.6)	40.1 (13.6)	21.2 (5.14)	12 (5.84)	9.14 (6.41)	43	57
GENDEP	265	106	170 (64.2)	42.3 (11.8)	28.3 (6.16)	18.7 (8.2)	14.2 (8.89)	34	50

*score rating is HAMD for AMPS-1 and AMPS-2 and MADRS for GENDEP.

For neuroticism PRS in GENDEP, AMPS-1 and AMPS-2 the number of SNPs in each risk score were similar between cohorts across all p-value cut-offs (S1 Table). For the MDD risk scores the number of SNPs were similar between cohorts in the lower p value thresholds but diverged at the higher p value cut-offs. These differences arise mainly due to the differences in imputation coverage and the differing ethnicities and their impact on LD block estimation.

Individual study analyses

The results of all the individual study analyses can be found in S2–S4 Tables. Two of the models returned nominally significant results, both of which were in the AMPS-2 cohort (Table 2). They were neuroticism p < 0.5 PRS at four weeks (β = -0.04, p = 0.02) and neuroticism p < 0.5 at eight weeks (β = -0.039, p = 0.03). Of particular note is the R2 of the PRS term of the significant models which accounts for approximately 10% of the variance. Note, however, that these results would not pass correction for multiple testing.

Table 2

Nominally significant individual PRS models (AMSP-2 cohort).

predictor	Time point(weeks)	p	Beta	SE	T Test stat	R2
Neuroticism p<0.5	4	0.019	-0.044	0.018	-2.42	0.1
Neuroticism p<0.5	8	0.029	-0.039	0.017	-2.26	0.08

Although we were unable to reject the null hypothesis in the rest of the models, a clear majority (56 of 72 models) identified beta coefficients in the same direction of effect (greater loading for MDD or neuroticism associated with a smaller percentage drop in depression score). Of the 16 positive beta coefficient models, ten were from GENDEP MDD PRS models, three were from GENDEP neuroticism PRS model, two were from AMPS-1 neuroticism PRS models and one was an AMPS-2 MDD PRS models (S2–S4 Tables).

Two of the 24 meta-analyses were nominally significant: MDD p < 5*10−5 PRS at four weeks (β = -0.02, p = 0.009, I2 = 0); and neuroticism p<0.1 PRS at eight weeks (β = -0.017, p = 0.03, I2 = 0) (Fig 1). Neither of these results would survive correction for multiple testing. The direction of effect in all of the meta-analyses was negative (greater genetic loading for MDD and neuroticism associated with a smaller percentage drop in depression score at both four and eight weeks; S5 Table. The forest plots of all other meta-analyses are provided as supplementary material (S1–S4 Figs).

Fig 1

Forest plot of nominally significant meta-analyses.

A) p < 5*10−5 MDD-PRS at 4 weeks, B) p<0.1 Neuroticism-PRS at 8 weeks.

Discussion

Our goal was to assess the extent to which PRS for MDD and PRS for neuroticism were associated with response to SSRIs in patients with MDD. Although most of the findings were null, there was a direction of effect where higher PRS for MDD and higher PRS for neuroticism were associated with less favourable response to SSRIs. It is likely that our analyses were under-powered–replication in larger datasets will therefore be of interest. We estimate that a training sample of approximately 10,000 and a target sample of 5,000 individuals would give 60% power in a PRS of 100,000 SNPs that explain 10% of the variance in the training sample [22]. For the two AMPS-2 nominally significant results the R2 values of approximately 10%, suggesting that these PRSs could potentially be useful clinically.

This work diverges from previous analyses in these cohorts which have focused on GWAS and candidate gene analyses to identify genetic loci that associate with antidepressant response with the exception of Garcia-Gonzalez et al[23]. However, the outcome is markedly different to the outcome used here. It is possible that the use of PRS is advantageous for clinical use over these methods as it allows for a whole-genome approach instead of focusing on specific SNPs, genes or regions. An individual’s response to antidepressants is likely to be influenced by many genetic factors and, as such, candidate gene methodologies will fail to capture polygenic influences. An additional strength of this work is that all three cohorts systematically assessed treatment response at comparable time-points and in the context of the use of the same class of antidepressants, namely SSRIs.

Limitations

Apart from the issue of low power, our methodology was one in which only the extreme ends of genetic loadings were considered. This makes it difficult to translate the findings into a general population setting and routine clinical practice. Further work is needed to assess genetic loadings for MDD and neuroticism within the general population and how these relate to the clinical cohorts described here. The use of different depression rating scales between GENDEP and the AMPS-1/AMPS-2 may have had some impact on the results as they may have captured different aspects of the depressive phenotype and symptom changes induced by antidepressants. However, I2 was low in the meta-analyses that achieved nominal significance. Using a consistent depression rating in future would aid in keeping heterogeneity consistently low.

Another limitation was in the estimation of LD blocks in the GENDEP cohort. Due to the cohort being composed of individuals across Europe, treating the group as a whole for estimating which SNPs are in LD may have led to inaccuracies. This could explain why many of coefficients in the GENDEP models showed as positive correlation unlike the models from AMPS-1 and AMPS-2. Principal component analysis of treatment centres showed overlapping clusters but they were not distinct enough to warrant calculating LD in each centre separately. Further work in this area should capture more detail on ethnicity and ancestral background, to allow for more robust determination of LD clumps and more informed decisions on the most appropriate inclusion criteria.

Finally, the result may have been impeded by the use of a single PRS predictor. Recent research has shown that the use of multiple scores covering a variety genetic loadings can explain significantly more variance that that of a single score [24]. As such, incorporation of multiple genetic risk scores for outcomes as complex as antidepressant response may prove more fruitful.

Conclusion

Stratified medicine in psychiatry is still in its infancy. Genotyping is not currently routine practice in clinical settings and the use of PRS to guide the use of SSRIs in MDD remains a long-term goal.

However, with increasingly large and well-phenotyped cohorts available for analysis and more powerful GWAS outputs being produced, we tentatively conclude that more targeted prescribing of anti-depressants in MDD based on genetic profiles is a realistic prospect for the future.

Number of SNPs in each PRS for MDD and Neuroticism.

(XLSX)

Click here for additional data file.

Results of the individual regressions in the AMPS-1 cohort.

(XLSX)

Click here for additional data file.

Results of the induvial regressions in the AMPS-2 cohort.

(XLSX)

Click here for additional data file.

Results of the induvial regressions in the GENDEP cohort.

(XLSX)

Click here for additional data file.

Results of the meta-analyses.

(XLSX)

Click here for additional data file.

MDD PRS meta-analysis results at 4 weeks. A) p < 5*10−8, B) p < 5*10−5, C) p <0.01, D) p < 0.05, E) p < 0.1, F) p < 0.5.

(TIF)

Click here for additional data file.

MDD PRS meta-analysis results at 8 weeks A) p < 5*10−8, B) p < 5*10−5, C) p <0.01, D) p < 0.05, E) p < 0.1, F) p < 0.5.

(TIF)

Click here for additional data file.

Neuroticism PRS meta-analysis results at 4 weeks A) p < 5*10−8, B) p < 5*10−5, C) p <0.01, D) p < 0.05, E) p < 0.1, F) p < 0.5.

(TIF)

Click here for additional data file.

Neuroticism PRS meta-analysis results at 8 weeks A) p < 5*10−8, B) p < 5*10−5, C) p <0.01, D) p < 0.05, E) p < 0.1, F) p < 0.5.0.

(TIF)

Click here for additional data file.