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

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.¹

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 pain¹⁶ and as high as 70% for low back pain.¹⁷

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)²¹ and the "Human Pain Genetics Database".²² 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.²⁷ 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

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

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.²⁸ 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”.²⁹ For PolyPhen‐2, scores of >0.85 were classified as “probably damaging”, and scores of 0.15–0.85 were considered as “possibly damaging”.³⁰

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.³¹

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/)³² 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).

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×.

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

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

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

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

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.⁴³ In addition, serotonergic synapse⁴⁴ and bile secretion⁴⁵ 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 polyneuropathy⁴⁸ and muscle pain.⁴⁹ 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.⁵⁰ 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.⁵¹

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.