Abstract To elucidate the genetic underpinnings of the antidepressant efficacy of S-ketamine (esketamine) nasal spray in major depressive disorder (MDD), we performed a genome-wide association study (GWAS) in cohorts of European ancestry (n = 527). This analysis was followed by a polygenic risk score approach to test for associations between genetic loading for psychiatric conditions, symptom profiles and esketamine efficacy. We identified a genome-wide significant locus in IRAK3 (p = 3.57 × 10^–8, rs11465988, β = − 51.6, SE = 9.2) and a genome-wide significant gene-level association in NME7 (p = 1.73 × 10^–6) for esketamine efficacy (i.e. percentage change in symptom severity score compared to baseline). Additionally, the strongest association with esketamine efficacy identified in the polygenic score analysis was from the genetic loading for depressive symptoms (p = 0.001, standardized coefficient β = − 3.1, SE = 0.9), which did not reach study-wide significance. Pathways relevant to neuronal and synaptic function, immune signaling, and glucocorticoid receptor/stress response showed enrichment among the suggestive GWAS signals. Subject terms: Genome-wide association studies, Depression Introduction Esketamine nasal spray has been shown to have rapidly-acting antidepressant effects in patients with treatment resistant depression (TRD) and in patients with major depressive disorder (MDD) at imminent risk for suicide^[32]1–[33]8. Predictors for conventional oral antidepressant treatment outcome including sociodemographic, symptom profiles, genetics, and clinical comorbidities were systematically reviewed by Perlman et al.^[34]9 In a small clinical study assessing the antidepressant efficacy of ketamine, a racemate consisting of two enantiomers, R- and S-ketamine, it was recently reported that body mass index (BMI) was associated with the remission rate, with greater BMI being associated with greater remission rate^[35]10. BMI and clinical comorbidities are influenced by both genetic and environmental factors. Genetic loading of such traits provides an objective way of measuring the relationship between these and other predictors with antidepressant treatment response. In studies assessing individual genetic factors, the brain-derived neurotrophic factor (BDNF) Val66Met allele was reported to impair basal and ketamine-stimulated synaptogenesis in prefrontal cortex in vitro^[36]11, and a significant genetic association between Val66Met and ketamine treatment outcome at 4 h post treatment was reported in a candidate gene study of small sample size^[37]12. A more recent study further suggested that the BNDF Val66Met polymorphism may influence the improvement in suicide ideation following ketamine infusion in a sample of depressed participants from Taiwan^[38]13. In general, however, genetic associations with MDD disease susceptibility outcome reported in relatively small candidate gene studies have proven difficult to replicate in studies of larger samples^[39]14. Therefore, studies of genetic effects influencing antidepressant treatment outcome may particularly benefit from the use of genome-wide association analysis (GWAS) approaches in clinical trials of larger patient samples. Here, we assessed the genetic contributions to esketamine treatment response from patients with TRD who participated in two Phase III trials testing the efficacy and safety of esketamine, using both a genome-wide association analysis and a polygenic risk score (PRS) approach. Results Esketamine treatment response outcome was assessed at the 4 week study endpoint using one continuous variable (percent change from baseline in the Montgomery–Asberg Depression Rating Scale (MADRS) score ) and two dichotomized variables (responder status, defined by a reduction of ≥ 50% on the MADRS, and remission status, defined by achieving a final MADRS score of < 12). The demographic and clinical characteristics of study participants are summarized in Table [40]1 and Supplemental Table [41]1. Participants of the randomized TRANSFORM-3 study were recruited from an elderly population and had a lower remission rate than participants of the open-labelled SUSTAIN-2 study. Gender and concomitant medication proportions were comparable between remitters and non-remitters. After controlling for study, the baseline demographic characteristics (age and baseline BMI) were comparable between remitters and non-remitters. As expected from the clinical literature, remitters had lower baseline depression symptom severity score than non-remitters. Table 1. Characteristics of study participants comparing remitters from non-remitters. Remitters (n = 255) Non-remitters (n = 272) p-value Mean (SD) Age* 50.6 (13.8) 53.4 (13.5) 0.424 Baseline BMI* 28.1 (5.6) 28.3 (5.8) 0.742 Baseline MADRS score* 29.7 (4.7) 33.0 (4.7) 6.36E-13 N (%) Gender, female 153 (60.0) 175 (64.3) 0.349 Study 7.34E−05 TRANSFORM-3 10 (3.9) 39 (14.3) SUSTAIN-2 245 (96.1) 233 (85.7) Concomitant antidepressant medications 0.782 DULOXETINE 90 (35.3) 87 (32.0) ESCITALOPRAM 78 (30.6) 80 (29.4) SERTRALINE 42 (16.5) 52 (19.1) VENLAFAXINE XR 45 (17.6) 52 (19.1) None 1 (0.4) [42]Open in a new tab *p-value reported is based on type III test statistics controlling for study. The genome-wide association analysis revealed one genome-wide significant association between an exonic synonymous variant (rs11465988, p = 3.57 × 10^–8) in the interleukin 1 receptor associated kinase 3 (IRAK3) gene and the percent change in MADRS score (Table [43]2, Fig. [44]1 for Manhattan plot, Fig. [45]2A for regional plot and Supplemental Fig. [46]1A for QQ plot, Genomic Control lambda (λ) = 0.986). SNPs (e.g. rs115989442, rs144324167, rs79138866, rs116371327, rs150373274, and rs144520864) in linkage disequilibrium (r^2 = 0.64) with rs11465988 are part of the regions engaging in intra-chromosomal loop (Fig. [47]2B, Supplemental Table [48]2) and could potentially be regulatory elements. rs144520864 is in fact located in a region with an annotated enhancer. An additional regional plot using rs17767394 as index SNP is also shown as Supplemental Fig. [49]2. Table 2. Variants with association p-value less than 1 × 10^–6 in GWAS. rsID Chr pos A1 A2 FRQ INFO Beta/OR SE p Func.refGene Gene.refGene GeneDetail.refGene Percentage change of MADRS from baseline rs11465988 12 66641813 C T 0.9898 0.51 − 51.6 9.2 3.57E−08 Exonic IRAK3 rs17767394 12 66636086 C A 0.9843 0.78 − 32.7 6 8.68E−08 Intronic IRAK3 rs4739050 8 64034747 G A 0.3376 1.02 7.5 1.4 6.06E−08 Intergenic TTPA;YTHDF3-AS1 dist = 36135; dist = 45537 rs151184257 4 105714757 A G 0.9888 0.61 − 40.5 7.9 4.51E−07 Intergenic CXXC4-AS1;TET2 dist = 96008; dist = 352275 rs115141868 2 70816605 A C 0.9898 0.48 − 46.8 9.3 7.65E−07 Intergenic TGFA;ADD2 dist = 35458; dist = 72611 Response status rs10957273 8 6.4E+07 T C 0.3028 1 0.3 0.2 8.07E−07 Intergenic TTPA;YTHDF3-AS1 dist = 30479; dist = 51193 [50]Open in a new tab Note that beta coefficient is reported for percentage change of MADRS score from baseline and OR is reported for responder status. Figure 1. [51]Figure 1 [52]Open in a new tab Manhattan plot of the esketamine efficacy endpoint (percentage change of MADRS score at endpoint compared to baseline) generated via FUMA^[53]51 v1.3.5e ([54]https://fuma.ctglab.nl/). The red dotted line indicates genome-wide significance threshold of 5 × 10^–8. Figure 2. [55]Figure 2 [56]Open in a new tab Genome-wide significant locus IRAK3. (A) Regional association plot; (B) circos plot. For the regional association plot generated via LocusZoom^[57]52 v1.4 ([58]https://locuszoom.sph.umich.edu/), SNPs in genomic risk loci are color-coded as a function of their r^2 to the index SNP rs11465988 in the locus, while SNPs with missing LD information are shown in grey. For the circos plot generated via FUMA^[59]51 v1.3.5e ([60]https://fuma.ctglab.nl/), the outer most layer is Manhattan plot and the middle layer highlights genomic risk loci (as defined by FUMA^[61]51 using minimum P-value of lead SNPs of 1 × 10^–5 and default values for other parameters) in blue, while the inner most layer highlights eQTLs and/or chromatin interactions. Only SNPs with p < 0.05 are displayed in the outer ring. SNPs in genomic risk loci are color-coded as a function of their maximum r^2 to the one of the independent significant SNPs in the locus. The rsID of the top SNPs in each risk locus are displayed in the most outer layer. For the inner most layer, if the gene is mapped only by chromatin interactions or only by eQTLs, it is colored orange or green, respectively. It is colored red when the gene is mapped by both. The other two GWAS for responder and remission status, respectively, did not yield any genome-wide significant finding (Supplemental Figs. [62]3A,B for Manhattan plots; Supplemental Figs. [63]1B,C for QQ plots, λ = 1.028 and 0.997, respectively). Nevertheless, a suggestive signal that merits comment was identified in chromosome 8 (rs4739050, nominal p = 6.06 × 10^–8, β = 7.5, SE = 1.4 for percentage change in MADRS score; rs10957273, nominal p = 8.07 × 10^–7, OR = 0.3, SE = 0.2 for responder status) from both the continuous endpoint GWAS and the responder status GWAS. Rs4739050 is an expression quantitative trait locus (eQTL) for gamma-glutamyl hydrolase (GGH) based on eQTLGen (p[eQTL] = 4.24 × 10^–9). A full list of suggestive associations with p-values less than 1 × 10^–4 is provided in Supplemental Table [64]3. Gene-level association analysis revealed one significant gene NME/NM23 family member 7 (NME7, p = 1.73 × 10^–6, Supplemental Fig. [65]4A) for the percentage change in MADRS score. In the percent change in MADRS score GWAS, a pathway enrichment analysis revealed suggestive enrichments of genes involved in the negative regulation of glucocorticoid metabolic process (nominal p = 3.53 × 10^–5) and neuronal action potential (nominal p = 0.0001). Pathway enrichment analysis also revealed suggestive (p-values listed are nominal) enrichments of genes involved in synaptic vesicle clustering (p = 4.33 × 10^–5), negative regulation of glucocorticoid metabolic process (p = 5.48 × 10^–5), regulation of synaptic vesicle clustering (p = 6.13 × 10^–5), anterior posterior axon guidance (p = 0.0002), and netrin mediated repulsion signals (p = 0.0002) in the responder status GWAS, and in the negative regulation of extrinsic apoptotic signaling pathway (p = 4.04 × 10^–5), NF-κB canonical pathway (p = 5.90 × 10^–5), stress pathway (p = 0.0002), and TNFR1 induced proapoptotic signaling (p = 0.0003) in the remission status GWAS (Supplemental Table [66]4). We did not identify an association between the change in MADRS score and the BDNF Val66Met polymorphism in the current study (p > 0.05). After applying corrections for multiple testing, none of the associations between esketamine’s antidepressant efficacy and the PRS genetic loading for psychiatric conditions or symptom profiles was significant at the study-wide level (which required p < 0.0004 for significance). In Table [67]3 and Supplemental Fig. [68]5 we list suggestive associations observed in these analyses, however, along with their nominal p-values. Thus we observed suggestive (i.e., p-values listed are nominal) negative correlations between the depressive symptom PRS^[69]15 (p = 0.001, Table [70]3 and Supplemental Fig. [71]5) and the esketamine treatment response outcome as measured by percentage change from baseline in the MADRS score at the end of four week treatment period (Table [72]3). In addition, the depressive symptom PRS (p = 0.004) displayed suggestive positive correlations with esketamine responder status. Lastly, depressive symptoms PRS (p = 0.002) and insomnia^[73]16 PRS (p = 0.003) exhibited suggestive positive correlations with esketamine remission status. Table 3. Polygenic Risk Score association with esketamine treatment outcome. Threshold r^2[PRS] r^2[Full] r^2[Null] Standardized coefficient Standard error p Number of SNP Base GWAS References