|Home | About | Journals | Submit | Contact Us | Français|
Heroin addiction is a chronic complex disease with a substantial genetic contribution. This study was designed to identify gene variants associated with heroin addiction in African Americans. The emphasis was on genes involved in reward modulation, behavioral control, cognitive function, signal transduction, and stress response. We have performed a case-control association analysis by screening with 1350 variants of 130 genes. The sample consisted of 202 former severe heroin addicts in methadone treatment and 167 healthy controls with no history of drug abuse. Single-SNP, haplotype and multi-SNP genotype pattern analyses were performed. Seventeen SNPs showed point-wise significant association with heroin addiction (nominal P < 0.01). These SNPs are from genes encoding several receptors: adrenergic (ADRA1A), arginine vasopressin (AVPR1A), cholinergic (CHRM2), dopamine (DRD1), GABA-A (GABRB3), glutamate (GRIN2A), and serotonin (HTR3A), as well as alcohol dehydrogenase (ADH7), glutamic acid decarboxylase (GAD1 and GAD2), the nucleoside transporter (SLC29A1), and diazepam binding inhibitor (DBI). The most significant result of the analyses was obtained for the GRIN2A haplotype G-A-T (rs4587976-rs1071502-rs1366076) with protective effect (P uncorrected = 9.6E-05, P corrected = 0.058). This study corroborates several reported associations with alcohol and drug addiction as well as other related disorders, and extends the list of variants that may affect the development of heroin addiction. Further studies will be necessary to replicate these associations and to elucidate the roles of these variants in drug addiction vulnerability.
Heroin addiction is a chronic, relapsing brain disease that is characterized by drug dependence, tolerance and compulsive seeking and use despite harmful consequences. This complex disorder is a worldwide major public health problem. The relatively high heritability (40-60%) of this disorder (Tsuang et al., 1996) indicates that genetic variants may play a role in vulnerability to its development. Chronic drug use alters gene expression, which activates or attenuates biochemical pathways and produces neuroadaptive changes in signal transduction functions. Identifying these variants is important for the understanding of the possible causes of this disease and to improve its diagnosis and treatment, as well as for primary prevention purposes.
Studies have suggested association of a number of gene variants with heroin addiction or combined drug addiction. These include variants in genes encoding opioid, dopamine, serotonin, and GABA receptors and also cholinergic muscarinic 2, ACTH, and cannabinoid 1 receptors. In addition, association was suggested for catechol-O-methyltransferase (COMT), period circadian protein 3, proenkephalin, proopiomelanocortin, tryptophan hydroxylase 2, and brain-derived neurotrophic factor (Kreek et al., 2005, Cheng et al., 2005, Loh et al., 2007, Nielsen et al., 2008, Proudnikov et al., 2008, Proudnikov et al., 2006, Szilagyi et al., 2005, Xu et al., 2004, Zou et al., 2007, Zuo et al., 2007)}. We have recently conducted a candidate gene association study in European Americans (EA), in which we detected associations with SNPs in OPRM1, OPRK1 OPRD1, GAL, HTR3B and CSNK1E, with severe heroin addiction (Levran et al., 2008)
Specific disease-associated genetic factors may vary in frequency among populations because of random drift or natural selection (Lohmueller et al., 2006). The degree and pattern of linkage disequilibrium (LD) also differ among populations. In addition, some functional variants may have different effects in different ethnic groups (Tate & Goldstein, 2004).
We performed a hypothesis -driven case-control association study, using an “addiction-focused” SNP array (Hodgkinson et al., 2008) in African Americans (AA) that extends our original study in EA (Levran et al., 2008). The two populations were separated for analysis based on significant differences in allele frequencies and the vulnerabilities of these studies to population stratification. The cases were selected from the extreme margin of the specific phenotype range (e.g., severe heroin addicts in methadone maintenance treatment) and the controls were selected by detailed personal interview and stringent criteria.
Three hundred and sixty nine subjects participated in the study (202 cases and 167 controls). There were 54% percent females in the control group and 38% in the case group. The mean age at recruitment was 31.9 (±11.7) years in the control group and 47.0 (±8.8) years in the case group. The patients (cases) were former severe heroin addicts from methadone maintenance treatment programs. Subjects were recruited at the Rockefeller University Hospital, the Manhattan Campus of the VA NY Harbor Health Care System, Weill Medical College of Cornell University, and the Dr. Miriam and Sheldon G. Adelson Clinic for Drug Abuse Treatment and Research, Las Vegas. Ascertainment was made by personal interview, using several instruments: the Addiction Severity Index (ASI) (Mclellan et al., 1992), the Kreek-McHugh-Schluger-Kellogg scale (KMSK) (Kellogg et al., 2003), and the Diagnostic and Statistical Manual of Mental Disorders, 4th Edition (DSM-IV). All subjects had a history of using heroin multiple times per day, every day for at least one year. The 167 control subjects were recruited at the Rockefeller University Hospital. Each of the following was used as an exclusion criterion from this category: a) At least one instance of drinking to intoxication, or any illicit drug use, in the previous 30 days. b) A past history of alcohol drinking to intoxication, or illicit drug use, more than twice a week, for more than six consecutive months. c) Cannabis use for more than 12 days in the prior 30 days or past use for more than twice a week for more than 4 years. All subjects completed a family history questionnaire and were self-identified as AA for three generations. Participants were excluded from the study if they had a relative in the study, or had a mixed ancestry. The Institutional Review Boards of Rockefeller University, VA NY and Cornell University approved the study. Rockefeller University IRB also reviews the Adeslon Clinic, Las Vegas.. All subjects signed informed consent for genetic studies.
Blood samples were taken and DNA was extracted using the standard salting-out method. DNA was quantified using PicoGreen (Invitrogen, Carlsbad, CA). 700 ng DNA (45 μL) was ethanol-precipitated for a second time and re-suspended in 7 μl Tris-EDTA (10 mM Tris-HCl, 1 mM EDTA, pH 8.0) as described (Levran et al., 2008). 5 μl (250-500 ng) were used for genotyping.
Genotyping was performed on a 1,536-plex GoldenGate Custom Panel (GS0007064-OPA, Illumina, San Diego, CA). 1350 SNPs were selected from 130 candidate genes (Hodgkinson et al., 2008). In addition, the array included 186 ancestry informative markers (AIMs) that were selected based on allele frequencies in the European, African and Chinese population of the HapMap project (Enoch et al., 2006, Hodgkinson et al., 2008). Genotyping was performed at the Rockefeller University Genomics Resource Center according to the manufacturer's protocol (Illumina). Analysis was performed using BeadStudio genotyping software (Illumina). Genotype data was filtered based on SNP call rates (>99.5%), MAF>0.01, cluster separation score, and deviation from Hardy-Weinberg Equilibrium (HWE). Random samples (~10%) were genotyped in duplicate and two identical samples were genotyped on each of the arrays for reproducibility control purposes.
Ancestry informative markers were employed to test for population stratification using the STRUCTURE v2.2 software (Pritchard et al., 2000). This software places individuals in K clusters of potential putative populations separated by marker frequencies. A no-admixture ancestry model with K set to 1 was initially used to infer the Dirichlet parameter lambda of the distribution of allele frequencies. The value of lambda thus obtained was later used under an admixture model with allele frequencies correlated using a burn-in length of 50,000 followed by 50,000 Markov Chain Monte Carlo iterations which was sufficient to yield stable results.
STRUCTURE was also used to compare the AA sample in this study with the EA sample analyzed in a recent study (Levran et al., 2008). Although two clear clusters were observed, K was set to 3 for better display of the results in the yielded triangle plot. The coordinates of any given point in this triangle are given by the relative distance to each edge (putative population). Individuals who are in one of the corners are therefore assigned completely to a particular population.
The genomic control method (Devlin et al., 2001) was used to estimate the degree of subdivision in our study sample. In this method, markers are used to estimate the background association which can then be used to infer the variance inflation factor, lambda, introduced when there is subdivision.
The frequency of marker alleles and genotypes in cases and controls were compared using Fisher's exact test, as implemented in R v2.7.0 (http://www.R-project.org). This test was preferred over a parametric test since there were numerous instances of small cell counts in the tables. This test was also used to estimate allelic odds ratios (OR) with their corresponding 95% confidence intervals. The most significant results were also evaluated using the chi-square test and no substantial differences were found. An uncorrected P-value of 0.01 was chosen as the threshold for point-wise significance for the association tests.
Marker genotype quality control was performed in the control sample by using the exact test of deviation from the Hardy-Weinberg (DHW) proportions as implemented in the R package Genetics Base v1.6.0 (http://rgenetics.org/). This test estimates allelic correlations occurring at the locus level which may be caused by genotyping errors. The DHW test is not powerful enough to detect deviations caused by small error rates but can be used to detect gross genotyping errors which can cause large DHW. For this reason, a conservative p-value of 0.05/1114 = 0.000045 was used as a cut-off for declaring significant departures from HW proportions.
To adjust for possible age and gender effects, a logistic regression analysis was run on the set of markers using PLINK v1.05 (http://pngu.mgh.harvard.edu/purcell/plink/)
The pattern of pair-wise LD between SNPs was measured by the standardized disequilibrium value (D') as implemented in Haploview v4.1 (Barrett et al., 2005). In order to capture the blocks that consist of the specific SNPs indicated in the single SNP association tests, LD haplotype blocks were identified using the Gabriel's method (Gabriel et al., 2002). When no block was identified by this method, the alternative solid spine of LD definition (Barrett et al., 2005) was employed.
Haplotypes were reconstructed using the accelerated expectation-maximization algorithm implemented in Haploview for genes that include at least one marker whose P-value was less than 0.01 in the single-locus association test. Haplotype association analysis was performed for each of the reconstructed haplotypes (against the rest) using the Fisher's exact test. Multi-locus genotype patterns (MLGP) of the markers in the blocks identified by Haploview were determined. For a set of a given number m SNPs, n observed m locus genotype pattern frequencies were stored in an n x 2 contingency table where the two columns corresponded to cases and controls. A Fisher's exact test for differences in MLGP frequencies between the two groups was performed. In addition to single-locus effects, MLGP analysis also evaluates epistatic interaction influencing the outcome. Since we are not interested in the effects of rare genotype patterns, in this analysis, pattern frequencies below 5% in cases and controls each were pooled into a single “rare” class.
Results were adjusted for multiple testing by using the R package QVALUE v1.1 (Storey, 2003, Storey & Tibshirani, 2003). Rather than dealing with the probability of rejecting a true null hypothesis (the family-wise error rate), the approach implemented in QVALUE estimates the expected proportion of false positives among all rejected hypotheses (the false discovery rate, FDR) (Benjamini et al., 2001). QVALUE takes a given set of P values and, for each test, estimates the minimum FDR that is incurred when calling that particular test significant (the q-value of the test). The q-value measures the significance of each of a family of tests performed simultaneously and holds under different forms of dependence. The smallest nominal P value of all tests performed (Pmin) was observed from the haplotype association tests. The QVALUE program was run on the list of P values created by adding Pmin to the set of P values obtained from the single-locus tests. The result is the estimated experiment-wise significance of Pmin. An FDR of 0.05 was used as the significance level.
Genotype and allele frequencies of 1194 SNPs from 130 candidate genes (Hodgkinson et al., 2008) were analyzed for association with heroin addiction in 369 AA samples (202 cases and 167 controls). One hundred fifty six SNPs were excluded from analysis because of inadequate quality (cluster separation score < 0.4), low variability (MAF < 0.01) (Table S1), or deviation from Hardy-Weinberg equilibrium in the control sample. Genotyping reproducibility was 99.9% and the mean marker call rate after exclusion was 99.8%.
Results of the most significant single SNP association tests are shown in Tables 1 and and2.2. Indication of association with heroin addiction was observed with seventeen SNPs (uncorrected P < 0.01) in the following genes: glutamate receptor, ionotropic, N-methyl D-aspartate (NMDA) subtype 2A, solute carrier family 29 (nucleoside transporters) member 1, dopamine receptor D1, alcohol dehydrogenase isozyme 7, 5-hydroxytryptamine (5-HT, serotonin) receptor, subtype 3A, glutamate decarboxylase isoforms 1 and 2, GABA-A receptor, subunit beta 3, diazepam binding inhibitor, cholinergic receptor, muscarinic 2, adrenergic receptor alpha-1A, and arginine vasopressin receptor subtype 1A. Listed in Table S2 are the alleles and genotype frequencies in cases and controls. Odds ratios were calculated for the minor allele and indicate a small effect (OR for risk effect range 1.54-1.94 and OR for a protective effect range 0.16-0.66 for the allelic test, Table 2). None of the tests were significant after correction for multiple testing. No significant effect was found for either age or gender.
Haplotypes were inferred from LD blocks, which include at least one SNP from the list of top signals on the single SNP analyses. Nominally significant (P < 0.008) association tests for haplotypes are listed in Table 3 and the relevant LD maps are shown in Fig. 1. Association was suggested for haplotypes of GRIN2A, GAD1, HTR3A, DRD1, ADH7 and DBI.
Four GRIN2A variants (rs1070487, rs6497730, rs4587976, and rs1650420), all located at a 32 kb area of intron 3 (5' to the translation site at exon 3), accounted for some of the strongest signals in the association test (P = 0.0006-0.0039, Tables 1 and and2).2). Two additional SNPS from the same block (rs1071502 and rs1366076) gave nominal significant P values for association (P < 0.05) but did not pass the threshold value. Eleven additional GRIN2A SNPs gave negative results. The LD map and haplotype block structure of this region are shown on Fig. 1e; SNPs rs1070487 and rs6497730, are in strong LD (D' = 0.88). SNP rs4587976 forms a 7 kb block with SNPs rs1071502 and rs1366076. SNP rs1650420 is in complete LD (D' = 1) with rs1366076 but is not part of a block, under this block definitions. Haplotype analysis of block 1 (rs4587976-rs1071502-rs1366076) revealed significant association of haplotypes GAT (protective) and CAT (risk), (uncorrected P =9.6E-05 and 0.0036, respectively, Table 3). The association test of the GAT haplotype was close to significance after correction for multiple testing (P = 0.058). The contributing SNP to this effect is rs4587976 (C as a risk allele, G as a protective allele) in concordance with the single SNP analysis. Multi-locus genotype pattern analysis of this block revealed a significant difference between cases and controls with a similar effect, uncorrected P = 0.0005, data not shown).
These four SNPs (rs1650420, rs6497730, rs4587976 and rs1070487) are common in both AA and EA (MAF > 0.33, Table 4), but the minor allele frequencies differ between these ethnic groups. The four SNPs are more common in EA and the difference in allele frequency of rs1650420 is significant after correction for multiple testing (P = 3.5E-06, Table 4). The minor alleles of SNPs rs1650420, rs1070487 and rs6497730 in AA are the major alleles in EA (Table 4).
STRUCTURE analysis using 174 AIMs with adequate quality excluded population stratification between cases and controls in this study (Fig. 2). It also shows clear distinction between the AA sample and the EA sample in our recent (Levran et al., 2008) (Fig. 3). The average extent of admixture (European ancestry proportion) in this AA sample population is 11%. Five AA individuals showed higher than expected EA admixture (0.64-0.84) and the analyses were repeated without them with no significant change in the results.
The allele frequencies (in controls) of the 28 SNPs that gave the most significant association with heroin addiction in the AA population (this study) and the EA population (Levran et al. 2008) are listed in Table 4. Significant differences, after correction for multiple testing (uncorrected P < 4.5E-05), between the two populations were observed for 12 SNPs. These include five SNPs in the EA study and seven SNPs in this study. Six additional SNPs, indicated in the current study, have point-wise significant difference in allele frequencies between the two populations (4.5E-05 < P < 2.3E-02). Four SNPs in this study and one SNP in the EA study have non significant differences in allele frequencies (Table 4). The allele frequencies of reference populations from HapMap (Yoruba in Ibadan, Nigeria (YRI) and Utah residents with Northern and Western European ancestry (CEU), (The International HapMap Consortium, 2003)) are listed in Table 4 for comparison. There is high similarity between the data from our populations and the HapMap data with an expected difference between the admixed AA population and the YRI population.
The goal in this study was to identify gene variations which contribute to vulnerability to develop heroin addiction in African Americans. The most significant results that remained close to significance after correction for multiple testing were obtained for the GRIN2A haplotype. Glutamatergic neurotransmission is the major excitatory system in human brain, and genes encoding glutamate receptors are candidate genes for neuropsychiatric disorders. The action of glutamate is mediated by a few receptor families including the NMDA receptor family that is involved in multiple cognitive processes including memory and learning (Paoletti & Neyton, 2007). NMDA receptor Ca2+ permeable channels are heteromers composed of subunit GRIN1 and one or more of the 4 subunits: GRIN2A-D. Expression of the GRIN2A subunit begins around puberty, and a reduced GRIN2A expression was suggested as a risk factor for schizophrenia (Watanabe et al., 1993). Repeated 3,4 methylenedioxymethamphetamine (Ecstasy) administration in rats induces neuroadaptive changes in expressions of glutamatergic NMDA subunits, in regions of the brain regulating reward-related associative learning, cognition, and memory (Kindlundh-Hogberg et al., 2008). Grin2a (Nr2a) knockout mice exhibit retarded discrimination learning (Brigman et al., 2008). Associations were reported for GRIN2A polymorphisms with attention-deficit/hyperactivity disorder (ADHD), autism, and worse chronic outcome in schizophrenia (Adams et al., 2004, Barnby et al., 2005, Itokawa et al., 2003, Tang et al., 2006, Turic et al., 2004). The biological functions of the GRIN2A SNPs indicated in this study are unknown but the finding is supported by associations with multiple SNPs at the same gene, and by the haplotype analysis that was close to significance after correction for multiple testing.
Although all the associations indicated in this study were not significant after correction for multiple testing and could have occurred by chance, some of them corroborate reported associations with alcohol and drug addiction as well as other related disorders. The DRD1 SNP rs5326, from this study is in strong LD (D' = 0.93) with the functional 3' UTR SNP rs686 that was associated with nicotine dependence in AA sample (Huang et al., 2008). The dopamine receptor D1 mediates many of the reinforcing and dependence producing properties of drugs (Koob, 1992), and DRD1 SNPs were associated with compulsive addictive behaviors, sensation seeking in alcohol dependents, and smoking (Comings et al., 1997, Kim et al., 2007, Limosin et al., 2003). The association of ADH7 SNP rs971074 in this study is in concordance with the results of Luo et al. for substance dependence and personality traits (Luo et al., 2006, 2007a, b) and may be related to the role of ADH enzymes in the development and maintenance of the dopaminergic system (Chambon, 1993). The 5' UTR AVPR1A SNP rs3759292, indicated in this study, is located 179 bp 5' to the microsatellite polymorphisms rs9325177 and rs11283312 that were linked to autism, personality traits and altruism (Knafo et al., 2008, Meyer-Lindenberg et al., 2008, Wassink et al., 2004). DBI SNP rs2289948, indicated in this study, is in complete LD with SNP rs8192506 that was associated with alcoholism in Japanese (Waga et al., 2007) and with anxiety disorders with panic attacks in EA (Thoeringer et al., 2007). The second DBI SNP rs12613135, indicated in this study, is located in the same haplotype block with SNPs rs2289948 and rs8192506. The diazepam binding inhibitor may play a role in both physical and psychological dependence on opiates as well as anxiety disorders, as indicated by rat studies (Liu et al., 2005). GAD2 SNP rs2058725, from this study, is in the same haplotype block and in complete LD with SNP rs701492 which was associated with alcoholism in Han Taiwanese (Loh El et al., 2006). Other GAD2 SNPs showed association with alcoholism in EA (Lappalainen et al., 2007). The CHRM2 5' UTR SNP rs2350780, indicated in this study, is in strong LD with SNP rs1455858 that was associated with alcohol dependence, drug dependence and affective disorders in EA and AA (Luo et al., 2005). This gene may be involved in memory, higher cognition and behavioral effects of drugs of abuse.
The candidate gene approach is limited in scope, but this study suggests that testing biological candidate genes in a well-characterized population with a well-defined phenotype may facilitate the identification of genetic variants that are associated with vulnerability to develop heroin addiction and are specific to an ethnic group. This approach is more economical than whole-genome screening and may increase the power of the study due to the smaller number of tests performed. We cannot determine if the negative findings in this study indicate true negative association, because the SNP coverage of some of the genes may be limited and due to the limited power of the study.
Population stratification and admixture are critical problems in case-control association studies. Spurious associations can arise when mixed populations, with different allele frequencies, are used. The AA population is a recently admixed population with estimated mean European ancestry proportion of 17-20% [e.g. (Guthery et al., 2007, Reiner et al., 2005). In contrast, in this study population, there was an average of 11% European ancestry admixture. To minimize the effect of admixture we applied two independent methods and excluded population stratification between cases and controls. The AIMs data also validated the accuracy of the three generations ethnic self-identification in our dataset. The few exceptions (~1%) may be due to the individuals' inaccurate knowledge of their own family histories.
It is intriguing that the SNPs with the lowest P values for association with heroin addiction in the AA populations, in this study, and the EA population in our recent study (Levran et al., 2008), are so distinct. One possible explanation for these findings is the differences in allele frequencies between the two populations that are observed for 18 SNPs out of the total of 26 SNPs indicated in both studies and are supported by the HapMap data. The results for the rest of the SNPs, if their association with heroin addiction is validated, can not be explained by allele frequency differences and may be influenced by other factors such as LD pattern or modifier genes that may be different in these populations.
It was documented that the risk to develop heroin addiction is similar among all ethnic groups. The possibility that the specific gene variants contributing to the vulnerability to develop heroin addiction differ among ethnic groups, is intriguing. This phenomenon was found for Crohn's disease (where the susceptibility CARD15 variants in EA were not present in Japanese), AIDS (where the protective CCR5 variant is common in EA and absent in other ethnic groups), and Alzheimer's disease (where the APOE E4 variant have differential effect between ethnic groups) (Tate & Goldstein, 2004). It points to a public health challenge to provide population-specific, and in the future individual-specific, diagnosis and treatment.
We thank all the clinical staff that enrolled and assessed subjects for this study, including Elizabeth Ducat, Brenda Ray, Dorothy Melia and Lisa Borg. We are grateful to D. Goldman and his group from the NIH/NIAAA, for the design of the micro-array; Connie Zhao and Bin Zhang from the Rockefeller Genomic Resource Center, for their excellent assistance in genotyping. We thank Ann Ho, David Nielsen, Vadim Yuferov and Susan Russo for discussion, comments and/or proofreading of the manuscript. We would like to express our profound gratitude to the late K. Steven LaForge for setting the foundation for this study.