|Home | About | Journals | Submit | Contact Us | Français|
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 (MDD). While no SNP associations reached the genome-wide level of significance, 14 SNPs of interest were identified for functional analysis. The rs11144870 SNP in riboflavin kinase (RFK) gene on chromosome 9 was associated with eight week treatment response (OR = 0.42, p = 1.04×10−6). The rs915120 SNP in the G protein-coupled receptor kinase 5 (GRK5) gene on chromosome 10 was associated with eight week remission (OR = 0.50, p = 1.15×10−5). 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 analytic strategy is to provide insights into the potential relevance of biologically plausible observed associations.
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.
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.
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.
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.
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).
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.
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 study7 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.
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.
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.
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.
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).
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.
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).
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 analyses7 were also evaluated for possible association with treatment response in the PGRN-AMPS samples. None of these associations were statistically significant (Supplemental Table S8).
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 treatment36 and GRK5 has been shown to regulate GPCR receptors such as the β1-adrenergic receptor37 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 study8 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
This research was supported, in part, by NIH grants RO1 GM28157, U19 GM61388 (The Pharmacogenomics Research Network), P20 1P20AA017830-01 (The Mayo Clinic Center for Individualized Treatment of Alcohol Dependence), and a PhRMA Foundation Center of Excellence in Clinical Pharmacology Award. Dr. Yuan Ji’s work was supported by a KL2 Mentored Career Development Award (NCRR Grant KL2 RR024151) and a Gerstner Family Mayo Career Development Award in Individualized Medicine. This publication was supported by NIH/NCRR/NCATS CTSA Grant Number UL1 RR024150. Its contents are solely the responsibility of the authors and do not necessarily represent the official views of the NIH. We thank the staff of the Mayo Clinic Rochester Department of Psychiatry and Psychology for their effort in recruiting the patients to the Mayo PGRN-AMPS, Dr. Ryan Abo for his assistance in gene annotation, and Lori Solmonson for her assistance with the preparation of the manuscript. The authors also wish to thank the STAR*D investigators for both their effort and generosity in sharing the clinical data and DNA samples.
Supplementary information is available at The Pharmacogenomics Journal’s website.
Conflict of interest
Dr. Mrazek has developed intellectual property that has been licensed by AssureRx Health, which has been subsequently incorporated into a physician decision support software product. He has also received research funding from AssureRx Health to create and maintain a bibliographic system designed to monitor the scientific literature related to studies addressing psychiatric pharmacogenomic relationships. No other author has any potential conflict of interest to report.