Introduction

The G protein-coupled mu opioid receptor (encoded by the OPRM1 gene) affects signal transduction pathways that mediate the effects of opioids and has an important role in opioid reward, tolerance, analgesia and dependence. The vulnerability to develop opioid addiction is partially inherited and individual differences in opioid effects may be attributed to genetic variations. Understanding the effects of this variability is of clinical importance.

The non-synonymous variant rs1799971 (118A>G, Asn40Asp) was shown to remove a potential N-glycosylation site in the extracellular domain, to be more potent in beta-endorphin binding and receptor activity, and to reduce receptor signaling efficacy.1 Healthy subjects with the 118G allele showed an increased basal level of cortisol and greater cortisol responses to opioid receptor blockade with naloxone in a population-specific manner.2, 3, 4 Several studies reported positive association of the 118G variant with opioid dependence and other substance dependencies in diverse populations,5, 6, 7 whereas other studies did not detect association.8, 9 The clinical relevance of this polymorphism for opioid analgesia and opioid adverse effects is still debatable.10, 11, 12 The 118G allele frequencies varied between populations, with high frequency in Asian populations (0.35–0.48), moderate frequency in European populations (0.1–0.17) and low frequency in African populations (<0.04).13

Another potential source of functional diversity is splice variants. At least half of the genome is expressed in alternatively spliced isoforms14 and gene functions may also be modified by single-nucleotide polymorphisms (SNPs) in alternatively spliced exons. A number of OPRM1 alternatively spliced variants have been reported in rodents and humans.15, 16, 17, 18, 19, 20 The receptors encoded by two human splice variants that retain different parts of intron 1 were shown to form heterodimers with the wild-type protein.19 An alternatively spliced exon with an alternative promoter, in the 5′ upstream region (∼28–30 kb), was identified in rodents and humans, along with cell-specific splice variants.15, 19, 21, 22, 23 A SNP in intron 1 showed an effect on OPRM1 gene expression in vitro and was associated with pain sensitivity.20

The proximal promoter region does not necessarily contain all elements required for tissue-specific gene expression, and regulatory elements can be located at great distances from the gene they regulate.14, 24 Only a few potential OPRM1 regulatory variants have been functionally characterized. SNP −554G>A decreased OPRM1 activity in a neuroblastoma cell line, and SNP −1320A>G showed increased OPRM1 activity.22, 25 Association studies of potential regulatory OPRM1 variants with substance dependence were limited to the proximal regulatory region. An association with substance dependence was shown with SNPs −1793T>A, −1699insT and −2044C>A.26, 27

In a recent association study of heroin addiction, we have detected association of two variants in intron 1 (rs510769 and rs3778151) in European Americans.9 Of note, the ‘addiction-focused’ array28 used in this study does not include OPRM1 SNPs upstream of exon 1. These results are supported by two studies in which intron 1 SNPs showed association with drug dependence in European Americans,29 or with positive response to heroin after first use in Chinese.30 We did not find a similar association in an African-American cohort.31 A 10K genome-wide association study of a subset of the same cohort from our laboratory identified association of a SNP rs1074287, located 11.6 kb upstream of exon 1, with heroin addiction.32

Our hypothesis is that there may be additional non-coding variants, possibly distant from OPRM1, that contribute to the individual phenotype (for example, response to opiates, vulnerability to develop opioid addiction and perception of pain). In addition, a combination of linked variants could contribute to the phenotypic difference in a different way than a single polymorphism.

Toward this goal, we performed an analysis of genetic variations and haplogroups in the OPRM1 100 kb 5′ flanking region, exon 1 and the 50 kb intron 1, and searched for regions in high linkage disequilibrium (LD) with the 118G and/or the intron 1 variants. Genotype data were obtained from the International HapMap Project and from tag SNPs in a sample of European Americans.

Materials and methods

Comparative genomic analysis

