Deep sequencing analysis to identify novel and rare variants in pain-related genes in patients with acute postoperative pain and high morphine use.

PURPOSE: Most of the genetic variants that are reported to be associated with common pain phenotypes and analgesic use are common polymorphisms. The objective of our study was to identify new variants and investigate less common genetic variants that are usually not included in either small single-gene studies or high-throughput genotyping arrays. PATIENTS AND METHODS: From a cohort of 1075 patients who underwent a scheduled total abdominal hysterectomy, 92 who had higher self-rated pain scores and used more morphine were selected for the re-sequencing of 105 genes.
RESULTS: We identified over 2400 variants in 104 genes. Most were intronic with frequencies >5%. There were 181 novel variants, of which 30 were located in exons: 17 nonsynonymous, 10 synonymous, 2 non-coding RNA, and 1 stop-gain. For known variants that are rare (population frequency <1%), the frequencies of 54 exonic variants and eight intronic variants for the sequenced samples were higher than the weighted frequencies in the Genome Aggregation Database for East and South Asians (P-values ranging from 0.000 to 0.046). Overall, patients who had novel and/or rare variants used more morphine than those who only had common variants.
CONCLUSION: Our study uncovered novel variants in patients who reported higher pain and used more morphine. Compared with the general population, rare variants were more common in this group.

Introduction

Acute pain is the body’s mechanism to signal tissue injury and danger. Although pain helps to protect against further tissue damage by altering host behavior, prolonged and persistent pain has little biological value. On the contrary, it has an adverse impact on a person’s psychosocial well-being. The anticipation of pain can also influence the patient’s willingness to undergo potentially beneficial medical treatments that may be perceived as painful. Any pain that persists after surgery or injury carries adverse health and socioeconomic impacts, reduces the quality of life, increases health care cost and decreases work attendance.1

Pain perception is highly subjective with wide inter-individual variability in its sensitivity and tolerance. Known biological factors that impact this perception include age, race, gender, physiology, and social and psychological status.2–5 Pain is also a heritable phenotype, with multiple lines of evidence from Mendelian pain disorders, twin studies and increased risk for chronic pain conditions in individuals with family history. For instance, mutations in SCN9A and related genes have been identified in autosomal recessive congenital indifference to pain (MIM #243000) and autosomal dominant Marsili syndrome (MIM#147430).6–9 For less extreme and more complex pain phenotypes, the genetic contribution to sensitivity variation for different types of pain varies from 22% to 60%.10–12 For chronic pain conditions, twin studies suggest heritability of 39–58% for neuropathic pain,13–15 46% for chronic pelvic pain16 and as high as 70% for low back pain.17

Candidate gene studies have uncovered the contribution of variants of genes in the pain pathways across different types of pain in multiple populations. With the advent of genome-wide association studies (GWAS), the number of variants and chromosomal loci associated with pain has been further expanded.18–20 Published results from various pain studies are captured in several online databases such as the "Pain Genes Database of pain-related transgenic knockout studies" (PainGenes db)21 and the "Human Pain Genetics Database".22 Due to the study design and limitations in statistical power, most of the identified variants have been common genetic polymorphisms. These common variants tend to have only small to moderate impact on the difference in quantitative measures of pain. Furthermore, rare and low-frequency variants have been suggested to account for the remaining heritability.23–26

To uncover novel and rare variants that might be enriched in individuals who experienced more intense pain, we re-sequenced 105 genes in 92 patients who self-reported higher postoperative pain or used more morphine. They were selected from a cohort of patients who underwent total abdominal hysterectomy in our hospital. Our results showed that these patients had higher frequencies of rare variants in pain-related genes compared with those from population databases.

Patients and methods

Our study was approved by the SingHealth Central Institutional Review Board and conducted in accordance with the Declaration of Helsinki. Written informed consent for genetic study was obtained from all patients prior to surgical procedure.

Subjects characteristics, pain assessment and sample collection

The study protocol for this prospectively recruited cohort of 1075 women who underwent planned total hysterectomy at the KK Women’s and Children’s Hospital has been described previously.27 Briefly, pain sensitivity and tolerance were determined preoperatively using the blood pressure cuff of a sphygmomanometer. The cuff was placed around the patient’s upper arm and inflated until she indicated pain. The mercury reading (in mm) at that point was taken as the pain threshold. Pain tolerance was recorded as the mercury reading at which the patient requested for the deflation of the cuff. Immediately after surgery completion, the patient was fitted with a patient-controlled analgesia pump (PCA) that was set to deliver an intravenous bolus of 1 mg morphine on demand, with lockout interval of 5 mins, no basal infusion and a maximum hourly dose of 10 mg morphine. At 4-hourly intervals, patients were asked to rate their pain according to the VAS (0=no pain, 10=worst pain imaginable), as well as pruritus and nausea on a scale of 0–3 (0=none, 1=mild, 2=moderate, 3=severe).

For sequencing analysis, we selected from 1047 patients who had complete morphine data and pain scores for the 24-hr postoperative period. Tukey fence analysis was applied to select patients with outlier acute pain profiles. Since our interest was on higher pain scores, we only selected the upper fence. Fourteen outliers were selected based on acute pain scores at 4 hrs and average pain scores. To increase the sample size, we selected additional 50 patients whose acute pain scores at 4 hrs were greater than the third quartile. Despite not reporting pain scores greater than the third quartile, additional 4 patients were included based on higher outlier morphine consumption. Lastly, we also included 41 patients with 4-hr pain scores in the third quartile and had 8-hr pain scores that were less than the first quartile. We assumed that these patients had higher acute pain but also rapid resolution. Our final list had 109 patients arranged according to the date of surgery. Of these, the first 92 on the list with adequate good quality DNA were used for preparing sequencing libraries. The demographic and clinical characteristics of the 92 patients who were sequenced and those who were not are shown in Table 1. From the medical record, 35 of the 91 patients had one or more chronic conditions, of which the most common was hypertension (19 patients) followed by diabetes (7 patients). Only one patient had a pain condition (migraine). None of the patients were on opioid medication.

Table 1

Characteristics of samples selected and not selected for sequencing

Variable	Sequenced	Not sequenced	P-value
Age	(n=92)	(n=955)
Mean (SD)	47.4 (6.0)	47.8 (5.3)	0.427
Median	47.0	48.0
Min, Max	34, 76	30, 78
Ethnicity n (%)	(n=92)	(n=955)
Chinese	69 (75.0)	686 (71.8)	0.263
Malay	18 (19.6)	166 (17.4)
Indian	5 (5.4)	103 (10.8)
BMI	(n=92)	(n=955)
Mean (SD)	23.88 (4.01) 23.89 (4.02)	24.76 (4.17)	0.053
Median	23.62 23.62	24.32
Min, Max	15.94, 36.67	15.56, 38.22
Pain threshold (mmHg)	(n=85)	(n=875)
Mean (SD)	240.66 (43.31)	245.19 (43.95)	0.364
Median (IQR)	250.00	250.00
Min, Max	80, 300	100, 300
Pain tolerance (mmHg)	(n=85)	(n=875)
Mean (SD)	275.51 (27.55)	282.43 (23.49)	0.011
Median	290.00	290.00
Min, Max	180, 300	170, 300
Time-averaged VAS	(n=92)	(n=955)
Mean (SD)	1.55 (0.89)	1.20 (0.86)	0.000
Median	1.50	1.00
Min, Max	0.00, 3.83	0.00, 9.33
PCA morphine	(n=92)	(n=955)
Mean (SD)	20.93 (12.49)	16.31 (12.19)	0.001
Median	21.50	14.00
Min, Max	1, 50	0, 71
PCA morphine/weight (mg/kg)	(n=92)	(n=955)
Mean (SD)	359.60 (219.06)	271.80 (200.55)	0.000
Median	343.85	229.51
Min, Max	17.24, 917.43	0.00, 1116.67

DNA sequencing

Genomic DNA was extracted in batches from frozen whole blood samples in EDTA tubes using the Gentra Puregene Blood Kit (Qiagen, Hilden, Germany). DNA was checked for quantity and purity using the Quawell Q5000 UV-Vis Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA).

