Pharmacogenomics of selective serotonin reuptake inhibitor treatment for major depressive disorder: genome-wide associations and functional genomics.

A genome-wide association (GWA) study of treatment outcomes (response and remission) of selective serotonin reuptake inhibitors (SSRIs) was conducted using 529 subjects with major depressive disorder. While no SNP associations reached the genome-wide level of significance, 14 SNPs of interest were identified for functional analysis. The rs11144870 SNP in the riboflavin kinase (RFK) gene on chromosome 9 was associated with 8-week treatment response (odds ratio (OR)=0.42, P=1.04 × 10⁻⁶). The rs915120 SNP in the G protein-coupled receptor kinase 5 (GRK5) gene on chromosome 10 was associated with 8-week remission (OR=0.50, P=1.15 × 10⁻⁵). Both SNPs were shown to influence transcription by a reporter gene assay and to alter nuclear protein binding using an electrophoretic mobility shift assay. This report represents an example of joining functional genomics with traditional GWA study results derived from a GWA analysis of SSRI treatment outcomes. The goal of this analytical strategy is to provide insights into the potential relevance of biologically plausible observed associations.

Introduction

Major depressive disorder (MDD) is a serious and prevalent psychiatric illness.[1] Selective serotonin reuptake inhibitors (SSRIs) are an effective treatment for MDD, but treatment response is highly variable.[2] Consequently, the identification of pharmacogenetic predictors of variable drug response phenotypes has become a major objective of MDD research.[3-5]

Genome-wide association (GWA) studies provide an agnostic approach which can identify novel genetic variants that may contribute to variation in drug response phenotypes.[6] Three GWA studies of SSRI antidepressant outcomes have been reported, but none have identified associations of SNPs with treatment outcomes that reached genome-wide statistical significance or were validated in a separate population.[7-9] This failure to demonstrate genome-wide significance may be related to (1) insufficient statistical power to detect small effect sizes for associated variants, (2) inconsistently defined phenotypes, (3) confounding effects from non-genetic factors, or (4) inconsistencies in treatment protocols.[10, 11] However, the conduct of functional genomic studies of the SNPs identified by GWA analyses provides a novel strategy to identify genomic variants that may influence antidepressant treatment.

The present study is a report of GWA analyses of MDD patients who participated in a single-site pharmacogenomic study of SSRI therapy. Following the GWA analyses, a reporter gene assay and an electrophoretic mobility shift assay were performed to determine whether any identified candidate SNPs were associated with a change in function.

Materials/Subjects and Methods

Study Design and Samples

The Mayo Clinic Pharmacogenomic Research Network Antidepressant Medication Pharmacogenomic Study (PGRN-AMPS) was supported by the NIGMS-Pharmacogenomics Research Network (PGRN), which has been described elsewhere.[12, 13] The PGRN-AMPS is an ongoing eight week outpatient SSRI clinical trial that was performed at the Mayo Clinic in Rochester, MN. Patients enrolled in the study met diagnostic criteria for MDD without psychosis or mania and had a 17-item Hamilton Depression Rating Scale (HAMD-17) score ≥14. The study was designed with inclusion and exclusion criteria similar to those used in the Sequenced Treatment Alternatives to Relieve Depression study (STAR*D).[14] Potential study subjects taking an antidepressant, antipsychotic or mood stabilizing medication were not eligible for enrollment. Patients with MDD initially received either 10 mg of escitalopram or 20 mg of citalopram. SSRI efficacy was determined using the 16-item Quick Inventory of Depressive Symptomatology (QIDS-C16) scores after four weeks and then eight weeks of SSRI therapy. At four weeks after the initiation of treatment, the dose could be increased to 20 mg of escitalopram or 40 mg of citalopram after a clinical assessment of the subject. Unless there was a contraindication, dose was increased if the QIDS-C16 score at the follow-up visit was ≥9, and possibly following a clinical evaluation if the score was between 6 and 8. The dose could also be decreased, or treatment could be discontinued, if a patient developed persistent side effects. Blood samples were obtained at baseline for DNA extraction, and then again at weeks four and eight for assays of drug and metabolite levels. All patients provided written informed consent. The study protocol was approved by the Mayo Clinic Institutional Review Board.

Treatment Outcomes

The two primary outcome phenotypes in the GWA analyses were “response” (defined as ≥50% reduction in QIDS-C16 score from baseline to the last visit) and “remission” (defined as a QIDS-C16 score of ≤5 at the last visit). For each of these two outcomes, analyses were performed using two strategies. The primary analyses included only subjects that were evaluated at the eight-week visit. The secondary analyses were performed with outcomes based on the final visit QIDS-C16 scores, referred to as the “last visit” assessment. These analyses included subjects who had completed the full eight-week study as well as those who dropped out of the study prior to the eight-week assessment.

Genotyping and QC