Comparative genomic analysis was obtained from the UCSC Genome Browser (http://genome.ucsc.edu/) and is based on the Vertebrate Multiz Alignment & PhastCons Conservation. Regulatory potential data were based on the ESPERR (evolutionary and sequence pattern extraction through reduced representation) that is computed from alignments of human, chimpanzee, macaque, mouse, rat, dog and cow.33

HapMap data

HapMap genotype and phased haplotype data were obtained from the genome browser of The International HapMap Project Phase 2 (http://www.hapmap.org/). The sample populations used were CEU (Utah residents with Northern and Western European ancestry from the CEPH collection), YRI (Yoruba in Ibadan, Nigeria) and combined JPN/CHB (Japanese in Tokyo and Han Chinese in Beijing). Only the founders in these sample populations were included. The phased haplotype data were generated using a maximum likelihood algorithm. The chromosomal region was set to position 154 300–154 450 kb (build 37.1).

Subjects

An exploratory representative group of 103 subjects was selected, on the basis of their 118A>G genotype, from our European American cohort, which was collected for genetic studies of heroin addiction as described.9 To reduce the potential effect of population stratification, only subjects with a European ancestry proportion of >0.7, based on ancestry informative markers analysis, were included (see below). The group can be divided as follows: (1) ‘118G’/‘118G’ (n=5); (2) ‘reference’ group (‘REF’)/‘118G’ (n=52); (3) ‘IVS1’/‘118G’ (n=17); and (4) ‘REF’/‘REF’ (n=29). As the 118G allele was not found to be associated with heroin addiction in this cohort, both cases (former heroin addicts in methadone maintenance treatment) and controls were included, except for subgroup 4, which included only control subjects. Subjects were recruited at the Rockefeller University Hospital, the Weill Medical College of Cornell University, and the Dr. Miriam and Sheldon G. Adelson Clinic for Drug Abuse Treatment and Research in Las Vegas, NV, USA. The Institutional Review Boards of the Rockefeller University Hospital and Cornell University approved the study for the three institutions. All subjects signed informed consent for genetic studies.

Ancestry informative markers

A total of 174 ancestry informative markers with adequate quality were used to calculate the proportion of European ancestry using STRUCTURE,34 with the CEPH diversity panel of 1051 individuals as a reference.9, 28, 35

LD structure, haplotypes and multilocus genotype patterns

The pattern of pairwise LD between the three SNPs (no. 5, 6 and 8) in the original sample9 was measured by D′ and r2 metrics. Haplotypes were reconstructed using the accelerated expectation-maximization algorithm implemented in Haploview version 4.2.36 Multi-locus genotype patterns were generated manually.

SNP selection and genotyping

Three SNPs (rs1799971 (no. 5), rs510769 (no. 6) and rs3778151 (no. 8)) were genotyped using the Illumina custom array.9 Additional five tag SNPs were selected for genotyping from the 100 kb 5′ flanking region (blocks 1: rs1551808 (no. 1), rs7758009 (no. 2) and rs7760028 (no. 3), and block 4: rs1074287 (no. 4), and rs3778146 (no. 7)) (Table 1, Figure 1 and Supplementary Table 1). The selection was based on information from HapMap, chromosomal location and frequency (MAF> 0.1 in CEU).

Table 1 LD haplotype block distribution of OPRM1 SNPs in the HapMap CEU population
Full size table
Figure 1
figure 1

(a) Schematic representation of the OPRM1 gene region and the location of the eight SNPs genotyped: rs1551808 (1), rs7758009 (2), rs7760028 (3), rs1074287 (4), rs1799971 (118A>G, 5), rs510769 (6), rs3778146 (7), and rs3778151 (8). Black boxes indicate exons, open boxes indicate untranslated regions (UTR) and dotted boxes indicate alternatively spliced exons. (b) A customized view from the UCSC genome browser that covers 100 kb of the 5′ region, exon 1 and 50 Kb of intron 1, showing pair-wise alignment of each of four species with the human genome. The tracks from top to bottom show the positions of segments aligned with the indicated comparison species. Green bars indicate conserved regions.

Full size image

DNA was extracted from blood and genotypes were determined by TaqMan technology. The TaqMan Pre-designed SNP Genotyping Assay catalog numbers (Applied Biosystems (ABI), Foster City, CA, USA) are listed in Supplementary Table 1. PCR was performed in duplicate on a GeneAmp PCR 9700 using TaqMan universal PCR master mix with AmpErase uracil-N-glycosylate (ABI) according to the manufacturer's protocol. Briefly, TaqMan assay mix (40 × ) and universal PCR master mix (2 × ) were mixed and the volume was adjusted to 4 μl in a 384-well optical reaction plate; 10 ng of genomic DNA (1 μl) was added. PCR amplification lasted for 2 min at 50°C, for 10 min at 95°C, followed by 40 cycles of 15 s at 92°C and 1 min at 60°C. Genotype analysis was performed on an ABI Prism 7900 sequence detection system using SDS 2.1 software (ABI).

For verification of the TaqMan assay, PCR amplification was performed using Platinum PCR Supermix (Invitrogen, Carlsbad, CA, USA). Primers were designed using software Primer3.37 PCR amplification consisted of 2 min at 94°C; eight ‘touch-down’ cycles of 30 s at 94°C, 30 s at 63–56°C and 30 s at 72°C; 30 cycles of 30 s at 94°C, 30 s at 56°C and 30 s at 72°C; and a final step of 7 min at 72°C. PCR products were purified and sequenced on an ABI Prism 3700 capillary sequencer (ABI). Electropherograms were scored using the Sequencer 4.5 software (Gene Codes Corporation, Ann Arbor, MI, USA).

Results

The genomic organization of the OPRM1 gene region is shown in Figure 1a. OPRM1 is mapped to chromosome 6q25.2 at position 154.4–154.7 Mb (NCBI build 37.1) and is flanked at the 5′ end by three non-functional pseudogenes. The next gene upstream to OPRM1, the regulator of G-protein signaling 17 gene (RGS17), is located at a 1 Mbp distance at position 153.4. At the 3′ end, OPRM1 is overlapped by the phosphoinositide-binding protein PIP3-E gene (IPCEF1), in the opposite orientation.

In this study, we focused on the region spanning the ∼100 kb 5′ flanking region, exon 1 and the 50 kb intron 1. Sequence comparison of this region between humans, non-human primates and rodents is shown in Figure 1b. Analysis of the LD matrix of the selected region in the HapMap CEU sample revealed that 108 SNPs are located in four major blocks of high LD (Figure 2a, Table 1). The first block spans ∼10 kb with 28 SNPs, the second block spans ∼23 kb with 22 SNPs, the third block spans ∼20 kb with 15 SNPs and the fourth block spans intron 1 (∼44 kb) with 42 SNPs. SNPs analyzed in related association studies (except for 118A>G, which was studied in many studies) are indicated by their references (also see Discussion).

Figure 2
figure 2

(a) Plot of pair-wise linkage disequilibrium (LD) of the 150 kb OPRM1 genomic region, in the HapMap CEU sample. The plot is not drawn to scale and the arrows show the correct position of the SNPs. The color scheme indicates the magnitude of r2. Black blocks indicate r2=1 and white blocks indicate r2=0. (b) The pattern and frequency of the major haplogroups in this population. The height of the boxes is proportional to the frequency in the population. Blue boxes indicate variant alleles, yellow boxes indicate reference alleles and gray boxes indicate all the other haplogroups.

Full size image

The HapMap CEU haplotype data revealed four major haplogroups (Figure 2b): (1) the ‘REF’ (∼36%) with no variants (including the reference 118A allele and the reference intron 1 sequence); (2) the block 2 group (∼22%); (3) the 118G group (∼15%) with additional variants in block 1; and (4) the intron 1 variants of block 4 (‘INT1’, ∼11%). The minor haplogroups are block 3 (4%) and 12% of the haplotypes that could not be assigned to one of the major haplogroups (grouped as ‘Others’ in Figure 2b). Of the 59 HapMap CEU subjects, 10 (17%) were homozygous for the ‘REF’ haplogroup, two (3%) were homozygous for the 118G haplogroup and five (8%) were homozygous for the block 2 haplogroup. No individual was homozygous for the ‘INT1’ haplogroup and one was heterozygous for the ‘118G’/‘INT1’ groups. All the rest were heterozygotes with different allelic combinations.

The two SNPs (rs510769 (no. 6) and rs3778151 (no. 8)) identified in our recent study to be associated with heroin addiction9 are part of the 42-SNP block 4 that spans intron 1. The average minor allele frequency (MAF) of the SNPs of block 4 is 0.15 in European Americans. Interestingly, four SNPs (rs9478498, rs9478499, rs1074287 (no. 4) and rs6936615), located at the −20 kb 5′ region (block 3 region) are in high LD with the SNPs of block 4. Four SNPs (rs7748401, rs10457090, rs3778152 and rs563649) located at intron 1 are not part of block 4 and are in high LD with the SNPs of block 3. Haplogroup ‘INT1’ includes the reference 118A allele (Table 1). Seven SNPs (rs9384167, rs12174208, rs11966947, rs7760028 (no. 3), rs12527423, rs17084868 and rs17084870) are in high LD with SNPs in block 1, but are physically located in or adjacent to the region of block 2 (Table 1, Figure 2a). The region that spans block 1 (−100 kb to −90 kb) shows several regions with a regulatory potential (Supplementary Figure 1).

We searched for similar haplogroups in the other HapMap populations (the combined JPT/CHB and YRI). The ‘INT1’ haplogroup that includes the intron 1 variants of block 4 accounts for ∼4% in these populations. In contrast, the ‘118G’ haplogroup has a very different pattern. In concordance with the 118G allele frequency, the 118G haplogroup that includes the 5′ variants of block 1, is very frequent in the Asian population and can be divided into two subgroups: the major one (28%) includes the variants in block 1, similar to European Americans, and the second (14%) has other combinations (data not shown). The 118G allele was not found in the African sample but a similar haplogroup with block 1 variants (without the 118G variant allele) was represented by 7.5% of the haplotypes. Interestingly, there was a significant representation (∼30%) of the block 3 haplogroup, with variants at the proximal 5′ region, in the African sample, in contrast to the low frequency in Asian and Caucasian samples.

Three major haplogroups were identified in our previous study of European Americans (350 cases and 184 controls) on the basis of OPRM1 exon 1 SNP rs1799971 (118A>G, no. 5), and intron 1 SNPs rs510769 (no. 6) and rs3778151 (no. 8).9 The major haplogroups identified are: (1) ‘REF’ (65, 56%); (2) ‘118G’ (13, 11%); and (3) ‘INT1’ (13, 20%) in controls and cases. The higher frequency of the ‘REF’ group may be explained by the fact that the haplogroup of ‘block 2’ cannot be distinguished from the reference group by the three SNPs genotyped. These results demonstrate the higher frequency of the ‘INT1’ haplogroup in cases, which was the basis for the association of this haplogroup with heroin addiction. The allele frequencies in controls are compatible with the HapMap data from Caucasians.

To further explore the finding that the 118G haplogroup includes variant alleles in the distant 5′ region, a representative group of 103 subjects with and without the 118G allele was selected from the original cohort. This exploratory sample does not represent the general population and was selected to enable unambiguous phase determination. To reduce population admixture, only subjects with a European ancestry proportion of >0.7, based on ancestry informative markers analysis, were included. The exploratory representative sample is composed of the following subgroups (Table 2): 1) ‘118G’/‘118G’ (homozygotes for the 118G allele, n=5); 2) ‘REF’/‘118G’ (heterozygotes for 118G allele that are non-carriers of the IVS 1 variants, n=52); 3) ‘IVS1’/‘118G’ (double heterozygotes for 118A>G and the two IVS1 SNPs, n=17); and 4) ‘REF’/‘REF’ (represents the reference genotype group, n=29, all controls).