The 105 target genes (Table 2) were selected based on published literature and the maximum target size for the chosen sequencing platform. SureSelect and HaloPlex advanced wizards (Agilent Technologies, Santa Clara, CA, USA) were used to design the capture probes for target regions. Genomic coordinates for specified targets were obtained from RefSeq, Ensembl, CCDS, Gencode, VEGA, SNP, and CytoBand genome annotation databases, using the H. sapiens hg19 (GRCh37) as the reference sequence. All coding exons with minimum extensions of 10 bases from both 3ʹ and 5ʹ ends of each exon were included. The design covered 99.47% of the target region using 12,776 amplicons. The total size of the amplicons was 637,374 kilobases (kb), with total analyzable target of 234,538 kb.

Table 2

List of pain-related genes sequenced in this study

Gene	Full name	Chr	MIM#
ABCB1	ATP-BINDING CASSETTE, SUBFAMILY B, MEMBER 1	7	171050
ADORA1	ADENOSINE A1 RECEPTOR	1	102775
ADRB2	BETA-2-ADRENERGIC RECEPTOR	5	109690
ANKK1	ANKYRIN REPEAT- AND KINASE DOMAIN-CONTAINING PROTEIN 1	11	608774
ATP1A2	ATPase, Na+/K+ TRANSPORTING, ALPHA-2 POLYPEPTIDE	1	182340
ATP1A3	ATPase, Na+/K+ TRANSPORTING, ALPHA-3 POLYPEPTIDE	19	182350
CACNA1B	CALCIUM CHANNEL, VOLTAGE-DEPENDENT, N TYPE, ALPHA-1B SUBUNIT	9	601012
CACNG2	CALCIUM CHANNEL, VOLTAGE-DEPENDENT, GAMMA-2 SUBUNIT	22	602911
CCNJL	CYCLIN J LIKE	5	NA
CD4	CD4 ANTIGEN	12	186940
CHRNA4	CHOLINERGIC RECEPTOR, NEURONAL NICOTINIC, ALPHA POLYPEPTIDE 4	20	118504
CNR1	CANNABINOID RECEPTOR 1	6	114610
CNR2	CANNABINOID RECEPTOR 2	1	605051
COMT	CATECHOL-O-METHYLTRANSFERASE	22	116790
CREB1	cAMP RESPONSE ELEMENT-BINDING PROTEIN 1	2	123810
CYP19A1	CYTOCHROME P450, FAMILY 19, SUBFAMILY A, POLYPEPTIDE 1	15	107910
CYP2C19	CYTOCHROME P450, SUBFAMILY IIC, POLYPEPTIDE 19	10	124020
CYP2C9	CYTOCHROME P450, SUBFAMILY IIC, POLYPEPTIDE 9	10	601130
CYP2D6	CYTOCHROME P450, SUBFAMILY IID, POLYPEPTIDE 6	22	124030
CYP3A4	CYTOCHROME P450, SUBFAMILY IIIA, POLYPEPTIDE 4	7	124010
CYP3A5	CYTOCHROME P450, SUBFAMILY IIIA, POLYPEPTIDE 5	7	605325
DDAH1	DIMETHYLARGININE DIMETHYLAMINOHYDROLASE 1	1	604743
DLD	DIHYDROLIPOAMIDE DEHYDROGENASE	7	238331
DNM2	DYNAMIN 2	19	602378
DRD2	DOPAMINE RECEPTOR D2	11	126450
EPHX1	EPOXIDE HYDROLASE 1, MICROSOMAL	1	132810
ESR1	ESTROGEN RECEPTOR 1	6	133430
ESR2	ESTROGEN RECEPTOR 2	14	601663
FBXW7	F-BOX AND WD40 DOMAIN PROTEIN 7	4	606278
FKBP4	FK506-BINDING PROTEIN 4	12	600611
FLOT1	FLOTILLIN 1	6	606998
GCH1	GTP CYCLOHYDROLASE I	14	600225
GDAP1	GANGLIOSIDE-INDUCED DIFFERENTIATION-ASSOCIATED PROTEIN 1	8	606598
GRIK4	GLUTAMATE RECEPTOR, IONOTROPIC, KAINATE 4	11	600282
GRIN1	GLUTAMATE RECEPTOR, IONOTROPIC, N-METHYL-D-ASPARTATE, SUBUNIT 1	9	138249
GRIN2B	GLUTAMATE RECEPTOR, IONOTROPIC, N-METHYL-D-ASPARTATE, SUBUNIT 2B	12	138252
GRM1	GLUTAMATE RECEPTOR, METABOTROPIC, 1	6	604473
GRM5	GLUTAMATE RECEPTOR, METABOTROPIC, 5	11	604102
HINT1	HISTIDINE TRIAD NUCLEOTIDE-BINDING PROTEIN 1	5	601314
HLA-B	MAJOR HISTOCOMPATIBILITY COMPLEX, CLASS I, B	6	142830
HTR1A	5-HYDROXYTRYPTAMINE RECEPTOR 1A	5	109760
HTR2A	5-HYDROXYTRYPTAMINE RECEPTOR 2A	13	182135
HTR2C	5-HYDROXYTRYPTAMINE RECEPTOR 2C	X	312861
IFI30	INTERFERON-GAMMA-INDUCIBLE PROTEIN 30	19	604664
IL10	INTERLEUKIN 10	1	124092
IL18	INTERLEUKIN 18	11	600953
IL1A	INTERLEUKIN 1-ALPHA	2	147760
IL1B	INTERLEUKIN 1-BETA	2	147720
IL2	INTERLEUKIN 2	4	147680
IL6	INTERLEUKIN 6	7	147620
KCNIP3	POTASSIUM CHANNEL-INTERACTING PROTEIN 3	2	604662
KCNJ6	POTASSIUM CHANNEL, INWARDLY RECTIFYING, SUBFAMILY J, MEMBER 6	21	600877
KCNQ2	POTASSIUM CHANNEL, VOLTAGE-GATED, KQT-LIKE SUBFAMILY, MEMBER 2	20	602235
KCNQ3	POTASSIUM CHANNEL, VOLTAGE-GATED, KQT-LIKE SUBFAMILY, MEMBER 3	8	602232
KCNS1	POTASSIUM CHANNEL, VOLTAGE-GATED, DELAYED-RECTIFIER, SUBFAMILY S, MEMBER 1	20	602905
KIF5A	KINESIN FAMILY MEMBER 5A	12	602821
LTA	LYMPHOTOXIN-ALPHA	6	153440
MAOA	MONOAMINE OXIDASE A	X	309850
MAOB	MONOAMINE OXIDASE B	X	309860
MAPK1	MITOGEN-ACTIVATED PROTEIN KINASE 1	22	176948
MC1R	MELANOCORTIN 1 RECEPTOR	16	155555
MTCO2	COMPLEX IV, CYTOCHROME c OXIDASE SUBUNIT II	M	516040
MYPN	MYOPALLADIN	10	608517
NGF	NERVE GROWTH FACTOR	1	162030
NOTCH3	NOTCH, DROSOPHILA, HOMOLOG OF, 3	19	600276
NTRK1	NEUROTROPHIC TYROSINE KINASE, RECEPTOR, TYPE 1	1	191315
OPRD1	OPIOID RECEPTOR, DELTA-1	1	165195
OPRK1	OPIOID RECEPTOR, KAPPA-1	8	165196
OPRM1	OPIOID RECEPTOR, MU-1	6	600018
OR5F1	OLFACTORY RECEPTOR, FAMILY 5, SUBFAMILY F, MEMBER 1	11	608492
OXT	OXYTOCIN	20	167050
OXTR	OXYTOCIN RECEPTOR	3	167055
P2RX3	PURINERGIC RECEPTOR P2X, LIGAND-GATED ION CHANNEL, 3	11	600843
P2RX4	PURINERGIC RECEPTOR P2X, LIGAND-GATED ION CHANNEL, 4	12	600846
P2RX7	PURINERGIC RECEPTOR P2X, LIGAND-GATED ION CHANNEL, 7	12	602566
P2RY1	PURINERGIC RECEPTOR P2Y, G PROTEIN-COUPLED, 1	3	601167
PDHA2	PYRUVATE DEHYDROGENASE, ALPHA-2	4	179061
PMP22	PERIPHERAL MYELIN PROTEIN 22	17	601097
POLG	POLYMERASE, DNA, GAMMA	15	174763
PRRT2	PROLINE-RICH TRANSMEMBRANE PROTEIN 2	16	614386
PTGS2	PROSTAGLANDIN-ENDOPEROXIDE SYNTHASE 2	1	600262
RAMP1	RECEPTOR ACTIVITY-MODIFYING PROTEIN 1	2	605153
RHEB	RAS HOMOLOG ENRICHED IN BRAIN	7	601293
SCN10A	SODIUM CHANNEL, VOLTAGE-GATED, TYPE X, ALPHA SUBUNIT	3	604427
SCN11A	SODIUM CHANNEL, VOLTAGE-GATED, TYPE XI, ALPHA SUBUNIT	3	604385
SCN1A	SODIUM CHANNEL, NEURONAL TYPE I, ALPHA SUBUNIT	2	182389
SCN3A	SODIUM CHANNEL, VOLTAGE-GATED, TYPE III, ALPHA SUBUNIT	2	182391
SCN9A	SODIUM CHANNEL, VOLTAGE-GATED, TYPE IX, ALPHA SUBUNIT	2	603415
SLC1A3	SOLUTE CARRIER FAMILY 1 (GLIAL HIGH AFFINITY GLUTAMATE TRANSPORTER), MEMBER 3	5	600111
SLC2A1	SOLUTE CARRIER FAMILY 2 (FACILITATED GLUCOSE TRANSPORTER), MEMBER 1	1	138140
SLC6A2	SOLUTE CARRIER FAMILY 6 (NEUROTRANSMITTER TRANSPORTER, NORADRENALINE), MEMBER 2	16	163970
SLC6A3	SOLUTE CARRIER FAMILY 6 (NEUROTRANSMITTER TRANSPORTER, DOPAMINE), MEMBER 3	5	126455
SLC6A4	SOLUTE CARRIER FAMILY 6 (NEUROTRANSMITTER TRANSPORTER, SEROTONIN), MEMBER 4	17	182138
TAGAP	T-CELL ACTIVATION GTPase-ACTIVATING PROTEIN	6	609667
TBK1	TANK-BINDING KINASE 1	12	604834
TH	TYROSINE HYDROXYLASE	11	191290
TNF	TUMOR NECROSIS FACTOR	6	191160
TNFRSF1A	TUMOR NECROSIS FACTOR RECEPTOR SUPERFAMILY, MEMBER 1A	12	191190
TPH2	TRYPTOPHAN HYDROXYLASE 2	12	607478
TRPA1	TRANSIENT RECEPTOR POTENTIAL CATION CHANNEL, SUBFAMILY A, MEMBER 1	8	604775
TRPV1	TRANSIENT RECEPTOR POTENTIAL CATION CHANNEL, SUBFAMILY V, MEMBER 1	17	602076
TRPV3	TRANSIENT RECEPTOR POTENTIAL CATION CHANNEL, SUBFAMILY V, MEMBER 3	17	607066
TTR	TRANSTHYRETIN	18	176300
UGT2B15	URIDINE DIPHOSPHATE GLYCOSYLTRANSFERASE 2 FAMILY, MEMBER B15	4	600069
ZNF767P	ZINC FINGER FAMILY MEMBER 767, PSEUDOGENE	7	NA

