Prerpints.org logo
Preprint
Article

Genetic Variants Linked to Opioid Addiction: A Genome Wide Association Study

This version is not peer-reviewed.

Submitted:

05 October 2024

Posted:

07 October 2024

You are already at the latest version

A peer-reviewed article of this preprint also exists.

Abstract
Opioid addiction disorder (OAD) affects millions of people worldwide. While it is known that OAD originates from many factors, including social and environmental factors, the role of genetic variants in developing the disease has also been reported. This study aims to investigate the genetic variants that may be associated with the risk of developing OAD upon exposure. Twenty-three subjects who had previously been given opioid-based painkillers to undergo small surgical treatment were recruited at Prisma Health Upstate clinic and elsewhere. Eleven of them were considered non-persistent opioid users (controls), and 12 of them were persistent opioid users (cases) after an initial surgery and at the time of sample collection. The subjects were asked to provide saliva samples, which were subjected to DNA sequencing at Clemson University Center for Human Genetics, and variant calling was performed. The GWAS for genes known to be associated with OAD resulted in 13 variants (intronic or SNV) with genome-wide significance (raw p-value < 0.01) and two missense variants, rs6265 (Val66Met in BNDF isoform a) and rs1799971 (Asn40Asp) in OPRM1 previously reported in the literature. Furthermore, extending the GWAS to find all the genomic variants and filtering the variants to include only variants found in cases (persistent opioid users) but not in controls (non-persistent opioid users) resulted in 11 new variants (p-value< 0.005). Considering that OAD is a complex disease and the effect may come from different variants in the same genes, we performed a co-occurrence analysis of variants on the genes and identified four genes, LRFN3, ZMIZ1, RYR3, and OR1L6 with three or more variants in the case subjects but not in control individuals.
Keywords: 
Subject: 
Medicine and Pharmacology  -   Anesthesiology and Pain Medicine

1. Introduction

Opioid addiction disorder (OAD) affects millions of people worldwide, and it has immense and multifaceted adverse socio-economic impacts affecting individuals, families, communities and economies [1,2]. The impacts involve drug overdose deaths [3], increased healthcare costs [4], low productivity [5], heightened criminal justice involvement [6,7], and broader societal challenges. Human studies consistently show that the heritability of opioid addiction is substantial [8,9]. Twin studies estimate the heritability of opioid use disorder (OAD) to be between 40% and 60% [8,9]. Family studies also suggest that first-degree relatives of individuals with OAD have a higher risk of developing the disorder [10]. Identifying genetic variants and genes linked to elevated risks can be used to develop a better understanding of the biology of OAD. It can be further utilized in screening people for their genetic susceptibility to OAD and to avoid prescription of opioid-based painkillers to limit the risks of developing OAD. The two popular approaches to association studies are (a) Candidate gene studies, which are hypothesis-driven studies that analyze genes with known connections to a phenotype, and (b) Genome-wide association studies (GWAS), which are statistical analyses of polymorphisms across a genome for association studies. Both approaches have been used to study genetic association with OAD in the past [11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26].
Several candidate gene studies aim to identify the genetic variants related to OAD [11,12,13,14,15,16,17,18,19,20,21] have reported several genetic variants on genes. For example, the mu-opioid receptor gene (OPRM1), and the variant rs1799971 (A118G), have been associated with differences in opioid effects and addiction risk [14,15,20]. Studies show that individuals with the G allele may experience altered drug effects and a higher risk for OAD [27]. Likewise, the dopamine receptor D2 gene (DRD2) has been implicated in the reward system’s role in addiction. Variants like the Taq1A polymorphism are associated with an increased risk of addiction due to their impact on dopamine receptor availability and function [11,12,18,19]. The delta-opioid receptor (DOR), encoded by the OPRD1 gene, has been investigated in the context of opioid addiction, though it is less studied compared to the mu-opioid receptor [16]. Research indicates that DOR plays a role in modulating mood, anxiety, and stress responses, all of which are factors in addiction vulnerability. Variants in the OPRD1 gene have been linked to altered responses to opioids, suggesting that DOR could influence individual susceptibility to OAD [16,17]. Brain-derived neurotrophic factor (BDNF) is a neurotrophin critical for brain plasticity, learning, and memory [28]. Alterations in BDNF expression and function have been associated with various neuropsychiatric conditions, including addiction [29]. Studies suggest that BDNF plays a role in the brain’s reward system and its response to opioids. Genetic variations in the BDNF gene may influence the risk of developing OAD by affecting neural plasticity and resilience to stress [13,19].
A handful of case-control genome-wide association studies (GWAS) have also been reported using varied cohorts consisting of populations of varied ethnicity, sex, and whether the controls were exposed to opioids [22,23,24,25,26]. For example, Nielsen DA et al. performed a GWAS on subjects of European origin with 104 heroin-dependent patients and 101 controls [22]. They compared the frequencies of 10,000 SNPs but failed to find any significant associations. However, they extended the cohort consisting of subjects of African and European descent and a larger pool of 100,000 SNPs in methadone-maintained 325 heroin addicts and 250 controls [23]. After the multiple testing corrections, only a single intergenic variant, rs10494334, was found significant in European subjects, while none in subjects of African descent.
Gelernter J et al. reported a GWAS on 5432 African American and 6877 European American subjects [24]. Analysis using the Diagnostic and Statistical Manual of Mental Disorders 4th edition (DSM-IV) symptom counts for opioid dependence or case-control status was performed on sub-groups and meta-analysis on the entire cohort. In the final meta-analysis, a variant rs62103177 in KCNG2 was found to be genome-wide significant. Other variants as rs60349741 in KCNC1 and rs114070671 in APBB2 had genome-wide significance in combined analysis but only nominally in meta-analysis [24].
Defining controls merely based on the DSM definition has a caveat that the subjects not exposed to opioids cannot develop OAD, even though they carry the risk-associated variants. In a study of European Americans consisting of 1290 cases and 1768 opioid-exposed controls, an SNP, rs12442183, 110 kb downstream to RGMA, was found to be significantly associated with OAD [25]. Further, microarray data analysis suggested that rs12442183 is an expression quantitative trait locus (eQTL) for the RGMA gene. However, previously reported variants on KCNG2, KCNC1, and APBB2 didn’t have a significant association despite the study using the overlapping sample with the previous study [25,30].
In another study comparing the daily injection opioid users (cases: n=1167) vs. the opioid abusers who never injected opioids daily (controls: n=161), a variant rs1436175 in the gene CNIH3 showed genome-wide significance. Five other variants, rs10799590, rs12130499, rs298733, rs1436171, and rs1369846, in CNIH3, reached genome-wide significance in a meta-analysis of the discovery cohort and two independent populations [26,30].
The studies mentioned above [22,23,24,25,26], consolidated together, indicate that clear external replication of GWAS findings is rare. The following factors contribute to this situation: (a) OAD is a complex psychiatric disease with relatively low heritability, and there is no single variant with a large effect size that can be detected in small cohorts (b) previous OAD GWAS were relatively small compared with those for legal substance use disorders and (c) in published work relevant to opioid use, there was considerable phenotypic heterogeneity across samples [30], and (d) OAD is likely highly polygenic and individual loci contributions are too small for detection in studies with limited sample size.
With the above issues in mind, we present a case-control study of twenty-three individuals who had been previously given opioid-based painkillers to undergo small surgical treatment and were recruited at Prisma Health Upstate clinic and elsewhere. Eleven of them were considered non-persistent opioid (controls), and twelve of them were persistent opioid users (cases) after an initial surgery and at the time of the sample collection (details in Section 4.1). The subjects were asked to provide saliva samples, which were subjected to DNA sequencing at Clemson University Center for Human Genetics, and variant calling was performed (details in Section 4.2). The variants are analyzed to identify their association with OAD.