Table 2 Multi-locus genotype patterns (MLGPs) in a representative sample of the 118A>G genotype subgroups
Full size table

All subjects were genotyped for an additional five SNPs (rs1551808 (no. 1), rs7758009 (no. 2), rs7760028 (no. 3), rs1074287 (no. 4) and rs3778146 (no. 7)) (Table 1, Figure 1). The SNPs were selected from the 100 kb 5′ flanking region, and intron 1 on the basis of LD information from HapMap (CEU blocks 1 and 4) and frequency in CEU (MAF>0.1). As shown in Table 2, the five subjects with the ‘118G’/‘118G’ pattern (subgroup 1) carry two copies of the three variant alleles (no. 1–3) of the SNPs from the 5′ region (C-G-C, block 1), two copies of the reference A allele of the 5′ SNP rs1074287 (no. 4) and two copies for the reference allele of the three SNPs from intron 1 (C-T-T, no. 6–8). The majority of the heterozygotes for the 118G allele (‘REF’/‘118G’, subgroup 2) also carry at least one variant allele of the SNPs no. 1–3 (TC-AG-TC), and none of them carries the variant allele of SNP rs1074287 (no. 4) or of the three IVS1 SNPs (no. 6–8). The subjects in subgroup ‘IVS1’/‘118G’ (subgroup 4) all carry at least one variant allele of SNPs no. 1–4 (TC-AG-TC-AG). Although the phase in these subjects cannot be determined unambiguously, the genotypes most probably reflect a combination of two haplogroups: (1) SNP no. 4/IVS1 SNPs no. 6–8 and (2) 118G/5′ SNPs no. 1–3, based on the other subgroups in which the phase can be determined with certainty. This subgroup represents a small group of the original cohort (5% in controls and 8% in cases). The majority of the ‘REF’/‘REF’ subgroup (subgroup 3) does not carry a variant allele of any of the genotyped SNPs. This exploratory data corroborate the HapMap data that the 118G haplogroup includes variant alleles of SNPs in the distant 5′ region and is distinct from the ‘INT1’ haplogroup.