DNA from 529 patients was genotyped by the RIKEN Center for Genomic Medicine (Yokohama, Japan) using Illumina Human610-Quad BeadChips (Illumina, San Diego, CA). Taqman genotyping assays (Applied Biosystems, Foster City, CA) were used to perform genotyping for the replication study using STAR*D samples, and the Illumina Goldengate platform was used to genotype selected top SNPs identified in the published STAR*D GWAS [7] using Mayo PGRN-AMPS DNA samples.

Quality control assessments included overall genotype concordance rates based on duplicate sample genotyping and Mendelian inheritance checks based on genotyping of a CEPH trio of two parents and their child. For each SNP, the minor allele frequency, call rate, and departure from Hardy-Weinberg Equilibrium were evaluated. Observed call rates, total heterozygosity and inbreeding coefficients were assessed for each subject using PLINK.[15] Sex-checks based on X-chromosome heterozygosity were performed, and tests of identity-by-descent were used to identify potentially related subjects.

More than 97% of study participants were of white non-Hispanic (WNH) ancestry. Consequently, genotype-phenotype association analyses were restricted to the WNH subjects. A subset of 4,855 independent SNPs with low local linkage disequilibrium (r2 < 0.063) was used to verify self-reported ancestry of the subjects using the software STRUCTURE.[16] In addition to the study population, 287 DNA samples from the “Human Variation Panel” of lymphoblastoid cell lines (sample sets HD100CAU, HD100AA and HD100CHI) obtained from the Coriell Institute (Camden, NJ) were included in this analysis. Probabilities of membership in each of three known ancestral groups of the cell lines were calculated for each study subject. These analyses were used to verify self-reported race and to assign race to subjects with unknown race. 509 subjects identified themselves as WNH. A further five subjects who self-identified their race as “other” or “multiracial” were found to have >75% European ancestry based on the STRUCTURE analysis. Thus, 514 WNH subjects were included in subsequent analyses.

Statistical Analyses

Prior to genotype-phenotype association analyses, EigenStrat software was used to determine eigenvalues for the SNP correlation matrix.[17] Eigenvalues that differed statistically from zero were determined based on Tracy-Widom test p-values, and the corresponding eigenvectors were used as covariates in the genetic analyses to account for differences in ancestry within the WNH set of subjects.

Of the 514 WNH subjects, 15 patients were found to have very low blood drug levels, suggesting non-adherence, and were excluded from analysis, resulting in 499 subjects who were included in the primary GWA analyses. Logistic regression was used to test for associations between the binary outcome variables, response and remission, and the genotype of each SNP. Genotypes were coded as the “dose” of the minor allele. These analyses assume log-additive allele effects on response or remission. All GWA analyses were adjusted for eigenvectors that captured population stratification. Gender and age (categorized into four quartiles) were also considered as covariates for the GWA analyses. Age was not significantly associated with response or remission (p > 0.10 for both outcomes). Gender was not associated with response (p > 0.10 for both 8-week and final-visit response), but was marginally associated with remission (p = 0.057 for 8-week remission, p = 0.036 for final-visit remission). We therefore performed response and remission GWA analyses both with and without gender as a covariate. As the results were nearly identical, we report only findings from the analyses not adjusted for gender.