2. Results

The first task of our study is to check if previously reported variants can be found in our cohort. For this purpose, from the literature, we compiled a list of genes harboring variants reported to have associations with OAD (Table 1). Then, we present the results of variants found on the listed genes in our present study. Two of the variants with coding consequence (rs6265 in BDNF and rs1799971 in OPRM1) show some indication of the association to OAD, though the p-values are not small enough (raw p-values are < 0.3) to call these significant. Then the search is extended to variants with no direct coding consequences on these genes, and several such variants on DRD3, KCNG2, and NRXN3 are found with (significance p-values < 0.01) for the association with OAD (Section 2.1). The next task is to search for genetic variants across the whole genome, where the alternate alleles are exclusively found in cases but not in any of the control subjects (Section 2.2). We hypothesize that several variants on the same gene may be working together to compromise the function of a gene, and to test it, we performed a co-occurrence analysis of variants on the genes. In co-occurrence analysis, multiple variants on the same genes are assumed to have additive impairments on the function of the gene, and these variants are grouped and analyzed for their functional consequences. This analysis identified LRFN3, ZMIZ1, RYR3, and OR1L6 with three or more variants on them and alternate alleles frequency sum greater or equal to 14 (Section 2.2.1), which potentially may have an association with OAD.
The OAD is a complex phenotype, and it may not result from a single variant or single gene, but a network of genes-/protein-interactions is perturbed due to several variants on multiple genes. To investigate it, we also performed the network analysis of protein-protein interaction (PPI) and network analysis of the genes harboring single or multiple variants suggesting an association with OAD. We identified a large, connected component of PPIs involving more than half of the genes of interest (GOI) (Section 2.3). Finally, we perform a gene function enrichment and pathway analysis to identify the molecular processes, biological function and interaction pathways that may be affected due to the presence of variants in the cases but not in controls and discuss the results (Section 2.4).

2.1. Variants on Genes Reported in the Literature with Association to OAD

We begin our association analysis to OAD with the set of genes reported in the literature to be associated with OAD. Thus, we analyzed gene variants (Table 1) to test their association with OAD in our samples. We extracted the variants within extended gene loci, for each gene, an extended region with 5 kb upstream and 5 kb downstream, including the gene boundary, is defined. The effects of variants are predicted using the Variant Effect Predictor (VEP) version 111 [38] about human genome assembly GRCh38.
A total of 6903 variants, including 67 novel, variants effects were predicted, consisting of one stop_gained, one inframe_insertion, one inframe_deletion, 58 missence_variant and 56 synonymous_variant with coding consequences. Further association testing is performed using PLINK v1.90b7.2 [39] of these variants. However, none of the variants with coding consequences could reach the genome-wide association significance (raw p-value < 0.01). In contrast, 13 variants (intronic or SNV) showed genome-wide significance (Table 2). The validation of intronic variants’ association with OAD requires further study. Considering this, we focused on only variants with the most severe consequence (coding consequences), but none showed significance for genome-wide association to OAD. However, two variants, rs6265 (C>T at 11:27658369) on BDNF, a missense coding variant (Val66Met in BNDF isoform a, previously reported [19]), and rs1799971, a missense variant on OPRM1 (118A>G) changing Asn40Asp in mu-opioid receptor (also previously reported [14,15,20]), have a raw p-value < 0.30, after relaxing the p-value to call association. Here, we understand that a raw p-value of 0.30 is not significant to ascertain an association. However, this very weak indication aligns with the previously reported studies. The large p-value of these variants’ association with OAD can be attributed to the small cohort used here.
The variants with no direct coding consequence found significant in extended gene loci of genes listed in Table 1 are shown in Table 2. These variants may have an impact on the genes by altering the expression of the genes. We found some variants, though not the listed variants, in genes DRD3 [32,33], KCNG2 [24,31], and NRXN3 have been reported in the literature to have an association with OAD [36,37].

2.2. Genes with Alternate Allele Exclusively in Cases but Not in Controls

We extended our search to include all the genomic variants found in the samples. The variant effect is predicted using VEP version 111 [38] for 20,910 variants found in the sample. Afterward, all the variants are ranked based on the severity of the predicted consequence and filtered to satisfy two requirements: (a) only variants for which all controls have homozygous reference alleles, and (b) only variants in the OAD group have the alternative alleles. A variant’s association with phenotype OAD is considered significant when the raw p-value is less than 0.05. It resulted in 158 variants, with significance within 126 genes. In the case of these variants with significance, the alternate allele frequency ranged from 11 to 4 in case samples and zero in controls. A list of such variants with association to OAD with only high-significance (raw p-value < 0.005) is provided in Table 3.
So far none of the genes listed in Table 3 has been directly linked to addiction. However, in a study of attention-deficit hyperactive patients, gene MTUS2 showed the highest variation frequency [40]. Studies have suggested the association of MTUS2 with psychological disorders [40]. The suggested association of MTUS2 to the psychological disorder may indirectly relate to addiction-like traits, but further exploration in this regard is needed to test this possibility.

2.2.1. Co-Occurrence of Variants on Genes and Association to OAD

So far, our focus has been on finding common genetic variants that may be associated with OAD. However, the OAD phenotype may result from the altered function of one or more proteins, and alteration could be caused by different variants within the same proteins, not necessarily identical variants in the affected individuals. It is also possible that there are multiple mechanisms and disruptions of pathways that can result in OAD. This led us to focus on proteins or genes and their combinations, instead of individual variants. Considering it, we analyzed the co-occurrence of multiple variants within the same gene/protein. The alternate alleles count of all gene variants exclusively in cases is aggregated by summing it, among the total 126 genes harboring variants associated with OAD, considering the significance (raw p-value < 0.05). Eighteen genes have multiple variants exclusively in cases, including 12 genes with two variants, three genes with three variants and three genes with greater than three variants.
The following genes harbor the most severe consequence (transcript ablation, splice acceptor variant, splice donor variant, stop gained/lost, frameshift variant, start lost, transcript amplification, feature elongation, or feature truncation) variants with alternate allele frequencies sum greater or equal to 10, ZMIZ1, LRFN3, OR1L6, RYR3, PWWP2B, ZNF92, CYP4F12, FAM186A, SCUBE2, NUTM2D. The number of variants, count of alternate alleles in cases and list of predicted consequences for the variants are provided in Table 4. Considering the significant raw p-values of the genes listed in Table 3 or Table 4 consisting a total 18 genes is compiled, we will call these as genes of interest (GOI) hereon for further analysis.
So far, among the genes/proteins listed in Table 4, none has been directly associated with addiction. However, some of these have been known to be associated with functions relevant to OAD. For example, LRFN3 (Leucine Rich Repeat And Fibronectin Type III Domain Containing 3), also known as SALM4 (Synaptic Adhesion-Like Molecule 4), is part of the LRFN family, which is involved in synaptic adhesion and regulation of excitatory synapses [41]. SALMs, including LRFN3/SALM4, play important roles in synaptic development and plasticity [42], which are critical processes in learning, memory, and behavior—all relevant to addiction. Similarly, RYR3 (Ryanodine Receptor 3) is part of the ryanodine receptor family, which functions as intracellular calcium channel that releases calcium from the endoplasmic reticulum into the cytoplasm [43,44]. Calcium signaling is critical for various cellular processes, including neuronal activity, synaptic plasticity, and neurotransmitter release [45]. ZMIZ1 (Zinc Finger MIZ-Type Containing 1) is a transcriptional coactivator that interacts with the androgen receptor and other transcription factors [46], and the androgen receptor has been implicated in reward-related behaviors, including aggression and stress responses, which can influence addiction [47]. In summary, the genes listed in Table 4 can potentially influence addictive behaviors in direct or indirect ways and in individual- or cohort-specific manners. This led us to further analyze the interaction networks of these genes and investigate their collective behaviors and plausible linkage to opioid addiction.