Discussion

Genetic diseases can be caused by disruption of the regulatory elements, and identification of these regulatory elements is an important and challenging problem. Tissue-specific and development-specific expression requires sophisticated expression profiles that may involve multiple enhancer elements. Long-range cis elements and intronic variants may have a role in this regulation. For example, an evolutionarily conserved enhancer was identified 86 kb upstream of the peptidylarginine deiminases gene (PADI3), the expression pattern of which was not explained by previous studies limited to the proximal promoter.38

Recently, an alternatively spliced exon with a specific promoter was identified in humans, ∼28 kb upstream of exon 1,15 suggesting an extended OPRM1 gene. It is also of note that, except for three pseudogenes, there is no gene located at 1 Mbp 5′ to OPRM1.

In our hypothesis-driven association study of heroin addiction, two OPRM1 intron 1 variants (no. 6 and 8) were found at a significantly higher frequency in subjects with heroin addiction than in healthy controls, suggesting an association with heroin addiction.9 A 10K genome-wide association study from our laboratory that used a subset of the sample population used in the hypothesis-driven study indicated association of SNP rs1074287 (no. 4) in the 5′ flanking region (−11.6 kb from exon 1) with heroin addiction.32 Analysis of the HapMap CEU LD data in this study revealed a unique haplogroup (‘INT1’) of 42 SNPs that spans intron 1 and also includes at least four variants in the 20 kb 5′ flanking region, including rs1074287. These data indicate that the findings of the two studies are compatible and refer to the same ‘IVS1’ haplogroup. On the basis of these data, the intron 1 SNPs that show association with heroin addiction may be markers in high LD with regulatory variants in the proximal 5′ region. It is also possible that the SNPs in the 5′ flanking region are markers for functional SNPs in intron 1.