Abbreviations: M, mitochondria; N.A., not available.

The HaloPlex Target Enrichment System (version F1) was used to index the samples and amplify the target regions according to the manufacturer’s instruction (Agilent Technologies). Libraries produced from the 92 samples were sequenced using 250 bp paired-end sequencing (600-cycle) on one MiSeq Reagent Kit (v3) on a MiSeq System (Illumina, San Diego, CA, USA).

Data processing and analysis

Bases were called using the on-instrument MiSeq Reporter software (version 2.6). Alignment processing and variant calling were performed with reference to human genome GRCh37 (hg19). The variant call format file generated was annotated and prioritized using wANNOVAR.28 Variants were considered novel if they were not previously reported in Genome Aggregation Database (gnomAD), Exome Sequencing Project, Human Genetic Variation Database, ClinVar, 1000 Genomes, or Human Gene Mutation Database databases, and not documented in scientific literature.

Consequences of sequence changes were assessed using Alamut Visual software version 2.10 (Interactive Biosoftware, Rouen, France) that included in silico prediction algorithms for likely effect on amino acid substitutions (SIFT v6.2.0, and PolyPhen-2 v2.2.2r398). Nonsynonymous variants with SIFT scores of <0.05 were classified as “deleterious”.29 For PolyPhen‐2, scores of >0.85 were classified as “probably damaging”, and scores of 0.15–0.85 were considered as “possibly damaging”.30

Two programs (MaxEnt and NNSPLICE) were used to evaluate the potential effect on splicing. Variants were considered positive if one or both programs had variation in the splice site score greater than the cutoff value of 10% from that of the reference allele.31

For rare variants (population frequencies of <1%), only exonic variants that are not synonymous, and intronic variants with predicted splice effects were compared with corresponding frequencies in gnomAD r2.0.2. Since our sequenced samples comprised 69 Chinese (75.0%), 18 Malays (19.6%) and 5 Indians (5.4%), analysis on statistically significant difference was performed with weighted gnomAD frequencies calculated from both the East Asian and the South Asian populations with respective weightings of 94.6% and 5.4%.