2.3. PPI Network Analysis of Genes of Interest

To gain insight into the functions involving the proteins in GOI, we built a protein-protein interaction (PPI) network based on the STRING database [48] as detailed in Section 4.4. We observe a total of seven connected components. Connected components 2 through 7 are small and shown in Figure 1; all other genes form a large connected component, marked as component-1. Component-1 consists of six genes listed in Table 3 and five genes harboring multiple variants listed in Table 4. This indicates that genes in component-1 are functionally interdependent, and any perturbation to this network of PPIs can potentially be associated with altered function and, thereby, some phenotype. To further investigate this hypothesis of perturbation of the PPI network and association with OAD, we performed a functional enrichment analysis of the genes in the network.

2.4. Gene Ontology and Functional Enrichment Analysis

The functional enrichment analysis was performed using the g:Profiler’s webservice[49], and its results are summarized in Figure 3 and Table 5.
Preprints 120336 g002
Figure 2. Summary of Gene Ontology and functional enrichment of Genes in PPIs network of GOIs performed using g:Profile webserver. Highly significant (considering the adjusted p-value), but non-redundant terms in the enrichment are annotated. The description of annotation labels to terms is provided in Table 5.
Figure 2. Summary of Gene Ontology and functional enrichment of Genes in PPIs network of GOIs performed using g:Profile webserver. Highly significant (considering the adjusted p-value), but non-redundant terms in the enrichment are annotated. The description of annotation labels to terms is provided in Table 5.
Preprints 120336 g002
Preprints 120336 g003
Figure 3. The distribution of percent of total PPIs as a function of scores: (a) combined_score and (b) Normalized score(norm_score) listed in STRING database for humans. The kernel density of the percent of PPIs is shown in a logarithmic scale to highlight distributions in low percent regions (in inset) within respective plots. The thresholds in this work are marked with dashed-vertical lines in both cases.
Figure 3. The distribution of percent of total PPIs as a function of scores: (a) combined_score and (b) Normalized score(norm_score) listed in STRING database for humans. The kernel density of the percent of PPIs is shown in a logarithmic scale to highlight distributions in low percent regions (in inset) within respective plots. The thresholds in this work are marked with dashed-vertical lines in both cases.
Preprints 120336 g003
From the enrichment analysis, the calcium transport, regulation and binding, enzyme binding, transmembrane transporter binding, and signaling receptor binding emerge as the high-confidence molecular function from gene ontology. While, signaling and amino-acid modification as the top biological function. The GO:cellular-component is the membrane and channels, which host most of the signaling and transmembrane proteins to carry the signaling in cells.
The KEGG pathway analysis shows enrichment in Calcium signaling pathways, circadian entrainment, dopaminergic synapses, morphine addiction and alcoholism (Figure 2 and Table 5). This further supports that the variants of GOIs are impairing or abolishing the function of the proteins, which in turn compromises the protein’s normal function in the PPI network since these proteins are involved in crucial reward pathways and addiction. Therefore, the association of these genetic variants to OAD risk is highly plausible, and further study is needed to confirm this.

3. Discussion

Opioid addiction disorder (OAD) is a complex phenotype with significant heritability estimated to be between 40% and 60%. The heritability raises the genetic susceptibility of individuals carrying the associated variants in their genome. In the past, several candidate gene studies and genome-wide association studies have implicated the association of variants on genes to OAD. Though, they have rarely been replicated in other studies. It has been pointed out that this situation can be due to the following limitations of the studies: (a) the cohort sizes used in these studies are relatively small compared to those used for substance use addiction, (b) the DSM definitions have been used to label controls, which may not be appropriate and thus diluting the association results, and (c) the ethnicity of populations of cohorts varied in different studies, which brings another variability in the study. Another limitation of some of the previous studies is that the controls were not exposed to substances of interest. Such controls could potentially carry genetic variants, making them predisposed to OAD, but they have never been exposed to drugs.
To address the abovementioned deficiency, We performed the genome-wide association study of a cohort of 23 samples, including 12 cases and 11 controls, where all controls were given opioid-based pain killers. In our study, we could find a weak association for two variants, one on BDNF and the other on OPRM1, when our focus was on variants only with coding consequences. Upon extending the search to variants with no direct coding consequence, we identified variants of DRD3, KCNG2, and NRXN3 genes with significant association with OAD. Interestingly, these genes have been predicted to have an association with OAD. This indicates that these variants may be altering the expression of the protein and thereby indirectly affecting the phenotype OAD. Next, we expanded our search space for the variants to the entire genome. However, genome-wide search yields many variants. Thus, only variants with coding consequences are prioritized in association analysis. We chose only variants with alternate alleles in case subjects and predicted them to have an association with OAD with significance. Several variants of the same gene protein may have a cumulative effect on the gene’s function, thus altering its biological function. We used co-occurrence analysis to test this hypothesis. Our co-occurrence analysis identified LRFN3, RYR3, ZMIZ1, and OR1L6 as genes harboring multiple variants exclusively in cases. From the genes with predicted association to OAD with significance (genes of interest/GOI), a network of protein-protein interaction is extracted based on PPIs listed in the STRING database for humans. The network is further expanded by including genes directly interacting with two or more genes in the GOI list. Then, the pathway and functional enrichment of the genes in the network are done to find the pathways with significant enrichment for the genes and their association with the OAD phenotype. We observe considerable enrichment for terms like calcium signaling pathway, circadian entrainment, dopaminergic synapse, oxytocin signaling pathway, estrogen signaling pathway, morphine addiction, alcoholism, and opioid signaling, which are closely related to OAD or addiction in general.

4. Materials and Methods

4.1. Sample Collection

Local Institutional Review Board approval was obtained before subjects recruitment and enrollment, and written informed consent was obtained from all eligible participants. Eligible subjects were those over 18 who had undergone treatment for a small surgery in the past and provided informed consent. Subjects were excluded if they had an active oral lesion, had a history of opioid abuse before treatment of their injury, or had a history of opioid use within one year of treatment for their injury. Subjects who provided informed consent were subsequently asked to provide a saliva sample for DNA analysis. Twenty subjects were recruited and enrolled between November 2022 and June 2023. Twelve subjects were persistent opioid users at the time of enrollment and saliva sample collection, whereas the non-persistent opioid users were only eight. To make the set balanced, we collected additional samples from researchers from the corresponding author (E.A.) department who, in the past, have been sedated for small surgery but are not persistent users. Thus, the number of cases of non-persistent users is eleven, very similar to the number of cases of persistent users.
Age, gender, and race demographics were recorded at the enrollment visit, in addition to details of their orthopedic injuries. The average age of all enrolled participants was 66 years, with a minimum age of 36 and a maximum age of 84. There were 8 females and 3 males in the non-persistent opioid user group versus 7 females and 5 males in the persistent opioid user group.

4.2. DNA Sequencing and Genotyping

