Next-generation sequencing of human opioid receptor genes based on a custom AmpliSeq™ library and ion torrent personal genome machine.

BACKGROUND: The opioid system is involved in the control of pain, reward, addictive behaviors and vegetative effects. Opioids exert their pharmacological actions through the agonistic binding at opioid receptors and variation in the coding genes has been found to modulate opioid receptor expression or signaling. However, a limited selection of functional opioid receptor variants is perceived as insufficient in providing a genetic diagnosis of clinical phenotypes and therefore, unrestricted access to opioid receptor genetics is required.
METHODS: Next-generation sequencing (NGS) workflow was based on a custom AmpliSeq™ panel and designed for sequencing of human genes related to the opioid receptor group (OPRM1, OPRD1, OPRK1, SIGMA1, OPRL1) on an Ion PGM™ Sequencer. A cohort of 79 previously studied chronic pain patients was screened to evaluate and validate the detection of exomic sequences of the coding genes with 25 base pair exon padding. In-silico analysis was performed using SNP and Variation Suite® software.
RESULTS: The amplicons covered approximately 90% of the target sequence. A median of 2.54×106 reads per run was obtained generating a total of 35,447 nucleotide reads from each DNA sample. This identified approximately 100 chromosome loci where nucleotides deviated from the reference sequence GRCh37 hg19, including functional variants such as the OPRM1 rs1799971 SNP (118 A>G) as the most scientifically regarded variant or rs563649 SNP coding for μ-opioid receptor splice variants. Correspondence between NGS and Sanger derived nucleotide sequences was 100%.
CONCLUSION: Results suggested that the NGS approach based on AmpliSeq™ libraries and Ion PGM sequencing is a highly efficient mutation detection method. It is suitable for large-scale sequencing of opioid receptor genes. The method includes the variants studied so far for functional associations and adds a large amount of genetic information as a basis for complete analysis of human opioid receptor genetics and its functional consequences. Copyright Â

Introduction

Opioid signaling is triggered by endogenous opioid peptides comprising endorphins, enkephalins and dynorphins with preferential selectivities at the μ-, κ- and δ-opioid receptors, respectively [1]. The existence of opioid binding sites in the brain was established in 1973 [2] and subsequently located in various different physiological systems. In peripheral [3], [4] and central [5], [6], [7] parts of the nociceptive system, activation confers an effective nocifensive mechanism. Activation at the Area postrema causes respiratory depression [8], activation in brain areas belonging to the reward system [9] is involved in addiction [10], [11], [12], activation in the gastrointestinal tract causes constipation [13], and activation at the hypothalamic–pituitary–adrenal (HPA) axis is involved in stress responses [14].

Genetic variation of human opioid receptors is an active research topic that has identified several possible modulators of the individual response to endogenous and exogenous opioids [15]. Variants were associated with the modulation of opioid receptor expression [16], [17] or signaling [17], [18], [19] up to an almost complete functional loss [18], [20]. However, the collection of certain opioid receptor variants has provided only modest explanations of clinical phenotypes [21], [22]. Therefore, a limited selection of functional opioid receptor variants for which specific genetic assays had been established [23] is perceived as insufficient in providing a genetic diagnosis of the clinical phenotype [24].

Since the early 70th of the last century when the first processes of Sanger sequencing were introduced [25], [26], the access to the whole genomic information has been widely accepted as a valuable method in clinical research [27]. With the recent broader availability of next generation sequencing (NGS) [28] the limitation to known functional variants has therefore fallen in favor of unrestricted access to opioid receptor genetics, which has already been shown to provide a working genetic marker for opioid-related clinical phenotypes [29]. In this report, the evaluation of a new NGS method based on a custom AmpliSeq™ library and Ion Torrent sequencing for the fast detection of genetic variations in human opioid receptor genes is described.

Methods

Gene selection