Interaction network and enrichment analyses

For genes with identified rare and/or novel variants of functional consequence, their involvement in biological pathways was queried using STRING database (version 10.5) (https://string-db.org/)32 that contains known and predicted protein interactions. We used Kyoto Encyclopedia of Genes and Genomes (KEGG) to assess network representation and for biological interpretation of the network nodes. Pathways with P-values <0.001 after false discovery rate adjustment were considered statistically enriched.

Statistical analyses

One-way ANOVA test was used to compare quantitative variables between groups, with Tukey post hoc test for comparison of more than two groups. Chi-square or Fisher’s exact test was used to compare frequencies for categorical variables. Post hoc Bonferroni test for multiple comparisons was performed for the comparison of the 62 rare variants for P-value correction. All statistical analyses were performed using IBM SPSS Statistics 19, with P-values ≤0.05 considered as statistically significant. For association analysis with rare variants, the P-value cutoff would be 0.00083 after applying Bonferroni correction for multiple testing.

Results

Quality of next-generation sequencing

Of the 92 samples, one failed to produce sequence data output. For the remaining 91 samples, 97.63% of the reads aligned to the reference genome (GRCh37/hg19) and 95.11% of the reads mapped to the targeted regions, with mean region coverage depth of 157.1× (Table S1). The mean coverage of targeted bases was 88.15% and 67.22% at 20× and 50×, respectively (Table S2).

Table S1

Summary of sequencing quality and output for the 91 samples sequenced

Number of reads			Enrichment	Mean coverage
	Total	Aligned
Mean	741,168	97.63%	95.11%	157.1×
Median	714,188	97.60%	95.20%	148.5×
Lowest	502,860	93.00%	93.50%	109.5×
Highest	1,079,530	99.10%	96.60%	243.1×
Number of bases sequenced
	Total	Aligned	Enrichment	Q30
Mean	106,485,318	96.06%	96.42%	92.16%
Median	101,200,682	96.40%	96.50%	92.10%
Lowest	73,331,116	92.60%	95.60%	91.00%
Highest	163,227,534	97.40%	97.50%	92.80%

Table S2

Percentage of bases sequenced at the different read depths

	Target base coverage at read depths
1×	10×	20×	50×
Mean	97.84%	90.44%	83.15%	67.22%
Median	97.90%	90.60%	83.20%	66.80%
Lowest	93.80%	84.20%	76.90%	58.90%
Highest	98.50%	93.10%	88.50%	76.10%

At the gene level, all 105 targeted genes had mean coverage of at least 30× even for the gene with the lowest coverage. Eighty-one genes had a mean of >100×. Except for TBK1 which had the lowest mean coverage of ~49×, the remaining 104 genes had mean coverage of at least 62×. The mitochondrial gene MTCO2 had the highest mean coverage (>13,000×), followed by CNR2 (349×) (Table S3). Despite the high mean target gene coverage, amplification failed in at least one sample for 6 of 1014 target regions. Five genes (ADRB2, CHRNA4, HLA-B, TNFRSF1A, and TRPV3) had at least one region that was not amplified and therefore not sequenced. There were also 18 target regions from 13 genes (ADRB2, ATP1A3, CYP2C19, GRM5, RAMP1, SCN1A, SCN3A, SCN9A, SLC1A3, TBK1, TNFRSF1A, TRPA1, UGT2B15) with read depth of <20×.

Table S3

Coverage details for each gene

Gene	Mean	Median	Lowest	Highest
ABCB1	107.3	103.5	71.3	156.7
ADORA1	294.6	281.1	212.1	469.5
ADRB2	219.3	210.3	157.4	370.4
ANKK1	157.3	148.8	115.4	250.0
ATP1A2	162.5	153.5	110.0	253.2
ATP1A3	154.7	146.8	109.5	238.5
CACNA1B	164.1	153.6	116.8	260.1
CACNG2	209.6	199.0	138.3	331.9
CCNJL	170.0	160.6	120.8	268.4
CD4	163.2	157.0	116.6	257.7
CHRNA4	164.3	158.5	116.3	243.6
CNR1	246.0	235.6	174.2	377.6
CNR2	348.7	330.7	258.4	552.3
COMT	185.6	175.5	133.3	289.5
CREB1	90.9	90.8	60.4	125.9
CYP19A1	178.5	173.4	123.5	267.6
CYP2C19	154.9	147.1	107.7	240.9
CYP2C9	121.2	117.7	81.7	187.9
CYP2D6	182.9	173.0	62.3	442.4
CYP3A4	107.4	103.5	75.5	163.2
CYP3A5	95.0	92.0	63.6	148.7
DDAH1	100.8	98.5	65.9	155.9
DLD	79.8	78.1	48.9	113.9
DNM2	128.2	121.3	93.0	199.5
DRD2	146.0	137.7	102.7	245.0
EPHX1	168.5	161.1	116.4	269.8
ESR1	184.6	175.1	129.4	276.0
ESR2	122.8	116.7	84.5	189.2
FBXW7	88.3	86.3	60.3	126.7
FKBP4	128.5	123.3	88.2	199.1
FLOT1	158.6	152.2	108.4	252.8
GCH1	100.4	95.3	72.1	154.3
GDAP1	147.7	142.1	103.8	223.1
GRIK4	128.9	121.4	91.4	202.1
GRIN1	113.6	109.0	81.1	184.0
GRIN2B	165.5	156.0	115.1	260.9
GRM1	165.9	156.9	119.1	250.9
GRM5	123.0	118.8	82.5	186.0
HINT1	121.6	115.7	76.3	193.3
HLA-B	112.6	111.7	67.6	177.6
HTR1A	226.4	211.5	154.9	381.7
HTR2A	116.1	110.7	76.0	182.1
HTR2C	126.2	120.7	82.2	196.9
IFI30	180.9	169.2	128.7	297.4
IL10	195.4	187.0	139.4	305.7
IL18	62.2	61.1	36.0	99.7
IL1A	97.3	95.0	59.2	153.9
IL1B	155.2	149.7	111.7	243.1
IL2	64.8	63.9	42.9	94.4
IL6	113.4	107.8	81.5	183.7
KCNIP3	127.6	121.7	88.6	197.3
KCNJ6	166.7	158.0	106.0	265.5
KCNQ2	130.4	124.3	94.8	204.3
KCNQ3	132.1	123.7	93.1	209.4
KCNS1	149.1	144.6	109.3	242.5
KIF5A	168.1	160.8	120.5	267.6
LTA	253.6	240.1	178.3	402.3
MAOA	85.7	82.2	57.6	134.5
MAOB	92.5	88.9	65.7	142.2
MAPK1	102.8	99.4	67.1	151.6
MC1R	233.3	217.9	167.9	404.7
MTCO2	13,303.0	13,457.6	5704.4	24,023.5
MYPN	109.2	106.3	75.2	167.2
NGF	298.3	283.2	199.5	464.6
NOTCH3	153.6	144.6	110.5	250.7
NTRK1	174.4	164.9	126.9	282.6
OPRD1	135.7	132.0	96.9	199.2
OPRK1	127.4	118.2	92.0	201.9
OPRM1	137.8	131.4	95.4	208.2
OR5F1	194.3	190.2	127.0	292.7
OXT	134.4	127.9	89.6	207.5
OXTR	167.5	160.9	127.7	260.7
P2RX3	152.1	146.8	102.2	246.7
P2RX4	174.7	167.4	124.5	268.6
P2RX7	159.5	151.0	114.3	245.3
P2RY1	155.7	147.3	101.4	246.6
PDHA2	259.0	252.4	189.9	412.3
PMP22	132.9	129.3	98.4	190.6
POLG	173.7	165.9	124.1	270.3
PRRT2	256.9	246.0	174.3	406.6
PTGS2	101.9	100.1	69.2	143.0
RAMP1	158.9	153.2	92.4	255.8
RHEB	85.4	83.1	55.3	134.3
SCN10A	171.3	161.8	122.7	267.2
SCN11A	109.1	105.0	74.0	163.6
SCN1A	76.4	75.4	51.5	110.1
SCN3A	90.8	88.8	59.5	137.4
SCN9A	92.2	90.3	62.4	132.9
SLC1A3	130.1	122.5	89.3	206.1
SLC2A1	155.9	146.3	109.1	249.5
SLC6A2	182.4	172.8	125.8	291.6
SLC6A3	182.4	171.6	132.9	289.3
SLC6A4	171.2	166.9	123.0	255.9
TAGAP	154.0	145.8	102.7	235.8
TBK1	48.5	48.1	32.8	67.1
TH	161.4	154.1	113.6	245.2
TNF	246.5	235.9	159.6	382.6
TNFRSF1A	125.4	116.3	90.2	194.8
TPH2	140.9	136.2	93.7	221.4
TRPA1	73.9	72.6	49.4	112.3
TRPV1	154.2	146.7	110.4	239.5
TRPV3	149.1	141.2	103.4	238.6
TTR	180.9	173.0	125.2	277.9
UGT2B15	67.1	67.4	39.8	108.9
ZNF767P	141.5	134.0	99.2	218.9

Summary of genetic variants identified

In total, 2466 variants were identified from 104 genes. Only the mitochondrial gene MTCO2 had no variant. Most of the variants were common (population frequency >5%) and low-frequency polymorphisms (frequency 1–5%), the remaining comprised 608 rare (frequency <1%) and 181 novel variants (defined as those with no Reference SNP numbers and not documented in databases or published literature). In terms of location, the largest number of 1477 were found in introns, followed by 771 in protein-coding exons, 123 in 3ʹ untranslated regions (or trailer sequences), 45 in 5ʹ untranslated regions (or leader sequences), 35 in upstream regions of genes, 12 in the downstream regions, and the remaining 3 in intergenic regions. Overall, there was an average of 27.1 variants per patient.

For single-nucleotide substitutions located in the exons, 386 were synonymous while 350 were missense variants. There were also 5 stop-gain variants and 1 stop-loss variant. For changes involving multiple nucleotides, there were 2 non-frameshift insertions, 1 frameshift insertion, and 4 non-frameshift deletions. In addition, there were 21 exonic non-coding RNA variants. The position of a putative OPRK1 variant (chr8:54141824:C>T) within the gene could not be determined.

Analysis of novel and rare variants

There were 181 novel variants in 70 genes, most of which were in the introns. Of the 30 variants found in exons, 17 were missense, 10 synonymous, 1 stop-gain, and 2 were non-coding RNA. The list of 30 exonic variants and two intronic variants with their predicted consequences are listed in Table 3, along with the number of reads for novel/alternate alleles and their corresponding reference alleles. The 32 novel variants were from 28 patients. The numbers of reads for the 2 alleles were mostly balanced. Hence, we did not perform Sanger validation.

Table 3

List of novel exonic and intronic variants (with splicing effect) identified in this population

Gene	GenBank ref	Variant	Alamut visual prediction		# reads
Exonic	PolyPhen-2	SIFT	Alt, Ref
ATP1A2	NM_000702.3	c.2493G>A:p.(Arg831Arg)	-	-	118, 110
CACNG2	NM_006078.4	c.256G>A:p.(Asp86Asn)	Benign	Deleterious	172, 179
		c.349A>G:p.(Met117Val)	Benign	Deleterious	148, 159
CHRNA4	NM_000744.6	c.505C>T:p.(Pro169Ser)	Prob	Deleterious	205, 231
CNR1	NM_016083.4	c.786T>C:p.(Ser262Ser)	-	-	102, 92
CYP2C9	NM_000771.3	c.551A>T:p.(His184Leu)	Benign	Tolerated	91, 64
FBXW7	NM_033632.3	c.468A>C:p.(Gln156His)	Prob	Tolerated	14, 21
FKBP4	NM_002014.3	c.64G>A:p.(Gly22Arg)	Poss	Deleterious	25, 16
FLOT1	NM_005803.3	c.71T>G:p.(Val24Gly)	Prob	Deleterious	44, 53
GRIN2B	NM_000834.4	c.1263T>C:p.(Ser421Ser)	-	-	148, 152
		c.831C>G:p.(Leu277Leu)	-	-	71, 104
		c.346T>G:p.(Ser116Ala)	Prob	Tolerated	135, 148
GRM5	NM_001143831.2	c.3266T>C:p.(Val1089Ala)	Benign	Tolerated	55, 57
KIF5A	NM_004984.3	c.2079T>C:p.(Asp693Asp)	-	-	167, 199
NTRK1	NM_001012331.1	c.1395G>A:p.(Leu465Leu)	-	-	130, 124
P2RX4	NM_002560.2	c.427G>C:p.(Gly143Arg)	Prob	Deleterious	91, 89
POLG	NM_002693.2	c.47C>G:p.(Pro16Arg)	Benign	Tolerated	25, 39
		c.984A>T:p.(Gln328His)	Benign	Tolerated	159, 223
SCN1A	NM_001165963.2	c.3483A>G:p.(Ala1161Ala)	-	-	14, 14
		c.2301C>T:p.(Asp767Asp)	-	-	118, 154
SCN3A	NM_006922.3	c.1950C>A:p.(Cys650*)	-	-	64, 102
SCN9A	NM_002977.3	c.5052A>T:p.(Thr1684Thr)	-	-	265, 286
SLC2A1	NM_006516.2	c.43G>A:p.(Ala15Thr)	Benign	Deleterious	78, 99
SLC6A2	NM_001043.3	c.140C>G:p.(Ala47Gly)	Benign	Tolerated	91, 161
		c.1711A>C:p.(Ile571Leu)	Benign	Tolerated	75, 73
SLC6A3	NM_001044.4	c.1372C>T:p.(Leu458Phe)	Benign	Tolerated	15, 11
TH	NM_199292.2	c.1224G>T:p.(Gly408Gly)	-	-	176, 241
TRPV1	NM_080706.3	c.1867C>T:p.(Pro623Ser)	Benign	Tolerated	121, 155
ZNF767P	NR_027788.1	n.1263del	-	-	108, 198
		n.2781C>T	-	-	24, 8
Intronic			MaxEnt	NNSPLICE	Alt, Ref
CNR2	NM_001841.2	c.-45-9G>C	+14.8%	+32.9%	70, 60
P2RX4	NM_002560.2	c.1045-18A>T	+24.6%	+52.5%	129, 141

Abbreviations: Prob, probably damaging; Poss, possibly damaging; Alt, alternate allele; Ref, reference allele.

Rare variants were found in 102 genes. All but one (IFI30) of the 70 genes with novel variants also had rare variants. Three genes (CYP19A1, IL2, MTCO2) had no such variants, while another 5 (ADORA1, HINT1, HTR2A, OXT, TTR) had no variant in either the exonic or intronic regions. Five genes (ADRB2, HINT1, HLA-B, IL1B, and PRRT2) had only one such variant. The 2 genes with the highest number of rare variants were CACNA1B with 47 and POLG with 29. Both NTRK1 and SCN10A had 21 while DNM2, KIF5A, and NOTCH3 had 20 variants. The remaining genes had 2–19 rare variants.

All 91 patients had at least 3 novel or rare variants (inclusive of intronic variants), or an average of 7.5 each. The highest number was 26 (one patient), followed by 25 (one patient) and 24 (one patient). There were two patients with 23 and another two with 19 variants. Three patients had 18 and the remaining 81 (89.0%) had between 3 and 17 variants each.

Rare variants that were enriched in the study population

Among the identified rare exonic variants (frequencies <1%) that are not synonymous, 54 had frequencies that were statistically significantly (P-value ≤0.05) higher than the corresponding frequencies for East/South Asians in the Genome Aggregation Database (gnomAD). Two of the 54 were in-frame: a 3-nucleotide insertion and a 3-nucleotide deletion. Of the 52 missense variants, 21 were predicted by both Polyphen-2 and SIFT to have a significant consequence on the encoded proteins, while another 13 were predicted to have a damaging effect by one of the two programs (Table 4). Two of the exonic variants (NOTCH3 c.3141C>G and POLG c.2069C>T) were also putative splice variants. For intronic variants that were rare, there were 8 with higher frequencies than those in gnomAD, and all were predicted to affect splicing. After Bonferroni correction for multiple testing, statistically significant difference remained for one exonic (POLG c.125_127dupGGC:p.(Arg42dup); corrected P-value of 0.017) and one intronic variant (CYP3A5 c.433-1G>C; corrected P-value of 0.017).

Table 4

List of rare non-synonymous or intronic variants with frequencies significantly higher than expected data in gnomAD

Gene	GenBank ref	Variant	Alamut visual prediction		This study Alleles counts		Weighted gnomADaAlleles counts		Fisher’s exact test
Exonic			PolyPhen-2	SIFT	Alt	Ref	Alt	Ref	P-value
ABCB1	NM_000927.4	c.2222G>T:p.(Arg741Ile)	Benign	Tolerated	1	181	1	19,040	0.019
ADRB2	NM_000024.5	c.776G>A:p.(Arg259His)	Prob	Tolerated	1	181	2	19,054	0.028
ANKK1	NM_178510.1	c.2059G>A:p.(Ala687Thr)	Benign	Tolerated	1	181	4	18,638	0.047
		c.404A>C:p.(His135Pro)	Prob	Deleterious	2	180	9	20,132	0.004
CACNA1B	NM_000718.3	c.265A>G:p.(Lys89Glu)	Prob	Deleterious	1	181	1	17,530	0.020
CHRNA4	NM_000744.6	c.1087G>A:p.(Val363Met)	Benign	Deleterious	1	181	2	20,507	0.026
CNR1	NM_016083.4	c.919C>T:p.(Arg307Cys)	Poss	Deleterious	1	181	0	19,021	0.010
COMT	NM_000754.3	c.718G>A:p.(Glu240Lys)	Benign	Deleterious	2	180	5	19,048	0.002
CYP2C9	NM_000771.3	c.1004G>A:p.(Arg335Gln)	Prob	Tolerated	1	181	4	19,044	0.046
CYP2C19	NM_000769.2	c.518C>T p.(Ala173Val)	Poss	Deleterious	3	179	70	20,530	0.026
CYP3A4	NM_017460.5	c.1105A>G:p.(Ile369Val)	Benign	Tolerated	1	181	0	20,498	0.009
CYP3A5	NM_000777.4	c.160C>A:p.(Arg54Ser)	Benign	Tolerated	1	181	0	19,048	0.010
DNM2	NM_001005360.2	c.2293C>T:p.(Pro765Ser)	Benign	Deleterious	1	181	0	19,052	0.010
		c.316G>A:p.(Asp106Asn)	Prob	Deleterious	2	180	28	20,530	0.028
		c.958G>A:p.(Asp320Asn)	Benign	Tolerated	1	181	1	20,526	0.018
EPHX1	NM_000120.3	c.130G>C:p.(Glu44Gln)	Poss	Tolerated	2	180	32	20,174	0.037
GRIK4	NM_014619.4	c.1247C>T:p.(Thr416Ile)	Benign	Deleterious	1	181	4	19,050	0.046
GRIN2B	NM_000834.3	c.3421_3423delGAG:p.(Glu1141del)	-	-	1	181	0	19,054	0.010
		c.514G>A:p.(Val172Ile)	Prob	Deleterious	1	181	0	20,521	0.009
GRM1	NM_001278064.1	c.2630G>A:p.(Arg877Gln)	Prob	Deleterious	1	181	0	19,001	0.010
GRM5	NM_001143831.2	c.2584C>A:p.(Leu862Ile)	Poss	Deleterious	2	180	9	19,022	0.005
HTR1A	NM_000524.3	c.722G>A:p.(Arg241His)	Poss	Tolerated	1	181	0	19,035	0.010
HTR2C	NM_000868.3	c.1255A>G:p.(Thr419Ala)	Benign	Tolerated	1	181	0	15,075	0.012
IL6	NM_000600.3	c.477G>T:p.(Lys159Asn)	Prob	Tolerated	1	181	0	17,813	0.011
KCNQ3	NM_004519.3	c.2305C>T:p.(Pro769Ser)	Poss	Deleterious	1	181	3	20,507	0.035
KIF5A	NM_004984.2	c.1995C>G:p.(Ser665Arg)	Benign	Tolerated	1	181	NAb	NAb	-
MYPN	NM_032578.3	c.2093A>G:p.(Asn698Ser)	Benign	Tolerated	2	180	37	20,518	0.046
NOTCH3	NM_000435.2	c.3141C>G:p.(Ile1047Met)	Benign	Tolerated	1	181	1	12,574	0.028
		c.515G>A:p.(Gly172Asp)	Prob	Deleterious	1	181	3	18,765	0.038
OXTR	NM_000916.3	c.490T>G:p.(Cys164Gly)	Benign	Tolerated	1	181	1	18,352	0.020
P2RX4	NM_002560.2	c.842C>T:p.(Thr281Ile)	Poss	Deleterious	1	181	1	19,050	0.020
P2RX7	NM_002562.5	c.556G>A:p.(Glu186Lys)	Prob	Deleterious	1	181	1	19,052	0.020
PDHA2	NM_005390.4	c.1082A>G:p.(Glu361Gly)	Poss	Tolerated	1	181	3	20,530	0.035
POLG	NM_002693.2	c.125_127dupGGC:p.(Arg42dup)	-	-	2	180	1	17,332	0.000c
		c.1402A>G:p.(Asn468Asp)	Benign	Tolerated	1	181	2	20,526	0.026
		c.1898A>C:p.(Lys633Thr)	Benign	Tolerated	2	180	6	19,031	0.002
		c.2069C>T:p.(Thr690Met)	Benign	Deleterious	1	181	3	20,530	0.035
		c.3139C>T:p.(Arg1047Trp)	Prob	Deleterious	1	181	3	20,526	0.035
SCN10A	NM_006514.3	c.2972C>T:p.(Pro991Leu)	Prob	Deleterious	1	181	0	20,525	0.009
		c.4417G>A:p.(Val1473Met)	Prob	Deleterious	1	181	0	18,785	0.010
		c.4766C>T:p.(Ala1589Val)	Prob	Deleterious	1	181	2	19,045	0.028
		c.5089G>A:p.(Val1697Ile)	Benign	Tolerated	3	179	42	20,520	0.007
SCN11A	NM_014139.2	c.2804A>C:p.(Gln935Pro)	Benign	Tolerated	1	181	2	19,029	0.028
SCN1A	NM_001165963.1	c.3283T>C:p.(Tyr1095His)	Prob	Deleterious	1	181	4	20,522	0.043
		c.4834G>A:p.(Val1612Ile)	Poss	Deleterious	3	179	44	20,256	0.008
SCN9A	NM_002977.3	c.554G>A:p.(Arg185His)	Prob	Deleterious	5	177	77	20,123	0.001
SLC2A1	NM_006516.2	c.322G>A:p.(Val108Met)	Poss	Deleterious	1	181	0	20,475	0.009
SLC6A2	NM_001043.3	c.730G>A:p.(Val244Ile)	Poss	Deleterious	1	181	0	20,530	0.009
TAGAP	NM_054114.4	c.1747C>A:p.(Gln583Lys)	Benign	Tolerated	1	181	0	19,050	0.010
		c.1907C>A:p.(Pro636His)	Prob	Tolerated	2	180	28	20,524	0.028
TH	NM_199292.2	c.770C>A:p.(Ala257Asp)	Benign	Tolerated	1	181	1	14,018	0.026
UGT2B15	NM_001076.3	c.28C>G:p.(Leu10Val)	Benign	Deleterious	1	181	1	17,200	0.021
		c.1058G>A:p.(Arg353Gln)	Benign	Tolerated	1	181	0	19,052	0.010
		c.1553G>A:p.(Arg518Gln)	Poss	Tolerated	1	181	2	20,523	0.026
Notes: aWeighted gnomAD frequencies of 94.6% East Asian (EAS) and 5.4% South Asian (SAS) populations. bAllele counts not available for East Asians or South Asians. cSignificant after Bonferroni correction.Abbreviations: Prob, probably damaging; Poss, possibly damaging; Alt, alternate allele; Ref, reference allele.
Intronic (splice variants only)			MaxEnt	NNSPLICE	Alt	Ref	Alt	Ref	P-value
ATP1A2	NM_000702.3	c.496-14G>C	+12.8%	+7.0%	1	181	2	18,182	0.029
CYP3A5	NM_000777.4	c.433-1G>C	−100.0%	−100.0%	3	179	13	20,528	0.000a
DNM2	NM_001005360.2	c.1782-7C>A	−24.0%	−26.1%	1	181	0	20,488	0.009
KCNIP3	NM_013434.4	c.307-15G>A	−18.0%	−4.3%	1	181	3	19,039	0.037
NTRK1	NM_001012331.1	c.360-4G>A	−4.6%	+10.0%	1	181	0	19,054	0.010
POLG	NM_002693.2	c.1712+5G>A	−100.0%	−98.3%	1	181	4	20,522	0.043
SCN3A	NM_006922.3	c.1032-3T>C	+12.1%	−1.8%	1	181	0	19,048	0.010
SLC6A4	NM_001045.5	c.1651-4T>C	−4.3%	−20.9%	1	181	3	20,441	0.035

Note: aSignificant after Bonferroni correction.

Abbreviations: Alt, alternate allele; Ref, reference allele.

The 62 rare variants in Table 4 (comprising 54 exonic that are non-synonymous and eight intronic-splice variants) were identified from 46 patients; 18 of whom also had novel variants of functional significance (non-synonymous or splicing variants). The highest number per patient was 6 (one patient) while another patient had 5. There were 2 patients with 4 such variants and 6 patients with 3. The remaining 36 patients had either 1 or 2 rare variants while 10 patients had only novel variants. Thirty-five patients did not have any novel or rare variants of functional significance that had higher frequencies than the general population.

Pathway analysis

Twenty-one of the genes that had either novel or rare variants with higher frequencies in this high-pain population were found to be involved in 7 non-redundant pathways in the STRING database (P-values of <0.001 after correcting for false discovery rate). The significantly enriched pathways include neuroactive ligand-receptor interaction, dopaminergic synapse and cocaine addiction, metabolism of xenobiotics by cytochrome P450 and morphine addiction, serotonergic synapse, and bile secretion, all known to be pain related (Table 5).

Table 5

KEGG pathways identified for genes with novel or rare nonsynonymous or splice variants identified in the study population

Pathway	P-valuea	Matching genesb
Neuroactive ligand-receptor interaction	4.75e-19 ~ 0.000858	CHRNA4, P2RX4, P2RX7, GRIK4, SLC6A4, CACNG2, HTR1A, CNR1.
Dopaminergic synapse	6.76e-14 ~ 0.000764	SLC6A3, CACNG2, GRIN2B, COMT, SLC6A4, TH, CACNA1B
Cocaine addiction	8.38e-09 ~ 0.000287	SLC6A3, TH, GRIN2B, COMT, CNR1
Metabolism of xenobiotics by cytochrome P450	2.21e-22 ~ 6.86e-14	EPHX1, CYP2C9, CYP2C19, CYP3A4, CYP3A5, UGT2B15, UGT2B15
Morphine addiction	6.47e-10 ~ 0.000858	ADRB2, CACNA1B, HTR1A
Serotonergic synapse	2.34e-14 ~ 0.000751	SLC6A4, HTR1A, CACNA1B, UGT2B15
Bile secretion	2.35e-11 ~ 0.000136	ATP1A2, CYP2C9

Notes: aCorrected for false discovery rate. bNovel/rare nonsynonymous and splice variants have been identified in the matching genes.

Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes.

Association of morphine usage with the presence of novel and rare variants

The patients were further grouped based on whether they carried the novel (listed in Table 3) and/or rare exonic variants (listed in Table 4). Their morphine usage was further compared with those who only had common variants. Although there was statistically significant difference only for the 20-hr PCA morphine, the trend was similar across all time-points (Table 6). The group with novel variants used more morphine compared with the group carrying rare variants. This in turn resulted in higher mean morphine dosage than the group of 35 patients with only common variants. There was no statistically significant difference in terms of age, BMI and self-reported pain scores between the groups.

Table 6

Comparison of morphine usage for patients carrying variants of different frequencies

Variable/group	Novel (n=28)a	Rare (n=46)b	Common (n=35)	P-valuec	P-valued
PCA morphine @4 hrs (mg)
Mean (SD)	8.43 (5.51)	8.28 (4.74)	7.40 (4.27)	0.632	0.681
Median	7.00	7.00	7.00
Min, Max	0, 20	2, 20	1, 20
PCA morphine @8 hrs (mg)
Mean (SD)	15.75 (9.10)	15.04 (8.36)	12.20 (7.48)	0.181	0.186
Median	16.00	13.00	13.00
Min, Max	0, 37	2, 38	1, 28
PCA morphine @12 hrs (mg)
Mean (SD)	20.32 (10.60)	18.46 (10.45)	14.43 (8.61)	0.054	0.074
Median	18.50	16.00	15.00
Min, Max	0, 42	3, 51	1, 32
PCA morphine @16 hrs (mg)
Mean (SD)	22.25 (12.00)	21.17 (11.53)	16.34 (10.11)	0.074	0.079
Median	20.00	20.50	15.00
Min, Max	0, 49	3, 52	1, 38
PCA morphine @20 hrs (mg)
Mean (SD)	25.04 (13.91)	24.46 (13.65)	18.03 (11.39)	0.049	0.034
Median	22.50	21.50	18.00
Min, Max	0, 52	3, 61	1, 41
PCA morphine @24 hrs (mg)
Mean (SD)	27.68 (15.24)	26.72 (14.61)	21.06 (12.94)	0.121	0.075
Median	24.50	24.00	20.00
Min, Max	0, 56	6, 61	1, 45
PCA morphine (total in mg)
Mean (SD)	28.30 (14.97)	26.98 (14.85)	21.29 (13.53)	0.119	0.066
Median	24.00	24.00	19.50
Min, Max	6,56	6,61	1,45
PCA morphine adjusted to body weight (mg/kg)
Mean (SD)	0.464 (0.245)	0.452 (0.255)	0.352 (0.218)	0.119	0.051
Median	0.377	0.412	0.330
Min, Max	0.097, 0.935	0.073, 0.963	0.017, 0.783

Notes: aTotal number of patients carrying the novel variants listed in Table 3. bTotal number of patients carrying the rare variants listed in Table 4 (including 18 who also had the novel variants in Table 3). ANOVA with Tukey post hoc tests for comparison between the three groups. dANOVA between the group with common variants (n=35) and the group carrying either novel and/or rare variants (n=56). Bold values indicate statistically significant.

Discussion

The advent of high-throughput genotyping technologies has led to the identification of genetic variants associated with many complex diseases and traits. In particular, GWAS had uncovered many common variants associated with various phenotypes. However, it is not designed to detect association involving variants of very low frequencies. Since NGS has become more cost-efficient, it is now feasible to genotype by resequencing, thereby uncovering the rare variants that may be important. By resequencing 105 known genes related to pain in our cohort of high postoperative pain patients, we were able to detect variants that were either absent or reported at very low frequencies in the general population.

The most interesting novel variant was the stop-gain in SCN3A. Pathogenic mutations in this gene have been linked to focal epilepsy. However, there was no record of this condition in the patient. On the other hand, there were 10 synonymous variants that were novel. Although synonymous variants are generally well tolerated and most have no functional consequence, novel variants are important in the context of rate of replication error, position within the gene and base changes tolerated. In addition, not all synonymous variants are insignificant as they might change the secondary structure of the mRNAs and influence their stability. In cases where the synonymous substitution affects RNA–RNA interaction, it may alter translation efficiency which is important in the biological regulation of gene expression and transcriptome complexity. Some synonymous substitutions are also known to lead to aberrant splicing.33–35

Among the genes with novel and rare variants that had higher frequencies in this sample set, 21 were mapped to the pain-related KEGG pathways. These genes likely play key roles in processes related to pain response and signaling, as well as drug metabolism. Recent animal studies have confirmed that neuroactive ligand-receptor influenced disease-related pain and its severity, duration, and relief.36,37 It is not surprising that variants in the cytochrome P450 pathway are more common in this high-pain population as opioids are metabolized through two major enzyme systems: CYP450 and, to a lesser extent, UDP-glucuronosyltransferases. Although more than 50 CYP450 enzymes are known, CYP1A2, CYP2C9, CYP2D6, CYP3A4, and CYP3A5 account for the metabolism of up to 90% of drugs.38,39 Increased cytochrome P450 enzyme activities may result in faster metabolism and shorter half-life, which may potentially decrease a drug’s pharmacologic effect.

Five of the genes in Table 5 are involved in the cocaine addiction pathway. Based on experimental studies and some indirect clinical evidence, dopamine has been suggested to have anti-nociceptive effect.40–42 Cocaine increases the level of dopamine and cocaine addiction is related to pain.43 In addition, serotonergic synapse44 and bile secretion45 pathways are also linked to pain. Variants in any of these genes might have affected the functions of the encoded proteins, resulting in the enhanced and prolonged postoperative pain experienced by our study subjects.

Genes with novel or rare variants that had higher frequencies than the general population include POLG, SCN10A, and DNM2. POLG (DNA polymerase gamma) encodes a polymerase responsible for the replication of human mitochondrial DNA. Mutations in the gene have been linked to mitochondrial diseases, such as myocerebrohepatopathy spectrum disorders, Alpers-Huttenlocher syndrome, myoclonic epilepsy myopathy sensory ataxia, ataxia neuropathy spectrum, and progressive external ophthalmoplegia.46,47 POLG mutations are associated with peripheral neuropathy and a potentially painful, axonal/mixed, mainly sensory polyneuropathy48 and muscle pain.49 This gene had the highest number of identified variants (six rare and two novel), including one missense variant found in 10 patients. Another gene SCN10A had four rare variants which were more prevalent in the study population. It encodes a component of the Nav1.8 sodium channel and is associated with peripheral neuropathy.50 The other gene which had multiple variants with higher frequencies is DNM2. This gene codes for Dynamin-2, one of the subfamilies of GTP-binding proteins. DNM2 has been associated with pain flare in patients who received palliative radiation therapy for painful bone metastases.51

Although our study uncovered novel and rare variants from patients who reported higher pain and used more morphine, it has several limitations. First, sequencing was only performed in <10% of a patient cohort, on those with the highest pain burden (self-reported pain scores and high morphine use). Second, the frequency comparison was done with data from population databases. In addition, functional effects were based on in silico predictions, and no in vivo or in vitro studies were carried out for validation. Lastly, although the number of reads for reference and alternate alleles were similar, the variants were also not Sanger validated. Therefore, further studies are warranted to address these limitations.

In summary, our results showed that some rare variants were more common in patients who reported more pain and used more PCA morphine. We also identified several novel variants that were predicted to either result in amino acid substitutions or affect splicing. Carriers of such variants tend to use more morphine over the first 24 hrs of the postoperative period. Whether the novel variants affect the sensitivity and tolerance to pain remain to be investigated. The cost of genomic technologies has become more affordable, and the analysis of sequencing data is also amenable to automated pipelines. Thus, it is possible to incorporate genotyping or sequencing for a set of gene variants that account for a significant portion of the inter-individual variation. The genetic information could be combined with other predictive factors in patient risk stratification. This will enable early intervention and timely modulation of nociception that has been shown to reduce the incidence of persistent pain and improve patient recovery.

Supplementary materials