The samples’ DNA was extracted using a modified Zymo Quick-DNA Microprep Plus Kit protocol. Libraries were constructed using 500 ng of extracted DNA with the Tecan Revelo DNA-Seq Enz kit on the Tecan MagicPrep NGS system. The libraries were quantified for concentration with the Invitrogen 1x dsDNA HS assay kit on the Invitrogen Qubit 4 Fluorometer. Library size was quantified using Agilent High Sensitivity D1000 assay on the 4150 TapeStation system. Libraries were normalized to 3 nm, pooled, and sequenced on an Illumina Novaseq 6000 sequencer using an S4 flow cell and 2x150 sequencing chemistry at Clemson University Center for Human Genetics.
Raw reads were filtered for low-quality and short reads, then aligned to human reference genome version GRCh38 to produce alignment files using the GPU-accelerated fq2bam module in the NVIDIA Clara Parabricks suite. Base quality recalibration was performed using GPU-accelerated Genome Analysis Tool Kit [50] BQSR module from Parabricks suite and known variant information from Mills and 1000 Genomes Gold Standard Indel dataset and dbSNP v138 dataset. Variant calling was performed on the recalibrated alignment files using GPU-accelerated GATK’s haplotypecaller module from the Parabricks suite and GRCh38 reference genome. Individual sample GVCFs were combined, indexed and joint-called using the Short Variant Discovery workflow from GATK’s Best Practices.[51] Joint-called variants were hard-filtered using gold standard default values recommended by the Broad Institute [52]. After that, the effects for variants are predicted using the Variant Effect Predictor (VEP) version 111 [38] compared to standard human genome assembly GRCh38.

4.3. Variants Filtering, Effects Prediction and Association Analysis

A subset of all variants found after variant calling is extracted into a vcf file using the bcftools module in SAMtools [53] for analyzing the association of a subset of genes listed in Table 1. The subset .vcf file is indexed and compressed using bcftools [53]. The variant effects are predicted using Ensembl Variant Effect (VEP) version 111 [38] with reference to human genome assembly GRCh38. We used PLINK v1.90b7.2 [39] to test the significance of the genotype-phenotype association.

4.4. PPIs Network Construction and Analysis

The Protein-Protein Interactions data for humans is downloaded from STRING database [48]. Considering the large number of human PPIs listed in STRING and varying degrees of confidence in PPIs, we started filtering in only high-confidence PPIs based on absolute combined_score, including only interactions with combined_score greater than a threshold.
However, we noticed that filtering based on absolute combined_score excluded all the interactions of proteins that had no partner with combined_score greater than the chosen threshold. However, such filtering excluded some of the GOIs altogether. To avoid such stringent absolute combined_score based filtering, we used normalization of combined_score, carried as follows. The combined_score is normalized by the max combined_score for the gene among all its interacting genes. Such normalization is carried out to ensure that high-confidence PPIs can be included in a subnetwork featuring only high-confidence interactions Figure 3. The percentage of PPIs keeps decreasing as a function of combined_score and reaches a minimum of around 800; after that, it starts increasing, implying there are high-confidence PPIs pairs with scores greater than 800; based on this, we choose combined_score=800 as the threshold for unnormalized scores (Figure 3). We observe a similar trend for normalized scores as well, and around norm_score=0.8, again the percent PPIs start increasing; thus, we choose 0.8 as a threshold for the norm_score to find the filtered subset of the STRING database consisting of 19,622 nodes and 204,572 edges (after removing self-loop and duplicate edges) for further analysis using Cytoscape-v3.10.2 [54].
We start with creating a subnetwork of PPIs for GOIs. Initially, the genes listed in Table 3 or Table 4 and any immediate neighbor to these are also selected. Further, any node within a distance of one and has at least two neighbors in the selected nodes set is also considered part of the subnetwork to extend the network of PPIs of GOIs. Finally, the subnetwork is created using selected nodes and all the edges. The created subnetwork consisted of 146 nodes and 390 edges, as shown in Figure 1.

4.5. Gene Ontology and Functional Enrichment

The g:Profiler [49] carries the enrichment of the genes/proteins list provided as input against the data sources: Gene ontology (molecular function, biological process and cellular component), biological pathways (KEGG, Reactome and WikiPathways), regulatory motifs in DNA (TRANSFAC, miRTarBase), protein databases (Human Protein Atlas and CORUM), and human phenotype ontology (HP). We used the list of genes present in the connected component-1 of the network for the enrichment analysis. The gene ids are: Gn148400, Gn132535, Gn139549, Gn198838, Gn157322, Gn145687, Gn089250, Gn105426, Gn103254, Gn140368, Gn172824, Gn109971, Gn166444, Gn110031, Gn129007, Gn133835, Gn137486, Gn203740, Gn188612, Gn141480, Gn141738, Gn160691, Gn237172, Gn175470, Gn050820, Gn139197, Gn153485, Gn167378, Gn152229, Gn221886, Gn206069, Gn088832, Gn120899, Gn213965, Gn204256, Gn132589, Gn197037, Gn177885, Gn116030, Gn221874, Gn176533, Gn067365, Gn214562, Gn119782, Gn178372, Gn189350, Gn198363, Gn078369, Gn169717, Gn165059, Gn137312, Gn169083, Gn244462, Gn197122, Gn162521, Gn184292, Gn142875, Gn203791, Gn143318, Gn143801, Gn198626, Gn142949, Gn196218, Gn139160, Gn166321, Gn180900, Gn168214, Gn188322, Gn188100, Gn169398, Gn166407, Gn108175, Gn141404, Gn186439, Gn146757, Gn185522, Gn183310, Gn183715, Gn080815, Gn184716, Gn100077, Gn178363, Gn173020, Gn175356, Gn072062, Gn171806, Gn171813, Gn169925, Gn169885, Gn170606, Gn164690, Gn163501, Gn162928, Gn161533, Gn160014, Gn156521, Gn150787, Gn146700, Gn146648, Gn144401, Gn141736, Gn140876, Gn089159, Gn138798, Gn141867, Gn118729, Gn137727, Gn136943, Gn136159, Gn135617, Gn130950, Gn130726, Gn127947, Gn127588, Gn108788, Gn126243, Gn112874, Gn004468, Gn103653, where prefix Gn is used in listing here in place of ENSG00000 for brevity. The enrichment query parameters are as follows: version e111_eg58_p18_f463989d; date 9/30/2024, 8:18:46 PM; organism hsapiens; query length 122; all results false; ordered false; no iea false; sources GO:MF, GO:CC, GO:BP, KEGG, REAC, TF, MIRNA, HPA, CORUM, HP, WP; multiquery false; numeric ns ENTREZGENE_ACC; domain scope annotated; measure underrepresentation false; significance threshold method g_SCS; user threshold 0.05; no evidences false; highlight results true. The results are summarized and discussed in the result section.

Author Contributions

Conceptualization, S.K.P., V.S. and E.A.; methodology, S.K.P. and V.S.; investigation, S.K.P.; writing—–original draft preparation, S.K.P.; writing—–review and editing, S.K.P., E.A. and V.S.; supervision, E.A.; project administration, E.A.; funding acquisition, E.A. All authors have read and agreed to the published version of the manuscript.

Funding

This work of S.P. and E.A. was supported by grants from National Institutes of Health (United States), grant numbers R01GM093937 and R35GM151964 and P20GM139769 awarded to the Clemson Center for Human Genetics

Institutional Review Board Statement

The study of twenty patients from Prisma Health Upstate clinic was conducted in accordance with the Declaration of Helsinki, and approved by the Institutional Review Board (or Ethics Committee) of Prisma Health Upstate clinic (protocol code Pro00084425 approval date November 8, 2019).

Informed Consent Statement

Informed consent was obtained from all twenty subjects from Prisma Health Upstate clinic.

Data Availability Statement

The data is not available due to privacy reasons.

Acknowledgments