Analyses were performed in R (http://www.R-project.org), SAS (SAS Institute Inc.) and PLINK.[15] Odds ratios (ORs), 95% confidence intervals (CIs), and p-values were calculated based on the logistic regression models. P-values < 10−7 were considered statistically significant at the genome-wide level.

Prior to performing the GWA analyses, the power to detect SNP effects was estimated, assuming 500 subjects with 30% failing to achieve remission. It was determined that with a type 1 error rate of 10−7, the data would provide >80% power to detect ORs of at least 1.8 for common alleles (MAF 0.10–0.20).

Functional Genomics

On review of the GWA analyses, 14 candidate SNPs were selected based on having some evidence of an association (p < 9.0 × 10−5) and being located within or near a candidate gene. SNPs in locations that could potentially influence gene transcription were preferentially selected. Reporter gene assays were used to demonstrate the effects of these SNPs on transcription and electrophoretic mobility shift (EMS) assays were used to assess DNA-protein binding.

For the reporter gene assays, 200 to 300 bp DNA sequences that included a candidate SNP were cloned into a pGL3-promoter luciferase reporter vector that contained an SV40 promoter upstream of the luciferase gene (Promega Corporation, Fitchburg, WI). One microgram of each reporter gene construct was co-transfected with 20 ng of the pRL-TK renilla luciferase vector as a control for transfection efficiency into a human neuroblastoma cell line, SK-N-BE(2) and two human glioblastoma cell lines, U-87 MG and U251 (American Type Culture Collection, ATCC, Manassas, VA), followed by dual-luciferase assays performed 24 hrs after transfection (Promega). Two independent transfections were performed for each reporter gene construct, with triplicate transfections for each construct in each experiment. Values for relative activity were expressed as percentages of the pGL3-promoter activity for vector without an insert. Comparisons were then made between pGL3 reporter gene constructs containing wild type and variant nucleotide for each candidate SNP. DNA samples that were used to amplify SNP regions were selected from the Coriell “Human Variation Panel” DNA samples (Camden, NJ).[18] Sequences of primers used to amplify genomic regions containing the selected SNPs are listed in Supplemental Table 1.

EMS assays were used to determine whether a nucleotide change altered the ability of oligonucleotides to bind nuclear proteins. Nuclear extracts were prepared from SK-N-BE(2), U87 MG and U251 cells using the Nuclear Extract Kit (Active Motif, Carlsbad, CA). Biotin-labeled and unlabeled oligonucleotides were synthesized by IDT Integrated DNA Technologies (Coralville, IA). EMS assays were performed using the LightShift® Chemiluminescent EMSA Kit (PIERCE, Rockford, IL). Oligonucleotide sequences used to perform the EMS assays are listed in Supplemental Table 2.

Replication

The 14 candidate SNPs included in the functional genomic experiments were genotyped in the STAR*D DNA samples. In addition, the 25 SNPs reported in the STAR*D GWA study[7] to have the greatest likelihood of being associated with treatment outcomes were genotyped in the PGRN-AMPS samples. The results of these analyses are described in the Supplemental Materials.

Results

Demographic and Clinical Characteristics of the Study Subjects

Of the 529 participants in the Mayo PGRN-AMPS study, 15 subjects who were not identified as white and non-Hispanic and 15 non-adherent subjects were excluded, resulting in 499 WNH subjects for the secondary “last visit” GWA analyses. The primary “eight week” GWA analyses were conducted using a subset of 398 WNH subjects who were protocol adherent. Demographic information, baseline clinical characteristics and treatment outcomes for both the eight week and the last visit cohorts are summarized in Table 1.

Table 1

Baseline clinical and demographic characteristics and outcomes for medication-adherent white non-Hispanic subjects included in the last-visit and eight-week analyses

Characteristic and Measure	Last visit (N = 499)	Eight weeks (N = 398)
Demographics
Age	39.91 (±13.8)	40.64 (±13.5)
Education (years)	14.68 (±2.4)	14.91 (±2.3)
Female gender	312 (62.5%)	255 (64.1%)
Clinical Characteristics
QIDS-C Baseline	15.08 (±3.47)	14.99 (±3.31)
QIDS-C Week Four	8.46 (±4.41)	8.34 (±4.41)
QIDS-C Week Eight	6.24 (±4.05)	6.24 (±4.05)
Age at Onset of First Episode of Depression	24.05 (±13.82)	24.70 (±14.11)
Recurrent Depression (>1 episode)	483 (96.8%)	386 (96.9%)
Baseline Medication
Citalopram	155 (31.2%)	124 (31.2%)
Escitalopram	342 (68.8%)	274 (68.8%)
Last known dose (mg)
Citalopram	28.10 (±10.02)	30.28 (±9.76)
Escitalopram	14.56 (±5.06)	15.75 (±4.96)
Outcomes
Remitter (QIDS ≤5)	206 (41.3%)	198 (49.7%)
Response (% reduction QIDS ≥50%)	287 (57.5%)	274 (68.8%)

Values are mean (± 1 standard deviation) or number of subjects, N (%).

Quality Control

Quality control analyses were performed for 592,236 genotyped SNPs. The genotypes had very high concordance rates in duplicate samples (> 99.99%) and very low rates of Mendelian inheritance errors (< 0.01%). 30,700 SNPs failed genotyping. Another 28,610 SNPs had allele frequencies below the preset threshold of 0.01 and, consequently, were excluded from analysis due to low power. An additional 49 SNPs were excluded due to significant departures from Hardy-Weinberg Equilibrium (p < 10−6). The remaining SNPs (n = 532,877) had call rates exceeding 98%. One person who was reported to be female appeared to be genetically male based on no heterozygosity of X-chromosome SNPs. Identity-by-descent analysis revealed several pairs of subjects that appear to be related (two likely parent-child pairs, one possible pair of siblings, and three pairs that may be more distant relatives). Genome-wide heterozygosity was within the expected range.

Genetic Association Analysis of SSRI Treatment Outcomes

Separate GWA analyses based on the two primary outcomes of response and remission were conducted for both the eight week protocol adherent sample (n = 398) and the last visit sample (n = 499). Four eigenvectors identified through EigenStrat analysis were used as adjusting covariates in the analyses. While none of these analyses demonstrated genome wide significance (p < 10−7), the 25 SNPs with the most significant associations with the four outcome phenotypes are listed in Table 2. Manhattan plots of p-values from the analyses for eight week outcomes are shown in Figure 1. Manhattan plots for last visit outcomes are included in Figure S1.

Table 2

Top 25 association results for each of the four outcomes: response at eight weeks and last-visit and remission at eight weeks and last visit (based on logistic regression adjusting for four eigenvectors)

				MAF		Association with Response/Remission
Rank	SNP	Chr.	Position	Responder/Remitter	Non-responder/Non-remitter	OR	P-value	Gene	Gene Region	Distance to Gene (bp)
Response at Eight Weeks
1	rs11144870a	9	78194033	0.22	0.39	0.43	1.04E-06	RFK	intron	−641
2	rs11144905a	9	78255119	0.26	0.43	0.45	1.43E-06	GCNT1	flanking_5UTR	−8847
3	rs9775457	9	139386750	0.33	0.17	2.63	3.15E-06	NA	NA	NA
4	rs785916a	9	78290826	0.21	0.36	0.46	7.50E-06	GCNT1	intron	−14826
5	rs1045074	6	4899709	0.47	0.30	2.13	1.18E-05	CDYL	3UTR	[419/1068]
6	rs3906767	2	12222673	0.13	0.25	0.40	1.30E-05	LPIN1	flanking_3UTR	−337687
7	rs3954862	2	12222204	0.13	0.25	0.40	1.30E-05	LPIN1	flanking_3UTR	−337218
8	rs731395	2	12214033	0.13	0.26	0.41	1.58E-05	LPIN1	flanking_3UTR	−329047
9	rs1869733	8	63970954	0.13	0.25	0.41	2.19E-05	FAM77D	intron	−22612
10	rs1428241	5	25579023	0.12	0.23	0.41	3.10E-05	CDH10	flanking_5UTR	−898355
11	rs2837754	21	40912929	0.05	0.14	0.32	3.72E-05	DSCAM	intron	−73677
12	rs4711589	6	39338549	0.25	0.40	0.50	3.82E-05	KCNK5	flanking_5UTR	−33345
13	rs2283859	22	28234519	0.20	0.33	0.47	3.99E-05	C22orf19	intron	−10
14	rs1182864	6	39338285	0.40	0.25	2.07	4.04E-05	KCNK5	flanking_5UTR	−33081
15	rs1455583	8	63975240	0.13	0.25	0.42	4.06E-05	FAM77D	intron	−18326
16	rs6493270	15	45400695	0.12	0.23	0.43	4.44E-05	SEMA6D	flanking_5UTR	−397283
17	rs4649163	1	24252188	0.24	0.38	0.50	4.70E-05	MYOM3	flanking_3UTR	−2934
18	rs2595372	18	27492355	0.16	0.30	0.49	5.07E-05	B4GALT6	intron	−305
19	rs680448	9	78316773	0.54	0.40	1.93	5.64E-05	GCNT1	flanking_3UTR	−4621
20	rs714679	11	44610240	0.38	0.54	0.53	5.71E-05	CD82	flanking_3UTR	−12325
21	rs4242365	8	125307549	0.25	0.40	0.51	5.90E-05	C8orf78	flanking_5UTR	−54605
22	rs10485285	6	81535326	0.29	0.45	0.52	6.75E-05	BCKDHB	flanking_3UTR	−422620
23	rs13014477	2	29649492	0.33	0.19	2.15	7.50E-05	ALK	intron	−41006
24	rs16987839	22	28313138	0.15	0.28	0.49	7.93E-05	NIPSNAP1	flanking_5UTR	−5810
25	rs3804024	21	40915577	0.06	0.14	0.35	8.10E-05	DSCAM	intron	−71029
Response at Last Visit
1	rs2248399a	13	105440474	0.12	0.23	0.46	1.39E-05	DAOA	flanking_3UTR	−499090
2	rs6501039	16	7920533	0.34	0.22	1.90	2.06E-05	A2BP1	flanking_3UTR	−218033
3	rs2831007	21	27863838	0.31	0.19	1.92	2.08E-05	C21orf94	flanking_5UTR	−443734
4	rs13014477	2	29649492	0.34	0.22	1.89	2.64E-05	ALK	intron	−41006
5	rs1751915	10	21484609	0.15	0.26	0.49	2.90E-05	C10orf113	flanking_5UTR	−9176
6	rs3811541	2	127129605	0.30	0.43	0.55	3.14E-05	GYPC	flanking_5UTR	−549
7	rs12450499	17	10109041	0.11	0.20	0.45	3.36E-05	MYH13	flanking_3UTR	−35867
8	rs2248174a	13	105442451	0.11	0.21	0.46	3.48E-05	DAOA	flanking_3UTR	−501068
9	rs11144870a	9	78194033	0.22	0.33	0.55	3.70E-05	RFK	intron	−641
10	rs7378714	5	111921446	0.29	0.18	1.96	4.64E-05	FLJ11235	flanking_3UTR	−136878
11	rs4620377	9	14623299	0.37	0.50	0.58	4.89E-05	ZDHHC21	intron	−6595
12	rs12328118	2	82811461	0.26	0.15	2.02	4.93E-05	LOC388965	flanking_5UTR	−1559857
13	rs4835490	4	149309776	0.24	0.14	2.07	4.95E-05	NR3C2	intron	−14274
14	rs3762876	4	95808500	0.42	0.30	1.75	5.08E-05	PDLIM5	flanking_3UTR	−100
15	rs3761555	23	122144118	0.22	0.38	0.52	5.16E-05	GRIA3	flanking_5UTR	−1659
16	rs2549151	16	81417750	0.40	0.53	0.58	5.21E-05	CDH13	intron	−31718
17	rs6493270	15	45400695	0.13	0.23	0.50	5.58E-05	SEMA6D	flanking_5UTR	−397283
18	rs1902763	12	50943932	0.21	0.12	2.19	5.62E-05	KRT7	flanking_3UTR	−14956
19	rs6104012	20	43167985	0.15	0.25	0.51	5.91E-05	WFDC5	flanking_3UTR	−3522
20	rs11692714	2	82835165	0.25	0.14	2.03	6.20E-05	LOC388965	flanking_5UTR	−1536153
21	rs17022383	2	82845146	0.25	0.14	2.03	6.20E-05	NA	NA	NA
22	rs9828629	3	71613036	0.41	0.29	1.78	6.36E-05	FOXP1	intron	−12231
23	Rs1945300	12	50953187	0.22	0.12	2.16	6.80E-05	KRTHB1	flanking_3UTR	−12777
24	rs1998560a	13	105441526	0.15	0.26	0.51	6.80E-05	DAOA	flanking_3UTR	−500134
25	rs2608017	12	50962184	0.21	0.12	2.20	6.81E-05	KRTHB1	flanking_3UTR	−3780
Remission at Eight Weeks
1	rs1379887a	6	78407561	0.27	0.42	0.49	9.05E-06	HTR1B	flanking_5UTR	−177722
2	rs915120a	10	121180103	0.32	0.48	0.50	1.15E-05	GRK5	intron	−163
3	rs10899655	11	78994647	0.11	0.23	0.41	1.62E-05	NARS2	flanking_5UTR	−1031280
4	rs7738598a	6	78671581	0.52	0.36	1.87	1.67E-05	HTR1B	flanking_5UTR	−441742
5	rs898040a	6	78408635	0.27	0.41	0.51	1.69E-05	HTR1B	flanking_5UTR	−178796
6	rs790575	12	91055739	0.17	0.30	0.47	2.03E-05	BTG1	flanking_3UTR	−5291
7	rs6793050	3	66468459	0.09	0.19	0.40	3.58E-05	NA	NA	NA
8	rs978814	6	81016631	0.37	0.52	0.54	3.82E-05	BCKDHB	intron	−22940
9	rs12413418	10	42493575	0.03	0.10	0.23	3.89E-05	ZNF11B	flanking_5UTR	−39577
10	rs16840617	4	7661440	0.40	0.26	1.91	4.05E-05	SORCS2	intron	−29515
11	rs12254134a	10	121161872	0.35	0.49	0.52	4.24E-05	GRK5	intron	−10796
12	rs12963758	18	7814091	0.10	0.03	4.35	4.44E-05	PTPRM	intron	−49822
13	rs1924939	13	71998057	0.29	0.44	0.54	4.67E-05	FLJ22624	flanking_5UTR	−201986
14	rs1359873	13	74194810	0.36	0.24	2.01	4.76E-05	LOC400145	flanking_5UTR	−303557
15	rs9436797	1	68340006	0.16	0.28	0.47	4.81E-05	GPR177	intron	−2982
16	rs3812121	6	81025446	0.37	0.52	0.55	4.94E-05	BCKDHB	intron	−14125
17	rs6414150	2	196603015	0.30	0.18	2.01	4.98E-05	DNAH7	intron	−1999
18	rs9343638	6	78438553	0.29	0.43	0.53	4.99E-05	HTR1B	flanking_5UTR	−208714
19	rs757240	4	7658657	0.39	0.26	1.89	5.33E-05	SORCS2	intron	−32298
20	rs2252523	16	1937615	0.25	0.39	0.53	5.49E-05	RPL3L	intron	−233
21	rs9542990	13	71993272	0.22	0.35	0.52	6.20E-05	FLJ22624	flanking_5UTR	−206771
22	rs9352817	6	81082598	0.38	0.52	0.55	6.70E-05	BCKDHB	intron	−27502
23	rs1936158	6	78241094	0.48	0.35	1.87	6.89E-05	HTR1B	flanking_5UTR	−11255
24	rs351294	5	154428932	0.05	0.13	0.33	7.40E-05	MRPL22	flanking_3UTR	−102208
25	rs17777705	6	78551483	0.15	0.27	0.49	7.62E-05	HTR1B	flanking_5UTR	−321644
Remission at Last Visit
1	rs7439567a	4	138684292	0.32	0.46	0.51	3.11E-06	PCDH18	flanking_5UTR	−11213
2	rs1379887a	6	78407561	0.27	0.41	0.52	7.88E-06	HTR1B	flanking_5UTR	−177722
3	rs9542990	13	71993272	0.22	0.35	0.53	1.60E-05	FLJ22624	flanking_5UTR	−206771
4	rs7738598a	6	78671581	0.52	0.38	1.76	1.66E-05	HTR1B	flanking_5UTR	−441742
5	rs898040a	6	78408635	0.27	0.40	0.54	1.81E-05	HTR1B	flanking_5UTR	−178796
6	rs9761827a	4	138635961	0.30	0.43	0.55	2.27E-05	PCDH18	flanking_3UTR	−24324
7	rs1924939	13	71998057	0.30	0.44	0.56	2.36E-05	FLJ22624	flanking_5UTR	−201986
8	rs7058787	23	45583663	0.45	0.27	1.94	2.44E-05	LOC100128442	gene	0
9	rs12963758	18	7814091	0.10	0.03	3.36	2.56E-05	PTPRM	intron	−49822
10	rs1945300	12	50953187	0.24	0.14	2.11	3.28E-05	KRTHB1	flanking_3UTR	−12777
11	rs12809631	12	129611143	0.17	0.28	0.50	3.33E-05	RIMBP2	flanking_5UTR	−42780
12	rs2385296	8	124655486	0.53	0.40	1.75	3.42E-05	FBXO32	flanking_5UTR	−32859
13	rs1902763	12	50943932	0.23	0.13	2.11	3.46E-05	KRT7	flanking_3UTR	−14956
14	rs3752183	19	61804990	0.46	0.33	1.78	3.54E-05	ZNF71	intron	−219
15	rs4947035	6	110352448	0.14	0.06	2.58	3.69E-05	GPR6	flanking_5UTR	−54539
16	rs3807736	7	77571209	0.42	0.55	0.58	4.46E-05	MAGI2	intron	−21738
17	rs17541104	15	91664471	0.37	0.25	1.80	4.48E-05	UNQ9370	flanking_3UTR	−112190
18	rs6580873	12	50968192	0.21	0.11	2.14	4.61E-05	KRTHB1	coding	[157/7]
19	rs9294091	6	78712418	0.37	0.25	1.80	5.73E-05	HTR1B	flanking_5UTR	−482579
20	rs2085144	8	99022252	0.16	0.27	0.51	6.09E-05	MATN2	intron	−929
21	rs1877531	6	78618795	0.40	0.27	1.76	6.24E-05	HTR1B	flanking_5UTR	−388956
22	rs1655645	15	27479972	0.50	0.38	1.73	6.62E-05	NDNL2	flanking_5UTR	−130663
23	rs2336413	2	205839562	0.42	0.30	1.78	6.65E-05	ALS2CR19	intron	−20771
24	rs9343638	6	78438553	0.29	0.42	0.57	7.17E-05	HTR1B	flanking_5UTR	−208714
25	rs12722898	1	26593874	0.48	0.35	1.72	7.28E-05	LIN28	flanking_5UTR	−15982

Abbreviations: SNP, single-nucleotide polymorphism; Chr., chromosome; MAF, minor allele frequency; OR, odds ratio

SNPs selected for functional genomic studies

Figure 1

Manhattan plots of −log10(p-values) from logistic regression adjusted for eigenvectors versus chromosomal position of SNPs. (a) Manhattan plot of −log10(p-values) for GWA analysis of eight-week SSRI response outcome. (b) Manhattan plot of −log10(p-values) for GWA analysis of eight-week SSRI remission outcome. Red dots represent SNPs with p-values ≤ 10−5 and blue dots represent SNPs with p-values ≤ 10−4.

Eight Week Response Analyses

GWA analysis of the eight week response outcome identified a region on chromosome 9 containing three interesting SNPs with p-values < 8.0 × 10−6. The rs11144870 SNP is located in intron 2 of the riboflavin kinase(RFK) gene. The rs11144905 and rs785916 SNPs were mapped to introns in the adjacent glucosaminyl N-acetyl transferase 1 (GCNT1) gene and were in linkage disequilibrium with rs11144870 (r2 = 0.78 and 0.73, respectively).

Eight Week Remission Analyses

Three SNPs, rs1379887, rs7738598 and rs898040, that were located in the 5′ UTR region of the 5-hydroxytryptamine serotonin receptor 1B gene (HTR1B) were associated with eight week remission at the p < 2 × 10−5 level. The rs915120 and rs12116187 SNPs in the G protein-coupled receptor (GPCR) kinase gene 5 (GRK5) on chromosome 10 were also identified as candidate SNPs of interest.

Analysis of Secondary Outcomes: Last Visit Response and Remission

The most significant association in the last visit response analysis was obtained for rs2248399, which is located in an intergenic region on chromosome 13 (p = 1.0 × 10−5). This SNP is located approximately 500 kb from the D-amino acid oxidase activator gene (DAOA). The three SNPs near the 5-hydroxytryptamine (serotonin) receptor 1B gene (HTR1B), rs1379887, rs7738598 and rs898040, that were associated with remission in the eight week analysis, were also identified in the last visit remission analysis (p-values = 7.9 × 10−6, 1.7 × 10−5 and 1.8 × 10−5, respectively). The last visit remission analysis also identified the rs7439567 (p = 3.1 × 10−6) and the rs9761827 (p = 2.3 × 10−5) SNPs in the protocadherin-18 gene (PCDH18) on chromosome 4.

As described previously, 14 candidate SNPs were chosen for functional studies based on the significance level of their association with treatment outcomes and their location within or near a candidate gene of interest (Figure 2A). The reporter gene assays identified significant differences in the function of WT and variant SNP sequences for five of these SNPs based on differences in relative luciferase activity between pGL3 constructs that contained the sequences in all three CNS-derived cell lines (Figure 2B). When cloned into the pGL3-Promoter vector, sequences containing the SNPs of interest appeared to function as either an enhancer element that significantly increased promoter activity when compared with activity for the pGL3-Promoter construct alone (p-values < 0.05 in both cell lines) or a silencer element depending on the cell line used. Results for the neuroblastoma SK-N-BE(2) cells are shown in Figure 2b, and results for U-87 MG and U251 cell lines in Figure S2a. The patterns of differences in luciferase activity between WT and variant alleles were consistent for SNPs rs11144870 (RFK), rs11144905 (GCNT1), rs2248399 (DAOA), rs915120 (GRK5) and rs12254134 (GRK5) in all 3 cell lines tested (p-values < 0.01).

Figure 2

Functional genomics of candidate SNPs. (a) A table to summarize the 14 SNPs that were selected for functional genomic studies. (b) Results of dual-luciferase reporter gene assays for SNPs rs11144870, rs11144905, rs2248399, rs915120 and rs12254134, performed in the neuroblastoma SK-N-BE(2) cell line. Each bar represents the average of relative luciferase activity reported as a % of the pGL3-Promoter construct activity obtained during 6 independent transfections (mean ± S.E.M). ** represents p-values < 0.01 compared to values for pGL3-WT activity. (b) Electrophoretic mobility shift (EMS) assays for SNPs rs11144870, rs785916 (GCNT1), rs2248399 (DAOA), rs2248174 (DAOA), rs1998560 (DAOA), rs915120 (GRK5) and rs12254134 (GRK5) with nuclear extract prepared from SK-N-BE(2) cells.

EMS assays identified seven SNPs that altered the ability of these sequences to bind nuclear proteins in SK-N-BE(2) cells. Figure 2c shows the results for EMS assays for these seven SNPs, including rs11144870, rs785916 (GCNT1), rs2248399, rs2248714, rs1998560 (DAOA), rs915120 and rs12254134. All of these SNPs showed significant differences in nuclear protein binding patterns between WT and variant sequence oligonucleotides in SK-N-BE(2) cells. Similar patterns were observed for all of these SNPs in U-87 MG and U251 cells (Figure S2b). HTR1B SNP rs7738598 displayed altered DNA-protein binding only in the two glioblastoma cells (Figure S2b), but not in the neuroblastoma cells.

Four candidate SNPs, rs11144870 (RFK), rs915120 (GRK5), rs12254134 (GRK5) and rs2248399 (DAOA) showed differences between WT and variant SNP sequences in both reporter gene and EMS assays performed with all three cell lines.

The SNPs selected for functional studies based on GWA analyses of the PGRN-AMPS outcomes data were also genotyped using DNA from the STAR*D study. Association analyses were performed for the same four treatment outcomes with 12 of the 14 SNPs, as two SNPs failed genotyping, using DNA from WNH STAR*D subjects who had HAMD-17 scores ≥14 prior to SSRI therapy. None of the associations for these SNPs were statistically significant in the STAR*D samples (Supplemental Table S7). The top 25 SNPs reported for the STAR*D GWA analyses[7] were also evaluated for possible association with treatment response in the PGRN-AMPS samples. None of these associations were statistically significant (Supplemental Table S8).

Discussion

Four GWA studies of antidepressant treatment outcomes have been conducted including the PGRN-AMPS study.[7-9] None of these GWA analyses was able to identify associations that reached genome-wide statistical significance. Furthermore, no specific SNPs or genes have been associated with outcomes across any two of these studies. Given that replication has not been reported, functional validation studies provide an alternative strategy to take the next step toward deciphering the underlying biology of candidate SNPs that might influence antidepressant response.

The association of SNP rs11144870 in RFK with SSRI outcomes is intriguing as riboflavin kinase is an essential enzyme that catalyzes the phosphorylation of riboflavin (vitamin B2) to form flavin mononucleotide (FMN). This is a critical step for both vitamin B2 metabolism and flavin cofactor synthesis. Riboflavin and its flavin cofactors influence the folate and methionine cycles since riboflavin functions as a cofactor for methylene tetrahydrofolate reductase (MTHFR). MTHFR is the enzyme that converts 5,10-methylenetetrahydrofolate to 5-methyltetrahydrofolate.[19] Previous studies have suggested that components of the folate and methionine cycles may be involved in increasing the risk for developing MDD and might influence treatment outcomes.[20-24] Since insufficient dietary B vitamins, including riboflavin, have been associated with depressive symptoms,[25] an alternation in the transcription of RFK might result in elevated levels of RFK protein which might indirectly influence the intensity of depressive symptoms and the effects of SSRI therapy.

The GCNT1 gene is a member of the beta-1,6-N-acetylglucosaminyltransferase gene family that has not previously been shown to be associated with MDD treatment outcomes. However, GCNT1 is highly expressed in brain, and further research focusing on variants in GCNT1 would be of interest (see the Nervous System database, http://www.itb.cnr.it/gncdb/).

A SNP in an intergenic region near the DAOA gene (rs2248399) was identified by these analyses and may be potentially functionally significant based on the reporter gene assay. This SNP and two other DAOA SNPs were also shown to have the potential to affect the binding of nuclear proteins. None of these SNPs have been included in previous candidate gene studies of DAOA for schizophrenia or bipolar disease.[26-31]

Six SNPs that mapped to an intergenic region ~150 to 500 kb distant from the HTR1B gene were associated with eight-week and last visit remission (see Table 2 and Supplemental Tables S5 and S6). The HTR1B gene has been reported to be associated with psychiatric phenotypes as well as response to SSRI treatment.[32-35] EMS assays were performed with these three SNPs and with rs7738598 were found to display a difference between WT and variant sequences in nuclear protein binding in the two glioblastoma cells that were tested.

Two of the intronic SNPs in GRK5 (rs915120 and rs12254134) that were associated with remission also altered function. A different member of the G protein-coupled receptor kinase family, GRK2, has been found to be upregulated during antidepressant treatment[36] and GRK5 has been shown to regulate GPCR receptors such as the β1-adrenergic receptor[37] and the dopamine D1A receptor.[38] GRK5 is highly expressed in many tissues, including human heart and brain.[39, 40] A single functional GRK5 polymorphism, rs17098707, that results in a Gln41Leu change in amino acid sequence has been reported to regulate cardiac function.[41, 42] The intronic SNPs identified in this GWA study are not in linkage disequilibrium with the Gln41Leu polymorphism, suggesting that these novel GRK5 SNPs might function independently from the Gln41Leu polymorphism. Since more than 90% of GPCRs are expressed in the brain, the identification of functional GRK5 SNPs may provide novel directions for future studies of variation in antidepressant response.

In the present study, the selection of the SNPs for functional assessment was limited to those identified during our GWA analyses. While other experimental approaches are available to assess the functional consequences of genetic variants, the application of these two commonly used in vitro functional assays has highlighted nine candidate SNPs that may be worthy of further mechanistic pursuit using alternative methods. While many GWA studies have been performed with psychiatric phenotypes, few have identified genomic loci that were replicated and could successfully be used as robust “biomarkers” in clinical psychiatric practice. A lack of reliable model systems for functional genomic studies of the biological mechanism underlying the association is among the factors that have prevented the translation of genomic research to psychiatric practices. The use of pluripotent stem (iPS) cells represents a novel and promising tool for functional validation and mechanistic studies of genomic loci identified through GWA studies.[43]

Limitations of our study also include the fact that detailed information on certain clinical factors, such as comorbid psychiatric diagnoses, was not available. In addition, the influence of potentially important covariates such as drug dose and blood levels has not yet been fully explored. However, important associations between genetic variants and clinical outcomes can be missed by adjusting for covariates such as blood drug levels that serve as intermediate factors. Subsequent analyses will be focused on blood drug levels and their association with genetic variation and treatment outcomes. Our top findings were not replicated by an analysis of samples from the STAR*D study and we also failed to replicate the best association findings from the STAR*D GWA analyses. Similarly, the PGRN-AMPS analyses did not replicate the rs1126757 in the IL11 gene that was reported to have some association with escitalopram response in the GENDEP project.[9] Finally, the association of the rs6989467 SNP reported for the MARS GWA study[8] was not replicated in the PGRN-AMPs analyses. Several important differences between the studies may have contributed to the lack of replication, including differences in baseline clinical characteristics of participants in the two studies.

GWA studies have been shown to be a powerful approach for the identification of novel genomic markers for disease risk and pharmacogenomic phenotypes.[44-46] However, GWA studies have not been as informative when applied to psychiatric phenotypes of disease risk or drug response. There have now been four GWA studies conducted for antidepressant outcomes in MDD patients. None of them has demonstrated genome-wide significant associations nor have these studies replicated each other. The PGRN-AMPS is unique in that it is the only analysis that used blood drug assays to evaluate treatment adherence and is the only study that analyzed the functional significance of identified SNPs of interest. Given the underlying phenotypic heterogeneity (which is typical of studies of psychiatric illnesses) it is particularly challenging to identify homogenous treatment samples. A possible solution is to subclassify patients using biomarkers such as pharmacometabolomic characteristics to identify more biologically similar patients.[12]