Several lines of evidence suggest an additional functionality to intron 1 beyond the basic splicing of the constitutive exons. Several alternatively spliced exons were described in this intron.19, 20 The finding of a high LD between SNPs spanning the entirety of intron 1 and a few SNPs in the proximal 5′ flanking region in different populations is intriguing. This haplogroup is represented at a relatively low frequency (4%) in Asian and African populations but at a much higher frequency (15%) in the European population, suggesting the possibility of an ancient haplotype with a specific selective advantage in Europeans. It is also possible that this population's difference occurs by random, as was recently suggested in a genome-wide study of genes with highly divergent allele frequencies between populations.39

A potential functional SNP in intron 1 (rs563649) was recently indicated in association with pain perception.20 The HapMap LD data indicate that, although this SNP is not a part of the ‘INT1’ haplogroup, it is in high LD with SNPs of haplogroup 3 that includes variants in the proximal 5′ region (block 3). Notably, this haplogroup is rare in European Americans, but frequent in the African sample. It is possible that this SNP is a marker of regulatory variant/s in the proximal 5′ region.

Many studies provided evidence for the importance of the OPRM1 118A>G polymorphism in addiction, pain, HPA activation and treatment response (for example, 2-9). These studies have shown mixed results regarding the association between the 118G variant and opioid dependence, opioid analgesia and opioid adverse effects.40 Most studies have used single SNP analyses and were generally limited to the coding region or to the proximal upstream region. Several studies employed haplotype analysis but were also limited to specific regions. Hoehe et al.26 performed a comprehensive haplotype analysis of OPRM1 but the study did not include SNPs from intron 1 and was limited to the proximal 5′ region (−2.4 kb). Lue et al.27 studied five variants from the proximal 5′ flanking region (−2 kb) and two SNPs from exon 1. Zhang et al.29, 30 examined 13 SNPs spanning the coding sequence, including several SNPs from intron 1, but none in the 5′ flanking region. Xuei et al.41 genotyped 18 SNPs including one in the proximal 5′ region (−6 kb). A recent study used a haplotype-based approach to predict response to naltrexone treatment on the basis of 10 OPRM1 tag SNPs but was limited to the proximal 5′ region (−11 kb).42 Haplotype-based analysis is generally considered more powerful than single-marker analysis and may help to determine whether a specific polymorphism occurs on a common genetic background.