We acknowledge the Palmetto Cluster, Clemson University, for providing computational resources.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:
DSM-IV Diagnostic and Statistical Manual of Mental Disorders 4th edition
GWAS Genome-Wide Association Study
DSM Diagnostic and Statistical Manual of Mental Disorders
OAD Opioid Addiction Disorder
GOI Genes of Interest
SNP Single Nucleotide Polymorphism
PPI Protein-Protein Interaction
GO Gene Ontology

References

  1. Florence, C.; Luo, F.; Rice, K. The economic burden of opioid use disorder and fatal opioid overdose in the United States, 2017. Drug Alcohol Depend. 2021, 218. [Google Scholar] [CrossRef] [PubMed]
  2. Wall, R.; Rehm, J.; Fischer, B.; Brand’s, B.; Gliksman, L.; Stewart, J.; Medved, W.; Blake, J. Social costs of untreated opioid dependence. J. Urban Health 2000, 77, 688–722. [Google Scholar] [CrossRef] [PubMed]
  3. Hedegaard, H.; Miniño, A.M.; Spencer, M.R.; Warner, M. Drug overdose deaths in the United States, 1999–2020. NCHS Data Brief. 2021, 426, 1–8. [Google Scholar] [CrossRef]
  4. Florence, C.S.; Zhou, C.; Luo, F.; Xu, L. The Economic Burden of Prescription Opioid Overdose, Abuse, and Dependence in the United States, 2013. Med. Care 2016, 54, 901–906. [Google Scholar] [CrossRef]
  5. Cheung, A.; Marchand, J.; Mark, P. Loss of Life and Labor Productivity: The Canadian Opioid Crisis. Annals Am. Acad. Pol. & Soc. Sci. 2022, 703, 303–323. [Google Scholar] [CrossRef]
  6. Robertson, A.G.; Easter, M.M.; Lin, H.J.; Frisman, L.K.; Swanson, J.W.; Swartz, M.S. Associations between pharmacotherapy for opioid dependence and clinical and criminal justice outcomes among adults with co-occurring serious mental illness. J. Subst. Abuse Treat. 2018, 86, 17–25. [Google Scholar] [CrossRef]
  7. Grella, C.E.; Ostile, E.; Scott, C.K.; Dennis, M.; Carnavale, J. A Scoping Review of Barriers and Facilitators to Implementation of Medications for Treatment of Opioid Use Disorder within the Criminal Justice System. Int. J. Drug Policy 2020, 81, 102768. [Google Scholar] [CrossRef]
  8. Tsuang, M.T.; Lyons, M.J.; Meyer, J.M.; Doyle, T.; Eisen, S.A.; Goldberg, J.; True, W.; Lin, N.; Toomey, R.; Eaves, L. Co-occurrence of abuse of different drugs in men: the role of drug-specific and shared vulnerabilities. Arch. Gen. Psychiatry 1998, 55, 967–972. [Google Scholar] [CrossRef]
  9. Kendler, K.S.; Karkowski, L.M.; Neale, M.C.; Prescott, C.A. Illicit psychoactive substance use, heavy use, abuse, and dependence in a US population-based sample of male twins. Arch. Gen. Psychiatry 2000, 57, 261–269. [Google Scholar] [CrossRef]
  10. Tsuang, M.T.; Lyons, M.J.; Eisen, S.A.; Goldberg, J.; True, W.; Lin, N.; Meyer, J.M.; Toomey, R.; Faraone, S.V.; Eaves, L. Genetic influences on DSM-III-R drug abuse and dependence: a study of 3,372 twin pairs. Am. J. Med. Genet. 1996, 67, 473–477. [Google Scholar] [CrossRef]
  11. Le Foll, B.; Gallo, A.; Le Strat, Y.; Lu, L.; Gorwood, P. Genetics of dopamine receptors and drug addiction: a comprehensive review. Behav. Pharmacol. 2009, 20, 1–17. [Google Scholar] [CrossRef] [PubMed]
  12. Chen, D.; Liu, F.; Shang, Q.; Song, X.; Miao, X.; Wang, Z. Association between polymorphisms of DRD2 and DRD4 and opioid dependence: Evidence from the current studies. Am. J. Med. Genet., Part B Neuropsychiatr. Genet. 2011, 156, 661–670. [Google Scholar] [CrossRef] [PubMed]
  13. Koo, J.W.; Mazei-Robison, M.S.; Chaudhury, D.; Juarez, B.; LaPlant, Q.; Ferguson, D.; Feng, J.; Sun, H.; Scobie, K.N.; Damez-Werno, D.; Crumiller, M.; Ohnishi, Y.N.; Ohnishi, Y.H.; Mouzon, E.; Dietz, D.M.; Lobo, M.K.; Neve, R.L.; Russo, S.J.; Han, M.H.; Nestler, E.J. BDNF is a negative modulator of morphine action. Science 2012, 338, 124–8. [Google Scholar] [CrossRef] [PubMed]
  14. Haerian, B.S.; Haerian, M.S. OPRM1 rs1799971 polymorphism and opioid dependence: evidence from a meta-analysis. Pharmacogenomics 2013, 14, 813–824. [Google Scholar] [CrossRef] [PubMed]
  15. Beer, B.; Erb, R.; Pavlic, M.; Ulmer, H.; Giacomuzzi, S. Association of Polymorphisms in Pharmacogenetic Candidate Genes (OPRD1, GAL, ABCB1, OPRM1) with Opioid Dependence in European Population: A Case-Control Study. PLoS ONE 2013, 8, 75359. [Google Scholar] [CrossRef]
  16. Crist, R.; Ambrose-Lanci, L.; Vaswani, M.; Clarke, T.; Zeng, A.; Yuan, C.; Ferraro, T.; Hakonarson, H.; Kampman, K.; Dackis, C.; Pettinati, H.; O’Brien, C.; Oslin, D.; Doyle, G.; Lohoff, F.; Berrettini, W. Case–control association analysis of polymorphisms in the delta-opioid receptor, OPRD1, with cocaine and opioid addicted populations. Drug Alcohol Depend. 2013, 127, 122–128. [Google Scholar] [CrossRef]
  17. Lutz, P.E.; Kieffer, B.L. The multiple facets of opioid receptor function: implications for addiction. Curr. Opin. Neurobiol. 2013, 23, 473–479. [Google Scholar] [CrossRef]
  18. Clarke, T.K.; Weiss, A.R.; Ferarro, T.N.; Kampman, K.M.; Dackis, C.A.; Pettinati, H.M.; O’brien, C.P.; Oslin, D.W.; Lohoff, F.W.; Berrettini, W.H. The dopamine receptor D2 (DRD2) SNP rs1076560 is associated with opioid addiction. Ann. Hum. Genet. 2014, 78, 33–39. [Google Scholar] [CrossRef]
  19. Bawor, M.; Dennis, B.B.; Tan, C.; Pare, G.; Varenbut, M.; Daiter, J.; Plater, C.; Worster, A.; Marsh, D.C.; Steiner, M.; Anglin, R.; Desai, D.; Thabane, L.; Samaan, Z. Contribution of BDNF and DRD2 genetic polymorphisms to continued opioid use in patients receiving methadone treatment for opioid use disorder: An observational study. Addict. Sci. Clin. Pract. 2015, 10. [Google Scholar] [CrossRef]
  20. Zhou, H.; Rentsch, C.T.; Cheng, Z.; Kember, R.L.; Nunez, Y.Z.; Sherva, R.M.; Tate, J.P.; Dao, C.; Xu, K.; Polimanti, R.; Farrer, L.A.; Justice, A.C.; Kranzler, H.R.; Gelernter, J. Association of OPRM1 Functional Coding Variant with Opioid Use Disorder: A Genome-Wide Association Study. JAMA Psychiatry 2020, 77, 1072–1080. [Google Scholar] [CrossRef]
  21. Strang, J.; Volkow, N.D.; Degenhardt, L.; Hickman, M.; Johnson, K.; Koob, G.F.; Marshall, B.D.; Tyndall, M.; Walsh, S.L. Opioid use disorder. Nat. Rev. Dis. Primers 2020, 6. [Google Scholar] [CrossRef] [PubMed]
  22. Nielsen, D.; Ji, F.; Yuferov, V.; Ho, A.; Chen, A.; Levran, O.; Ott, J.; Kreek, M. Genotype patterns that contribute to increased risk for or protection from developing heroin addiction. Mol. Psychiatry 2008, 13, 417–428. [Google Scholar] [CrossRef] [PubMed]
  23. Nielsen, D.A.; Ji, F.; Yuferov, V.; Ho, A.; He, C.; Ott, J.; Kreek, M.J. Genome-wide association study identifies genes that may contribute to risk for developing heroin addiction. Psychiatr. Genet. 2010, 20, 207–214. [Google Scholar] [CrossRef] [PubMed]
  24. Gelernter, J.; Kranzler, H.R.; Sherva, R.; Koesterer, R.; Almasy, L.; Zhao, H.; Farrer, L.A. Genome-wide association study of opioid dependence: Multiple associations mapped to calcium and potassium pathways. Biol. Psychiatry 2014, 76, 66–74. [Google Scholar] [CrossRef] [PubMed]
  25. Cheng, Z.; Zhou, H.; Sherva, R.; Farrer, L.A.; Kranzler, H.R.; Gelernter, J. Genome-wide Association Study Identifies a Regulatory Variant of RGMA Associated With Opioid Dependence in European Americans. Biol. Psychiatry 2018, 84, 762–770. [Google Scholar] [CrossRef]
  26. Nelson, E.C.; Agrawal, A.; Heath, A.C.; Bogdan, R.; Sherva, R.; Zhang, B.; Al-Hasani, R.; Bruchas, M.R.; Chou, Y.L.; Demers, C.H.; others. Evidence of CNIH3 involvement in opioid dependence. Mol. Psychiatry 2016, 21, 608–614. [Google Scholar] [CrossRef]
  27. Dunn, K.E.; Huhn, A.S.; Finan, P.H.; Mange, A.; Bergeria, C.L.; Maher, B.S.; Rabinowitz, J.A.; Strain, E.C.; Antoine, D. Polymorphisms in the A118G SNP of the OPRM1 gene produce different experiences of opioids: A human laboratory phenotype–genotype assessment. Addict. Biol. 2024, 29, e13355. [Google Scholar] [CrossRef]
  28. Miranda, M.; Morici, J.F.; Zanoni, M.B.; Bekinschtein, P. Brain-derived neurotrophic factor: a key molecule for memory in the healthy and the pathological brain. Front. Cell. Neurosci. 2019, 13, 472800. [Google Scholar] [CrossRef]
  29. Liu, Q.R.; Walther, D.; Drgon, T.; Polesskaya, O.; Lesnick, T.G.; Strain, K.J.; De Andrade, M.; Bower, J.H.; Maraganore, D.M.; Uhl, G.R. Human brain derived neurotrophic factor (BDNF) genes, splicing patterns, and assessments of associations with substance abuse and Parkinson’s Disease. Am. J. Med. Genet. Part B Neuropsychiatr. Genet. 2005, 134, 93–103. [Google Scholar] [CrossRef]
  30. Crist, R.C.; Reiner, B.C.; Berrettini, W.H. A review of opioid addiction genetics. Curr. Opin. Psychol. 2019, 27, 31–35. [Google Scholar] [CrossRef]
  31. Gaddis, N.; Mathur, R.; Marks, J.; Zhou, L.; Quach, B.; Waldrop, A.; Levran, O.; Agrawal, A.; Randesi, M.; Adelson, M.; Jeffries, P.W.; Martin, N.G.; Degenhardt, L.; Montgomery, G.W.; Wetherill, L.; Lai, D.; Bucholz, K.; Foroud, T.; Porjesz, B.; Runarsdottir, V.; Tyrfingsson, T.; Einarsson, G.; Gudbjartsson, D.F.; Webb, B.T.; Crist, R.C.; Kranzler, H.R.; Sherva, R.; Zhou, H.; Hulse, G.; Wildenauer, D.; Kelty, E.; Attia, J.; Holliday, E.G.; McEvoy, M.; Scott, R.J.; Schwab, S.G.; Maher, B.S.; Gruza, R.; Kreek, M.J.; Nelson, E.C.; Thorgeirsson, T.; Stefansson, K.; Berrettini, W.H.; Gelernter, J.; Edenberg, H.J.; Bierut, L.; Hancock, D.B.; Johnson, E.O. Multi-trait genome-wide association study of opioid addiction: OPRM1 and beyond. Sci. Rep. 2022, 12. [Google Scholar] [CrossRef] [PubMed]
  32. Gondré-Lewis, M.C.; Elman, I.; Alim, T.; Chapman, E.; Settles-Reaves, B.; Galvao, C.; Gold, M.S.; Baron, D.; Kazmi, S.; Gardner, E.; Gupta, A.; Dennen, C.; Blum, K. Frequency of the Dopamine Receptor D3 (rs6280) vs. Opioid Receptor µ1 (rs1799971) Polymorphic Risk Alleles in Patients with Opioid Use Disorder: A Preponderance of Dopaminergic Mechanisms? Biomedicines 2022, 10, 870. [Google Scholar] [CrossRef]
  33. Freiermuth, C.E.; Kisor, D.F.; Lambert, J.; Braun, R.; Frey, J.A.; Bachmann, D.J.; Bischof, J.J.; Lyons, M.S.; Pantalon, M.V.; Punches, B.E.; others. Genetic variants associated with opioid use disorder. Clin. Pharmacol. Ther. 2023, 113, 1089–1095. [Google Scholar] [CrossRef]
  34. Ho, A.M.; Tang, N.L.; Cheung, B.K.; Stadlin, A. Dopamine receptor D4 gene -521C/T polymorphism is associated with opioid dependence through cold-pain responses. Annals of the New York Academy of Sciences. Blackwell Publishing Inc., 2008, Vol. 1139, pp. 20–26. [CrossRef]
  35. Chung, P.; Logge, W.B.; Riordan, B.C.; Haber, P.S.; Merriman, M.E.; Phipps-Green, A.; Topless, R.K.; Merriman, T.R.; Conner, T.; Morley, K.C. Genetic polymorphisms on OPRM1, DRD2, DRD4, and COMT in young adults: lack of association with alcohol consumption. Front. Psychiatry 2020, 11, 549429. [Google Scholar] [CrossRef] [PubMed]
  36. Pedrosa, E.; Kaushik, S.; Lachman, H.M. ChIP-chip analysis of neurexins and other candidate genes for addiction and neuropsychiatric disorders. J. Neurogenet. 2010, 24, 5–17. [Google Scholar] [CrossRef]
  37. Stoltenberg, S.F.; Lehmann, M.K.; Christ, C.C.; Hersrud, S.L.; Davies, G.E. Associations among types of impulsivity, substance use problems and neurexin-3 polymorphisms. Drug Alcohol Depend. 2011, 119, e31–e38. [Google Scholar] [CrossRef] [PubMed]
  38. McLaren, W.; Gil, L.; Hunt, S.E.; Riat, H.S.; Ritchie, G.R.; Thormann, A.; Flicek, P.; Cunningham, F. The ensembl variant effect predictor. Genome Biol. 2016, 17, 1–14. [Google Scholar] [CrossRef]
  39. Purcell, S.; Neale, B.; Todd-Brown, K.; Thomas, L.; Ferreira, M.A.; Bender, D.; Maller, J.; Sklar, P.; De Bakker, P.I.; Daly, M.J.; others. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 2007, 81, 559–575. [Google Scholar] [CrossRef]
  40. Bogari, N.M.; Al-Allaf, F.A.; Aljohani, A.; Taher, M.M.; Qutub, N.A.; Alhelfawi, S.; Alobaidi, A.; Alqudah, D.M.; Banni, H.; Dairi, G.; others. The co-existence of ADHD with autism in Saudi children: an analysis using next-generation DNA sequencing. Front. Genet. 2020, 11, 548559. [Google Scholar] [CrossRef]
  41. Lie, E.; Ko, J.S.; Choi, S.Y.; Roh, J.D.; Cho, Y.S.; Noh, R.; Kim, D.; Li, Y.; Kang, H.; Choi, T.Y.; others. SALM4 suppresses excitatory synapse development by cis-inhibiting trans-synaptic SALM3–LAR adhesion. Nat. Commun. 2016, 7, 12328. [Google Scholar] [CrossRef]
  42. Lie, E.; Yeo, Y.; Lee, E.J.; Shin, W.; Kim, K.; Han, K.A.; Yang, E.; Choi, T.Y.; Bae, M.; Lee, S.; others. SALM4 negatively regulates NMDA receptor function and fear memory consolidation. Commun. Biol. 2021, 4, 1138. [Google Scholar] [CrossRef] [PubMed]
  43. Fill, M.; Copello, J.A. Ryanodine receptor calcium release channels. Physiol. Rev. 2002, 82, 893–922. [Google Scholar] [CrossRef] [PubMed]
  44. Santulli, G.; Nakashima, R.; Yuan, Q.; Marks, A.R. Intracellular calcium release channels: an update. J. Physiol. 2017, 595, 3041–3051. [Google Scholar] [CrossRef] [PubMed]
  45. Brini, M.; Calì, T.; Ottolini, D.; Carafoli, E. Neuronal calcium signaling: function and dysfunction. Cell. Mol. Life Sci. 2014, 71, 2787–2814. [Google Scholar] [CrossRef] [PubMed]
  46. Li, X.; Zhu, C.; Tu, W.H.; Yang, N.; Qin, H.; Sun, Z. ZMIZ1 preferably enhances the transcriptional activity of androgen receptor with short polyglutamine tract. PLoS One 2011, 6, e25040. [Google Scholar] [CrossRef] [PubMed]
  47. Tobiansky, D.J.; Wallin-Miller, K.G.; Floresco, S.B.; Wood, R.I.; Soma, K.K. Androgen regulation of the mesocorticolimbic system and executive function. Front. Endocrinol. 2018, 9, 279. [Google Scholar] [CrossRef]
  48. Szklarczyk, D.; Kirsch, R.; Koutrouli, M.; Nastou, K.; Mehryary, F.; Hachilif, R.; Gable, A.L.; Fang, T.; Doncheva, N.T.; Pyysalo, S.; others. The STRING database in 2023: protein–protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res. 2023, 51, D638–D646. [Google Scholar] [CrossRef]
  49. Raudvere, U.; Kolberg, L.; Kuzmin, I.; Arak, T.; Adler, P.; Peterson, H.; Vilo, J. g: Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 2019, 47, W191–W198. [Google Scholar] [CrossRef]
  50. Van der Auwera, G.A.; O’Connor, B.D. Genomics in the cloud: using Docker, GATK, and WDL in Terra; O’Reilly Media, 2020.
  51. Poplin, R.; Ruano-Rubio, V.; DePristo, M.A.; Fennell, T.J.; Carneiro, M.O.; Van der Auwera, G.A.; Kling, D.E.; Gauthier, L.D.; Levy-Moonshine, A.; Roazen, D. ; others. Scaling accurate genetic variant discovery to tens of thousands of samples. BioRxiv, 2011; 201178. [Google Scholar] [CrossRef]
  52. Caetano-Anolles, D. (How to) Filter variants either with VQSR or by hard-filtering. GATK [Internet], 2023. [Google Scholar]
  53. Danecek, P.; Bonfield, J.K.; Liddle, J.; Marshall, J.; Ohan, V.; Pollard, M.O.; Whitwham, A.; Keane, T.; McCarthy, S.A.; Davies, R.M.; Li, H. Twelve years of SAMtools and BCFtools. GigaScience 2021, 10. [Google Scholar] [CrossRef]
  54. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504. [Google Scholar] [CrossRef]