Exomic genotyping was performed for the human OPRM1, OPRK1 and OPRD1 genes (NCBI IDs 4985, 4988 and 4986), located on chromosomes 6, 8, and 1, and encoding for the three major (μ, κ, δ) opioid receptors. Their endogenous ligands are opioid neuropeptides comprising endorphins/endomorphines, dynorphins, enkephalins, respectively, and they are the main targets of opioid analgesics as listed in the DrugBank database (version 4.1; http://www.drugbank.ca [30]). In particular OPRM1 is addressed by most exogenous opioids. In addition, the SIGMAR1 gene (NCBI ID 10280) located on chromosome 9, was included that codes for σ1-opioid receptors. These are not activated or blocked by opioid peptides or naloxone, respectively, however, historically belong to the opioid receptors [31], are targets of some drugs classified as opioids such as pentazocine [32] or dextromethorphan [33] and their activation modulates peripheral μ-opioid analgesia [34]. Furthermore, the gene panel was extended to the opiate receptor-like 1 gene (OPRL1; NCBI ID 4987), located on chromosome 20 and coding for the nociceptin receptor that shows a high degree of structural homology to the classical opioid receptors although opioid peptides or morphine-like compounds have little or no affinity for it [35]. Its endogenous ligand is nociceptin, however, exogenous opioids comprising buprenorphine [36] and its glucuronide metabolites [37] have been shown to exert agonist activity that seems to contribute to their clinical effects, and the nonselective opioid agonist etorphine also binds at the OPRL1 gene product [38], which supports the inclusion of this gene in the present panel.

DNA template preparation and amplification

The investigation followed the Declaration of Helsinki on Biomedical Research Involving Human Subjects and was approved by the Ethics Committee of the Medical Faculty of the Goethe-University, Frankfurt, Germany. All subjects had provided informed written consent covering the genotyping. Genomic DNA was available from venous blood samples drawn from 79 pain patients in a previously reported [29] context of the pharmacogenomics of opioid dosing requirements.

DNA was extracted from 200 μl blood on a BioRobot EZ1 workstation applying the blood and body fluid spin protocol provided in the EZ1 DNA Blood 200 μl Kit (Qiagen, Hilden, Germany). A multiplex PCR amplification strategy for the coding DNA sequences was accomplished online (Ion Ampliseq™ Designer; http://www.ampliseq.com) to amplify the target region specified above (for primer sequences, see Supplementary Table 1) with 25 base pair exon padding. After a comparison of several primer design options, the design providing the maximum target sequence coverage was chosen. The ordered amplicons covered approximately 90% of the target sequence (Table 1). A total of 10 ng DNA per sample were used for the target enrichment by a multiplex PCR and each DNA pool was amplified with the Ion Ampliseq™ Library Kit in conjunction with the Ion Ampliseq™ “custom Primer Pool” - protocols according to the manufacturer procedures (Life Technologies, Darmstadt, Germany).

Table 1

AmpliSeq™ amplicons and coverage details of the human opioid receptor NGS assay.

Gene	Chr	Chr start	Chr end	Amplicons	Total bases	Covered bases	Coverage	Sum (total, covered, %)
OPRD1	Chr1	29138628	29139147	1	519	94	0.181	1924, 1336, 69.4%
	Chr1	29185440	29185840	2	400	292	0.730
	Chr1	29189228	29190233	6	1005	950	0.9453
OPRM1	Chr6	154331605	154331719	1	114	114	1.000	23787, 21859, 91.8%
	Chr6	154331605	154332178	3	573	573	1.000
	Chr6	154360211	154360465	2	254	254	1.000
	Chr6	154360270	154360465	2	195	195	1.000
	Chr6	154360520	154360994	4	474	474	1.000
	Chr6	154360211	154360994	6	783	783	1.000
	Chr6	154407616	154408945	8	1329	1329	1.000
	Chr6	154407616	154408964	8	1348	1348	1.000
	Chr6	154410935	154411338	3	403	403	1.000
	Chr6	154410871	154411338	3	467	467	1.000
	Chr6	154412061	154412632	4	571	571	1.000
	Chr6	154412918	154413098	2	180	180	1.000
	Chr6	154414379	154414680	2	301	301	1.000
	Chr6	154428725	154428812	1	87	87	1.000
	Chr6	154428574	154428812	2	238	238	1.000
	Chr6	154429125	154430178	6	1053	865	0.8215
	Chr6	154428857	154430178	9	1321	1123	0.8501
	Chr6	154431464	154431613	1	149	149	1.000
	Chr6	154439792	154453514	72	13722	12180	0.88763
	Chr6	154567801	154568026	2	225	225	1.000
OPRK1	Chr8	54138250	54142414	23	4164	3872	0.9299	5683, 5117, 90%
	Chr8	54147293	54147696	3	403	403	1.000
	Chr8	54155273	54155497	1	224	116	0.518
	Chr8	54163315	54163670	3	355	338	0.952
	Chr8	54163982	54164219	1	237	194	0.819
	Chr8	54163982	54164282	1	300	194	0.647
SIGMAR1	Chr9	34634693	34635662	6	969	969	1.000	3520, 3427, 97.3%
	Chr9	34634693	34635880	8	1187	1187	1.000
	Chr9	34636968	34637111	1	143	143	1.000
	Chr9	34636968	34637134	1	166	166	1.000
	Chr9	34637238	34637442	2	204	204	1.000
	Chr9	34637191	34637442	2	251	251	1.000
	Chr9	34637578	34637848	2	270	251	0.930
	Chr9	34637518	34637848	2	330	256	0.776
OPRL1	Chr20	62711445	62711730	2	285	230	0.807	3954, 3708, 93.7%
	Chr20	62716343	62716526	2	183	183	1.000
	Chr20	62717832	62717977	2	145	145	1.000
	Chr20	62723293	62723494	2	201	201	1.000
	Chr20	62724015	62724331	4	316	316	1.000
	Chr20	62729129	62729535	5	406	406	1.000
	Chr20	62729603	62732021	15	2418	2227	0.921

After each pool had undergone 17 PCR cycles, the PCR primers were removed with FuPa Reagent and the amplicons were ligated to the sequencing adaptors with short stretches of index sequences (barcodes) that enabled sample multiplexing for subsequent steps (Ion Xpress™ Barcode Adapters Kit; Life Technologies). After purification with AMPure XP beads (Beckman Coulter, Krefeld, Germany), the barcoded libraries were quantified with a Qubit® 2.0 Fluorimeter (Life Technologies, Darmstadt, Germany) and normalized for DNA concentration to a final concentration of 20 pmol/l using the Ion Library Equalizer™ Kit (Life Technologies, Darmstadt, Germany). Equalized barcoded libraries from 11 to 12 samples at a time were pooled. To clonally amplify the library DNA onto the Ion Sphere Particles (ISPs; Life Technologies, Darmstadt, Germany), the library pool was subjected to emulsion PCR by using an IT OneTouch template kit on an IT OneTouch system (Life Technologies, Darmstadt, Germany) following the manufacturer's protocol.

Sequencing

Enriched ISPs which carried many copies of the same DNA fragment were subjected to sequencing on an Ion 316 Chip to sequence pooled libraries with eleven to twelve samples. The 316 chip was chosen (instead of the low-capacity 314 or the high-capacity 318 chip) to obtain a mean sequencing depth of coverage of 50 × which means that, on average, each base has been sequenced 50 ×, when eleven samples were loaded. A larger number of samples could be analyzed simultaneously using the 318 chip, but this would increase the turnaround time for each sample, depending on the number of samples that are received by the laboratory. Sequencing was performed using the sequencing kit (Ion PGM 200 Sequencing Kit; Life Technologies, Darmstadt, Germany) as per the manufacturer's instructions with the 200-bp single-end run configuration.

Bioinformatics generation of sequence information

The raw data (unmapped BAM-files) from the sequencing runs were processed using Torrent Suite Software (Version 4.4.2, Life Technologies, Darmstadt, Germany) to generate read alignments which are filtered by the software into mapped BAM-files using the reference genomic sequence (hg19) of target genes. Variant calling was performed with the Torrent Variant Caller Plugin using as key parameters: minimum allele frequency = 0.015, minimum quality = 10, minimum coverage = 20 and minimum coverage on either strand = 3.

The annotation of called variants was done using the Ion Reporter Software (Version 4.4; Life Technologies, Darmstadt, Germany) for the VCF files that contained the nucleotide reads. The GenomeBrowse® software (Version 2.0.4, Golden Helix, Bozeman, MT, USA) was used to map the observed sequences to the reference sequences GRCh37 hg19 (dated February 2009). The SNP and Variation Suite software (Version 8.4.4; Golden Helix, Bozeman, MT, USA) was used for the analysis of sequence quality, coverage and also for variant identification.

Method validation

For method validation, two genomic regions were chosen at random from the whole analyzed region for validation by Sanger sequencing [25], [26] in an independent external laboratory (AGOWA, Berlin, Germany), which was performed in ten DNA samples randomly chosen from the n = 79 samples in the present cohort. Amplification of the respective DNA segments was done using PCR primer pairs (forward, reverse) of (i) 5´-ATGAAGACAGCAACCAACATTTAC-3´ and 5´-CCAGATGCAGATATTGATGATCTT-3´ and (ii) 5´-TTTAACCTGGCTTTGGCAGATG-3´ and 5´-ACATCGACGTCTTCCCTGACT-3´. The results of Sanger sequencing were aligned with the genomic sequence and analyzed using Chromas Lite® (Version 2.1.1, Technelysium Pty Ltd., South Brisbane, Australia) and the GenomeBrowse® (Version 2.0.4, Golden Helix, Bozeman, MT, USA) was used to compare the sequences obtained with NGS or Sanger techniques.

Results

The NGS assay of human opioid receptors was established on 79 genomic DNA samples [29]. As proposed previously [39], only exons and their boundary sequences for which read-depths > 20 for each nucleotide could be obtained were considered as successfully analyzed. Applying this criterion, complete or nearly complete coverage of the relevant sequences was obtained (Table 1; for details on missing variants, see Supplementary Table 2).

The sequencing of the whole cohort required seven separate runs with each 11–12 patients' samples. Coverage statistics (Table 1) were comparable among all runs and were in the range of accepted quality criteria [40], [41], [42]. During the seven runs, a median of 2.54 × 106 reads per run was generated. The mean depth was near from 200 reads, the mean read length evaluated 163 bases and average chip loading (Fig. 1) was 62%. To ensure a high density of ISPs on a chip and hence, a high sequencing output, the chip loading value should be ≥ 60% (Table 2). The observed NGS results matched with the results obtained with conventional sequencing of random samples (Fig. 2). In all validation samples, the correspondence between NGS and Sanger derived nucleotide sequences was 100%.

Fig. 1

Pseudo-color image of the Ion 316™ v2 Chip plate showing percent loading across the physical surface. This sequencing run had a 70% loading, which ensures a high ISP density. Every 316 chip contains more than 6 million wells and the color scale on the right side conduces as a loading indicator. Deep red coloration stays for a 100% loading which means that every well in this area contains an ISP (templated and non-templated) whereas deep blue coloration implies that the wells in this area are empty.

Table 2

Ion Torrent PGM Run Characteristics (n = 7) of the NGS of human opioid receptor genes in n = 79 DNA samples obtained from pain patients.

Chip	Samples on chip	Loading density (%)	Total bases	Total Reads	Mean read length (bp)
316_1	12	65	502 M	2,949,673	170
316_2	12	80	220 M	1,755,623	125
316_3	11	63	339 M	2,060,050	165
316_4	11	50	406 M	2,374,049	171
316_5	11	67	540 M	3,186,425	169
316_6	11	56	477 M	2,813,274	170
316_7	11	52	495 M	2,645,442	174

Fig. 2

The Figure shows the alignment of the ion torrent sequence of the κ-opioid receptor illustrated by Golden Helix Genome Browse® readout (A) versus the same sequence according to a Sanger electropherogram (B). The grey bars in panel A show the read depth of the single nucleotides (all nucleotide positions ≈ 300 reads). The green scale in panel A shows the amino acids from different transcripts in detail. The colored bars in Panel B indicate for a sequencing quality parameter, which means a high filled box means a trustworthy sequencing process.

From the NGS runs, a total of 35,447 nucleotides were read from each DNA sample. Following elimination of nucleotides agreeing with the standard human genome sequence GRCh37 hg19 (dated February 2009), the final result of the NGS consisted of a vector of nucleotide information for each individual DNA sample. This vector had a length equaling the number of chromosomal positions in which a non-reference nucleotide had been found in any probe of the actual study sample. In the 28 pain patients who required standard opioid doses for analgesia (for clinical details, see [29]) and could therefore substitute for a random sample with respect to opioid receptor genetics, these vectors had a length of 91 nucleotides (Fig. 3). They included the genetic information about known functional opioid receptor variants such as about the OPRM1 rs1799971 SNP (118 A > G; Chr6:154360797 in Fig. 3) as the so far scientifically most regarded variant possibly triggering reduced opioid efficacy or about the OPRM1 rs563649 SNP (Chr6:154407967 in Fig. 3) coding for the μ-opioid receptor splice variant MOR-1K probably triggering opioid-induced hyperalgesia [43]. However, the full NGS information comprises a whole matrix composed of the single vectors and the evaluation and interpretation of this information requires bioinformatics approaches able to unleash the advances of NGS over single variant approaches as shown elsewhere [29].

Fig. 3

Example opioid receptors' genetic pattern of 28 pain patients receiving standard opioid doses (extended from [29]). The plot shows the occurrence of variants (lines) per DNA sample (columns) as vectors of a length corresponding to the number of gene loci in which a non-reference nucleotide was found in any sample of the whole cohort. The vectors are composed of information about the number of non-reference alleles found at the respective locus in the respective sample, color codes as white, 0 non-reference alleles = wild type genotype, yellow, heterozygous, and red, 2 non-reference alleles. The opioid receptors are indicated by different color shading behind the vectors of individual genetic information.

Discussion

A comprehensive NGS assay for the exons and regulatory parts of human opioid receptors was developed that produced valid nucleotide sequences corresponding to those obtained with the classical Sanger sequencing technique. The NGS assay is suitable for small to large-scale experimental setups aiming at accessing the information about any nucleotide in a study cohort, with a selection of those that differ from the reference nucleotide.

Single opioid receptor variants have been causally involved in reduced responses to analgesic treatment [44], hyperalgesic responses to opioid analgesics [45], protection against opioid induced side effects [46], [47], the individual sensitivity to painful stimuli [48], [49], the risk for drug addiction [50], [51], the individual response to stress [12], the success of naloxone treatment of alcoholism [52], or even the susceptibility to breast cancer [53] or the incidence of sudden infant death syndrome [54]. At the molecular level, in particular μ-opioid receptor variants have been shown to reduce opioid receptor expression [16], [17], receptor signaling [17], [18], [19] up to an almost complete functional loss [18], [20], and alternative splicing with the appearance of a six-transmembrane segment μ-opioid receptor [45] that triggers excitatory effects [43]. However, despite this apparently successful research on opioid receptor genetics, results have neither entered clinical practice nor have they been included in main treatment guidelines for pain or addiction [24], [55]. This may have several causes such as too small effect sizes counterbalanced by concomitantly other genetic variants [56]. However, an important perception of the discrepancy between the pathophysiological importance of opioid receptors including their role as main drug targets and the modest utility of their genotyping is, that the limited selection of published functional variants with accessibility by specific genetic assays [23] is probably insufficient. Indeed, we have recently shown that NGS of human opioid receptor genes outperforms single opioid receptor variants as a genetic biomarker of opioid-analgesia related phenotypes [29].

Research interest in the complete genomic information dates back to the seventies of the last century when the selective incorporation of chain-terminating dideoxynucleotides by DNA polymerase during in vitro DNA replication had been introduced [25], [26]. Techniques significantly improved during the last decades [57] with the development of contemporary machines in the late 1990s released to the market around the year 2005. The term “next generation” DNA sequencing refers to high-throughput technologies capable of parallel analyses of large numbers of different DNA sequences in a single reaction [58]. NGS has been attributed the potential to accelerate biomedical research [28], [59], [60]. Current leading manufacturers of commercial NGS platforms are Illumina (Genome Analyzer, HiSeq, MiSeq) and Life Technologies (SOLiD System, Ion Torrent, Proton). The systems combine conceptually similar workflows, starting with the creation of the genetic sample, which commences library preparation involving fragmentation of genomic DNA, purifying to uniform and desired fragment size and ligation to sequencing adapters specific to the platform. Differences apply to the reaction biochemistry and the way how the sequencing information is read [58]. In the present ion semiconductor sequencing method, libraries are immobilized to beads and amplified in microdroplets of aqueous solution and oil using emulsion PCR. Individual nucleotide bases are incorporated via DNA polymerase, which in the case of success triggers the release of a proton. The semiconductor chip that acts as a pH meter [61] providing the final readout. Alternative techniques use the detection of light instead, i.e., from optical fluorescence signals in the case of successful nucleotide incorporation the DNA nucleotide sequence is assembled. The different techniques differ with respect to the obtained throughput and accuracy [62], [63].

The high throughput and comprehensive information about DNA sequences are presently reflected in the assay costs. The sequencing of the opioid receptor genes of 79 patients required € 4500 for the AmpliSeq™ custom panel, € 5800 for library preparation, € 980 for template preparation and € 1200 for sequencing. In addition, approximately € 800 was spent for consumables and supplies. With twelve barcoded samples loaded on each chip, analysis costs for a single patient's opioid receptor gene sequence were approximately € 168. NGS costs are expected to quickly fall in near future [64]. However, despite this rapid technological progress, the analysis of the generated large data sets remains challenging [65]. As the sequencing process is only the beginning of the procedure, the analysis of the resulting “big data” requires substantial computational power, bioinformatics expertise and “up to date” databases of genomic variations. NGS technologies seem to shift the workload essentially away from the laboratory sample preparation toward various data analysis processes.

Research on the genetic variation of human opioid receptors as the main targets of endogenous and exogenous opioids is of immediate interest for assessing the clinical effects of opioid analgesics and for studying the epidemiology of substance addiction. The present NGS method is suitable for large-scale sequencing of an extended set of human genes belonging or related to opioid receptors. By covering almost the complete relevant coding and regulatory parts of the genes, the method includes all variants studied so far for functional associations and adds a large amount of genetic information as a basis for complete analysis of human opioid receptor genetics and its functional consequences.

Conflict of interest statement

The authors have declared that no further conflicts of interest exist.

Author contributions

Conceived and designed the experiments: JL, DK. Performed the experiments: DK. Analyzed the data: JL, DK. Wrote the paper: JL, DK.

Funding

This work has been funded by the European Union Seventh Framework Programme (FP7/2007–2013) under grant agreement no. 602919 (“GLORIA”, JL). An additional contribution, in particular the funding of the next generation genotyping equipment, has been received from the Landesoffensive zur Entwicklung wissenschaftlich-ökonomischer Exzellenz (LOEWE), Zentrum: Translational Medicine and Pharmacology (JL). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.