The mixed single SNP association results for OPRM1 SNP 118A>G may be explained by different haplotype patterns. For example, the finding that about a third of the chromosomes with the 118G allele in the HapMap Asian populations do not show high LD with block 1 SNPs suggests the existence of a second 118G haplogroup that is specific to this population. It is feasible that single SNP association studies of the 118A>G SNP in this population may result in mixed results. The data obtained in this study indicate that, in European Americans, the 118G variant allele is positioned within a specific haplogroup and is in high LD with several distant variants that may possibly have a regulatory effect. It remains to be shown whether these variants are correlated with OPRM1 expression levels. The variant 118G allele was not found in the African HapMap population, suggesting that this is a relatively new variant that arose after the out-of-Africa migration of modern humans. The high prevalence of this allele in the Asian population and the lesser prevalence in European Americans may be explained by genetic drift and/or selective advantage. Interestingly, a similar haplogroup (block 1), without the 118G allele, was found in the African HapMap sample, suggesting that the change to 118G may have occurred on this ancestral haplotype.

In summary, the identification of LD blocks and haplogroups in the OPRM1 intron 1 and the 5′ flanking region may be of importance to the field of opioid response and addiction. This study suggests an extended OPRM1 locus and possible new regulatory regions. The haplotype architecture described in this study may provide direction for future genetic studies. These findings will require confirmation by functional expression studies.