Preprints 120336 g001
Figure 1. Protein-Protein interaction (PPI) network of high confidence (normalized interaction score > 0.8, details in Section 4.4) PPIs subset from STRING database. The genes harboring a single variant with a very high alternate allele count exclusively (Table 3) in cases are shown with orange oval nodes, while genes harboring multiple variants with a high alternate allele count exclusively in cases (Table 4) are shown in green oval nodes/genes. All other gray nodes are either immediate neighbors (directly interacted by PPIs) of green or orange nodes or nodes having at least two neighbors in the list of selected nodes/genes. Finally, the subnetwork is created using selected nodes and all the edges.
Figure 1. Protein-Protein interaction (PPI) network of high confidence (normalized interaction score > 0.8, details in Section 4.4) PPIs subset from STRING database. The genes harboring a single variant with a very high alternate allele count exclusively (Table 3) in cases are shown with orange oval nodes, while genes harboring multiple variants with a high alternate allele count exclusively in cases (Table 4) are shown in green oval nodes/genes. All other gray nodes are either immediate neighbors (directly interacted by PPIs) of green or orange nodes or nodes having at least two neighbors in the list of selected nodes/genes. Finally, the subnetwork is created using selected nodes and all the edges.
Preprints 120336 g001
Table 1. List of genes reported in literature harboring variants associated with OAD risk.
Table 1. List of genes reported in literature harboring variants associated with OAD risk.
Gene Chra Gene BP(beg)b Gene BP(end)c Variant (rsID) Sample population Opioid exposed controls Reference
OPRM1 chr6 154039240 154132356 rs1799971 EuAd, AfAe Yes [14,15,20]
OPRD1 chr1 28812170 28871267 rs2236861 EuA, AfA; EuA, AfA Yes; No [15,16]
DRD2 chr11 113409605 113475398 rs1799978 Euf, NaAmg, Asnh, Peri Yes [19]
BDNF chr11 27654893 27700455 rs6265 Euf, NaAmg, Asnh, Peri Yes [19]
APBB2 chr4 40810027 41214542 rs114070671 EuAmk, AfAmm No [24]
KCNG2 chr18 79797938 79900100 rs62103177 EuAm, AfAm No [24,31]
KCNC1 chr11 17734781 17783057 rs60349741 EuAm, AfAm No [24]
CNIH3 chr1 224616317 224740554 rs1436175, rs10799590, rs12130499, rs298733, rs1436171, and rs1369846 EuAun Yes [26]
RGMA chr15 93035271 93089211 rs12442183 EuAm Yes [25]
DRD3 chr3 114127580 114179052 rs324029 and rs2654754 EuAm, AfAm No [32,33]
DRD4 chr11 637269 640706 rs1800955 Chinese males No [34,35]
NRXN3 chr14 78170373 79868291 rs8019381 p Caucasians No [36,37]
a Chromosome; b Gene’s beginning base position; c Gene’s ending base position; d European ancestry; e African ancestry; f European; g Native South/North American; h Asian; i Persian; k European-American; m African-American; n European-Australian; p Associated to substance use disorder
Table 2. The non-coding consequence variants reaching significance (raw p-value < 0.01) with association to OAD in the extended genes loci that are listed in Table 1.
Table 2. The non-coding consequence variants reaching significance (raw p-value < 0.01) with association to OAD in the extended genes loci that are listed in Table 1.
Chra Variantb or rsID BPc Gene F_Ad F_Ue Conseqn.f χ 2 p-value ORg
3 TGAAA>T 114142739 DRD3 0.0909 0.5455 ivh 10.48 0.0012 0.0833
18 rs76838079 79873271 KCNG2 0.3182 0 iv 8.324 0.0039 NA
3 rs73232565 114124222 - 0 0.2727 - 7.527 0.0061 0
11 rs3051820 17785864 - 0.2727 0.6818 - 7.379 0.0066 0.175
4 rs1011069 41217734 - 0.2727 0 - 6.947 0.0084 NA
14 rs143010574 79227804 NRXN3 0.2727 0 gutvi, iv 6.947 0.0084 NA
4 rs7695309 41216892 - 0.3636 0.0455 - 6.844 0.0089 12
14 rs7145683 79241818 NRXN3 0.3636 0.0455 gutv, iv 6.844 0.0089 12
14 G>GA 79242068 NRXN3 0.3636 0.0455 gutv, iv 6.844 0.0089 12
14 rs12889183 79360638 NRXN3 0.5 0.1364 iv 6.705 0.0096 6.333
14 rs11625994 79364803 NRXN3 0.5 0.1364 iv 6.705 0.0096 6.333
14 rs8008332 79365491 NRXN3 0.5 0.1364 iv 6.705 0.0096 6.333
14 rs2167150 79367835 NRXN3 0.5 0.1364 iv 6.705 0.0096 6.333
a Chromosome; b rsID if found is listed otherwise the variant is listed in this column in form X>Y, where X is reference allele while Y represents alternate allele; c Base position; d minor allele frequency in cases; e minor allele frequency in controls; f variants’ consequence; g odds ratio; h intronic variant; i genic upstream transcript variant
Table 3. The variants across the genome are ranked by consequence and filtered as per significance associated with OAD.
Table 3. The variants across the genome are ranked by consequence and filtered as per significance associated with OAD.
Chra Variantb or rsID BPc Gene C_Ad Conseqn.e χ 2 p-value
12 rs773026868 50352078 FAM186A 11 fsvf 13.25 0.000272
11 rs60494098 9091455 SCUBE2 10 mvg 12.94 0.000321
1 rs10907376 223394461 CCDC185 9 mv 11.31 0.000769
13 C>A 29324631 MTUS2 8 mv 9.778 0.001766
7 rs3750050 77627396 PTPN12 7 mv 8.324 0.003912
7 rs1046515 140694787 ADCK2 7 mv 8.324 0.003912
16 rs308925 77735937 NUDT7 7 mv 8.324 0.003912
2 A>G 207613128 METTL21A 7 mv 8.324 0.003912
6 rs1048886 70579486 SDHAF4 7 mv 8.324 0.003912
16 rs3869427 69954416 CLEC18A 7 mv 8.324 0.003912
1 rs147489167 248574363 OR2T34 7 mv 8.324 0.003912
a Chromosome; b rsID if found is listed otherwise the variant is listed in this column in form X>Y, where X is reference allele while Y represents alternate allele; c Base position; d alternate allele count in cases; e varinats’ consequence; f frameshift variant; g missense variant
Table 4. Details of co-occurring variants on genes, where the sum of alternate allele count in cases are greater or equal to 10 but zero in controls.
Table 4. Details of co-occurring variants on genes, where the sum of alternate allele count in cases are greater or equal to 10 but zero in controls.
Chra Gene #variantsb list of C_Ac list of Consequencesd
10 ZMIZ1 7 5, 5, 5, 4, 4, 4, 4 fsve, mvf, fsv, fsv, fsv, fsv, fsv
19 LRFN3 6 5, 5, 5,5, 5, 5 idg, pavh, iik, fsv, fsv, fsv
9 OR1L6 4 4, 4, 4, 4 mv, mv, mv, mv
15 RYR3 3 5, 5, 4 fsv, fsv, sgm & fsv
10 PWWP2B 3 4, 4, 4 sg, mv, fsv
7 ZNF92 2 6, 6 mv, mv
19 CYP4F12 3 4, 4, 4 mv, mv, sdvn & ntvp
10 NUTM2D 2 5, 5 mv, mv
a Chromosome; b number of variants; c alternate allele counts in cases for each variant; d variants’ consequence; e frameshift variant; f missense variant; g inframe deletion; h protein altering variant; k inframe insertion; m stop gained; n splice donor variant; p nmd transcript variant
Table 5. Description of annotated terms shown in Figure 2. Terms relevant to addiction are shown with underlined-text.
Table 5. Description of annotated terms shown in Figure 2. Terms relevant to addiction are shown with underlined-text.
ID Source Term ID Term name P a d j (query)
1 KEGG KEGG:04020 Calcium signaling pathway 1.780 × 10 − 15
2 KEGG KEGG:04713 Circadian entrainment 8.371 × 10 − 14
3 GO:BP GO:0014808 release of sequestered calcium ion into cytosol by sarcoplasmic reticulum 1.334 × 10 − 12
4 KEGG KEGG:04728 Dopaminergic synapse 8.377 × 10 − 12
5 GO:BP GO:0051208 sequestering of calcium ion 2.017 × 10 − 11
6 KEGG KEGG:04340 Hedgehog signaling pathway 4.704 × 10 − 11
7 GO:CC GO:0033017 sarcoplasmic reticulum membrane 3.405 × 10 − 10
8 WP WP:WP3929 Chemokine signaling pathway 9.043 × 10 − 10
9 KEGG KEGG:04921 Oxytocin signaling pathway 1.489 × 10 − 9
10 KEGG KEGG:04915 Estrogen signaling pathway 4.481 × 10 − 9
11 REAC REAC:R-HSA-5578775 Ion homeostasis 5.472 × 10 − 9
12 GO:BP GO:0010646 regulation of cell communication 9.459 × 10 − 9
13 KEGG KEGG:05032 Morphine addiction 1.816 × 10 − 7
14 KEGG KEGG:05034 Alcoholism 2.403 × 10 − 7
15 GO:MF GO:0005219 ryanodine-sensitive calcium-release channel activity 4.909 × 10 − 7
16 REAC REAC:R-HSA-9006934 Signaling by Receptor Tyrosine Kinases 1.062 × 10 − 6
17 GO:MF GO:0140096 catalytic activity, acting on a protein 9.842 × 10 − 7
18 GO:BP GO:0006942 regulation of striated muscle contraction 2.257 × 10 − 6
19 GO:BP GO:0036211 protein modification process 9.419 × 10 − 6
20 GO:BP GO:1901564 organonitrogen compound metabolic process 9.464 × 10 − 6
21 GO:CC GO:0098797 plasma membrane protein complex 8.009 × 10 − 6
22 KEGG KEGG:05031 Amphetamine addiction 6.351 × 10 − 6
23 REAC REAC:R-HSA-180292 GAB1 signalosome 5.328 × 10 − 5
24 REAC REAC:R-HSA-111885 Opioid Signalling 2.349 × 10 − 4
25 GO:MF GO:0005509 calcium ion binding 2.788 × 10 − 4
26 GO:MF GO:0005102 signaling receptor binding 2.920 × 10 − 4
27 GO:BP GO:0009725 response to hormone 5.256 × 10 − 4
28 WP WP:WP3680 Physico chemical features and toxicity associated pathways 2.095 × 10 − 4
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Alerts
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2025 MDPI (Basel, Switzerland) unless otherwise stated