Front. PsychiatryFrontiers in PsychiatryFront. Psychiatry1664-0640Frontiers Media S.A.10.3389/fpsyt.2022.1078894PsychiatryOriginal ResearchCpH methylome analysis in human cortical neurons identifies novel gene pathways and drug targets for opioid use disorderNagamatsuSheila T.123RompalaGregory4HurdYasmin L.4Núñez-RiosDiana L.123Montalvo-OrtizJanitza L.123*Traumatic Stress Brain Research GroupAlvarezVictor E.BenedekDavidCheAliciaCruzDianne A.DavisDavid A.GirgentiMatthew J.HoffmanEllenHoltzheimerPaul E.HuberBertrand R.KayeAlfredKrystalJohn H.LabadorfAdam T.KeaneTerence M.LogueMark W.McKeeAnnMarxBrianMashDeborahMillerMark W.NollerCrystalJM-OScottWilliam K.SchnurrPaulaSteinThorUrsanoRobertWilliamsonDouglas E.WolfErika J.YoungKeith A.1Division of Human Genetics, Department of Psychiatry, Yale University School of Medicine, New Haven, CT, United States2VA Connecticut (VA CT) Healthcare Center, West Haven, CT, United States3Clinical Neurosciences Division, U.S. Department of Veterans Affairs National Center of Posttraumatic Stress Disorder, West Haven, CT, United States4Icahn School of Medicine at Mount Sinai, New York, NY, United States
Edited by: Chao Chen, Central South University, China
Reviewed by: Jean Lud Cadet, National Institute on Drug Abuse (NIH), United States; Caesar Li, Johnson & Johnson, United States
*Correspondence: Janitza L. Montalvo-Ortiz, janitza.montalvo-ortiz@yale.edu
This article was submitted to Behavioral and Psychiatric Genetics, a section of the journal Frontiers in Psychiatry
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Introduction
DNA methylation (DNAm), an epigenetic mechanism, has been associated with opioid use disorder (OUD) in preclinical and human studies. However, most of the studies have focused on DNAm at CpG sites. DNAm at non-CpG sites (mCpHs, where H indicates A, T, or C) has been recently shown to have a role in gene regulation and to be highly abundant in neurons. However, its role in OUD is unknown. This work aims to evaluate mCpHs in the human postmortem orbital frontal cortex (OFC) in the context of OUD.
Methods
A total of 38 Postmortem OFC samples were obtained from the VA Brain Bank (OUD = 12; Control = 26). mCpHs were assessed using reduced representation oxidative bisulfite sequencing in neuronal nuclei. Differential analysis was performed using the “methylkit” R package. Age, ancestry, postmortem interval, PTSD, and smoking status were included as covariates. Significant mCpHs were set at q-value < 0.05. Gene Ontology (GO) and KEGG enrichment analyses were performed for the annotated genes of all differential mCpH loci using String, ShinyGO, and amiGO software. Further, all annotated genes were analyzed using the Drug gene interaction database (DGIdb).
Results
A total of 2,352 differentially methylated genome-wide significant mCpHs were identified in OUD, mapping to 2,081 genes. GO analysis of genes with differential mCpH loci showed enrichment for nervous system development (p-value = 2.32E-19). KEGG enrichment analysis identified axon guidance and glutamatergic synapse (FDR 9E-4–2.1E-2). Drug interaction analysis found 3,420 interactions between the annotated genes and drugs, identifying interactions with 15 opioid-related drugs, including lofexidine and tizanidine, both previously used for the treatment of OUD-related symptoms.
Conclusion
Our findings suggest a role of mCpHs for OUD in cortical neurons and reveal important biological pathways and drug targets associated with the disorder.
opioidmethylationnon-CpG sitepostmortem human brainorbitofrontal cortexepigenetic1IK2CX002095-01A1R21DA050160U.S. Department of Veterans Affairs10.13039/100000738National Institute on Drug Abuse10.13039/100000026Introduction
Opioid use disorder (OUD) is a chronic lifelong disorder at epidemic levels in the United States (1). OUD has been associated with a high disease burden and overdose death rates (17.8 per 100,000 individuals) (2). While genetic risk factors have been identified in recent large-scale genome-wide association studies (GWAS) of OUD (3, 4), these explain only part of the variance. Environmental factors interplay with the genetic background to influence OUD risk, possibly via epigenetic mechanisms. Alterations in DNA methylation (DNAm), one of the epigenetic mechanisms, have been associated with OUD. DNAm generally occurs in cytosine-guanine dinucleotides linked by a phosphate group (CpG site) and is associated with gene repression when located in the CpG island of the promoter region (5). Our group and others have identified differentially methylated CpG sites associated with OUD in studies examining peripheral tissue (6–8).
In the human postmortem brain, a previous study evaluating DNAm using the Infinium HumanMethylation450K BeadChip identified 1,298 differential 5mC in the human postmortem orbitofrontal cortex (OFC) of heroin users (p < 0.001) (9). The OFC is known to be involved in decision-making processes related to drug addiction (10–12). A recent study evaluating DNAm using the Illumina MethylationEPIC BeadChip (EPIC) array in the dorsolateral prefrontal cortex (dlPFC) reported no significant 5mC CpGs associated with opioid intoxication (13). Another DNAm study in postmortem brain tissue using the same array conducted a co-methylation analysis and identified 6 co-methylated modules related to OUD, involved in response to organic substance and astrocyte and glial differentiation (14). Further, we recently examined neuronal-specific genome-wide 5mC profiles associated with OUD in the OFC, in parallel with 5-hydroxymethylcytosine (5hmC), and identified 397 and 1,740 5mC and 5hmC differential CpGs, respectively (15). We observed epigenetically dysregulated genes implicated in pain signaling and observed enrichment for the Wnt signaling pathway and the G-protein signaling pathway (15). Even though there has been evident progress in the study of DNAm in OUD (or related traits) in recent years, these studies have been limited to examining DNAm at CpG sites using array-based technology.
Compared to array-based technology, single-base resolution sequencing for epigenetic modification detection has increased accuracy, coverage, and resolution in DNAm assessment (16). Bisulfite sequencing has been widely applied in a targeted approach to assessing DNAm of opioid receptors in humans [e.g., child abuse (17), alcohol dependence (18)], and animal models [e.g., pain (19)]. Studies utilizing genome-wide approaches assessing mCpG have been performed in different tissues, including human postmortem brain in the context of schizophrenia (20) and autism (21). However, no previous work has focused on substance use disorders, including OUD.
DNA methylation can also occur in cytosines phosphate-linked to adenine, thymine, and cytosine, known as non-CpG sites (mCpHs). mCpHs is usually low in adult peripheral somatic cells, corresponding to only 0.02% (5). However, a high prevalence of mCpHs is observed in embryonic cells and brain tissue (5). In the brain, mCpHs are highly enriched in neuronal cells, and studies have suggested having a role in gene regulation and brain function (5). Despite mCpHs representing a promising target for the study of brain disorders, little work has focused on evaluating the role of mCpHs in disease etiology. Further, since most DNAm studies have been performed using array-based technology, which only assesses DNAm at CpG sites, very limited work has evaluated mCpH.
In this study, we examined mCpHs at the genome-wide level in the OFC neuronal nuclei of OUD subjects compared with non-OUD subjects. We identified 2,351 differential mCpHs associated with OUD, of which 1,513 were hypomethylated and 838 hypermethylated. Functional annotation analyses showed enrichment for cell-cell communication, glutamatergic synapses, and cholinergic synapses. Gene drug interaction analysis detected 5,690 interactions, including opioid-related drugs as potential drug targets for OUD treatment. Our study highlights the role of mCpH in OUD and reveals important biological pathways and drug targets associated with the disorder.
Materials and methodsSample description
Postmortem human brain OFC from 38 individuals were collected from the brain tissue repository at the National Post-Traumatic Stress Disorder (PTSD) Brain Bank (NPBB) (22), a biorepository from the U.S. Department of Veterans Affairs (VA). Demographics and clinical characteristics are depicted in Table 1.
Demographics of the postmortem OFC specimens obtained from the non-OUD and OUD groups.
OUD− (N = 26)
OUD + (N = 12)
Ancestry
African American European American
5
4
19
8
PMI (μ ± SD)
30.65 ± 8.15
29.6 ± 7.5
Age (μ ± SD)
43.1 ± 11.55
37.6 ± 8.9
Cigarette smoking
13
10
Alcohol dependence
5
3
PTSD
13
12
Fluorescence-activated nuclei sorting, DNA extraction, and reduced representation oxidative bisulfite sequencing
A total of 100–200 mg of OFC tissue was lysed in homogenization buffer composed of 0.1% Triton, 0.32 M sucrose, 5 mM CaCl2, 3 mM MgCl2, and 10 mM Tris-HCl. After, they were filtered, loaded onto a sucrose cushion, and ultracentrifuged. The nuclei were resuspended in bovine serum albumin with Anti-NeuN-PE being added 4′,6-diamidino-2-phenylindole (DAPI) before sorting. FANS procedure was conducted on a BD 5-laser cell sorting system at the Icahn School of Medicine Flow Cytometry CoRE. The NeuN + nuclei collected (0.5–1 M) were pelleted by centrifugation and processed for DNA extraction according to the DNeasy Blood and Tissue Kit (Cat. #69504, Qiagen) manufacturer’s protocol. Eluted samples were concentrated and stored at −80°C. A total of 400 ng of DNA was used to prepare the methylation library for bisulfite treatment using the NuGEN Ovation RRoxBS Methyl-Seq library preparation kit. The bisulfite treatment converts unmethylated cytosines into uracils. Sequencing was conducted using the Illumina NovaSeq6000 system. Library preparation and sequencing were carried out at the Weill Cornell Epigenomics Core (New York, NY). The complete protocol is described in Rompala et al. (15).
mCpHs differential analysis
Sequencing reads were mapped to GRCh38, and converted and unconverted cytosines were detected using the Bismark bisulfite read mapper (23). The % of unconverted cytosines was used to represent the beta values. mCpHs were filtered for a minimum of 10x coverage in all samples and further normalized for coverage variability correction. Differential analysis was conducted for single mCpHs using the methylkit R package (24). Age of death, cigarette smoking, ancestry, postmortem interval (PMI), and posttraumatic stress disorder (PTSD) were included as covariates. The false discovery rate (FDR) significance threshold was set as <0.05.
mCpH annotation
mCpH annotation was performed using genomation (25) and biomaRt (26) R packages. First, annotation of the closest Ensembl transcript ID to the mCpH site was performed using genomation. Briefly, after coercing mCpH sites to GRanges, we used the “annotateWithGeneParts” function to annotate mCpH located at promoters/introns/exons/intergenic regions, followed by the “annotateWithFeatureFlank” function, which performs the annotation of CpG island/shores/flanking regions. Further, the “getAssociationWithTSS” function in genomation was used to calculate the distance to the nearest TSS and annotate the mCpH region with the Ensembl transcript ID. During this step, no threshold was applied for the distance between the mCpH and the annotated gene (TSS distance is provided in Supplementary Table 1). The annotation of gene symbols, descriptions, and biotypes was performed using BiomaRt. By using a known gene dataset from the UCSC genome browser (27), we confirmed the distance between the BiomaRt annotated symbol and the mCpH, where we applied a threshold of 1,500 bp distance. The location of the mCpH was assigned by including 0 = inside the gene, negative number = upstream, and positive number = downstream (Supplementary Table 1). Gene symbols were used for the functional annotation, GWAS enrichment, and drug interaction analysis. For the input for functional annotation, we used two different approaches: All genes and a subset excluding mCpHs in the intergenic regions (i.e., not annotated in the column “Distance <1,500 Up and downstream,” Supplementary Table 1). For GWAS enrichment and drug interaction analysis, we evaluated all annotated genes.
Functional annotation analysis
Functional annotation was conducted using the enrichR package (28) via the web-based tool “Enrichr.” Enrichr performs enrichment analysis with gene sets built based on prior biological knowledge. We selected eight databases, including (1) gene ontology (GO) database, which contains information about gene function describing the process and cellular location (GO_Molecular_Function_2021, GO_Cellular_Component_2021, and GO_Biological_Process_2021), (2) COVID-19 database that includes a collection of genes related to COVID-19 (COVID-19_Related_Gene_Sets_2021), (3) KEGG database that has a collection of known pathways in each of the genes involved (KEGG_2021_Human), (4) Wikipathways, a collaborative platform that maintains a curated collection of biological pathways (WikiPathway_2021_Human), (5) Allen Brain Atlas, an atlas of gene expression integrating it with function and structure (Allen_Brain_Atlas_10x_scRNA_2021), (6) MAGMA database, which contains a collection of gene sets associated with drugs and diseases (MAGMA_Drugs_and_Diseases). Further, we conducted an enrichment analysis using STRING (29), a database of known protein-protein interaction (PPI), which also evaluates databases such as (1) InterPro, which uses protein families and predicted domains to provide the functional annotation, (2) Reactome pathways, a manually curated pathway database, and (3) Tissue expression. We conducted the functional annotation analysis for the whole set of annotated genes (Supplementary Table 2), and for a subset removing sites located in the intergenic regions with a distance greater than 1,500 bp up and downstream using the UCSC software (Supplementary Table 3). Functional annotation for the overlapped annotated genes between mCpHs and 5mC in CpG sites (Supplementary Table 4) was performed using STRING. PPI analysis was performed using the STRING software (29) with the highest confidence score (0.9) and allowing experimentally determined interactions and co-expression. Significant sites were defined as FDR < 0.05.
Cell-type enrichment analysis
Brain cell-type enrichment analysis was carried out using transcriptomic data from human postmortem brain samples to characterize the cell-type enrichment of OUD-associated differential mCpHs in neurons, microglia, astrocytes, oligodendrocytes, and endothelial cell type-specific markers (30).
Genome-wide association enrichment analysis
Genome-wide association studies enrichment was conducted using FUMA (31). We used the Entrez IDs as input to the GENE2Function web tool. The annotation was performed using the David software (32). The default adjusted p-value from FUMA was considered to define significance (adjusted P ≤ 0.05). The Bonferroni correction threshold for the P-value was set as ≤ 1.89E-4 and described in Supplementary Table 4. We further evaluated the mCpHs annotated with the genes enriched for opioid use-related traits by examining whether the significant mCpH located upstream or downstream of the nearest genes had a known regulatory role. For this, we used JASPAR Transcription Factor Binding Site Database (33) and Encode Database, both linked to the UCSC genome browser web tool (27), that predicts transcription factor binding sites and H3K27Ac marks in the assessed region, respectively.
Drug interaction analysis
The drug interaction analysis was performed using the Drug Gene Interaction Database (DGIdb) (34) by assessing the interaction between all differential annotated genes and drugs.1 Drug interaction was also evaluated for genes involved in synaptic pathways selected from functional annotation using STRING. Cytoscape (35) was used for the visualization of identified drug interactions.
Results
Neuronal nuclei from postmortem OFC tissue of 38 donors (nOUD = 12; ncontrol = 26) were analyzed. The samples included males of European and African ancestry (Table 1).
Differential mCpHs
Individual mCpHs were evaluated for OUD differential analysis. We observed 2,351 differential mCpHs associated with OUD (Supplementary Table 1) that included 1,513 hypomethylated mCpHs and 838 hypermethylated mCpHs. Assessing the gene features of all differentially methylated CpH loci, most were located in promoter regions (Figure 1A). We compared our differential mCpHs with findings from our previous study focused on 5mC at CpG sites (15) and found that 61 annotated genes overlapped with OUD-associated genes with differential 5mC (Figure 1B and Supplementary Table 4). No significant enrichment was identified for the overlapped genes. Further, cell-type enrichment analysis was conducted to evaluate mCpH sites on gene markers for unique cell types. We observed enrichment of differential mCpH loci on neuron-specific markers, followed by endothelial, astrocytes, oligodendrocytes, and microglia-specific markers (Figure 1C). Differential mCpHs showing hypermethylation were primarily located in the intergenic region, while the hypomethylated mCpHs were more prevalent in promoter regions (Figure 2).
Evaluation of differential mCpHs associated with OUD. (A) Genomic location of differential mCpHs. (B) Cell-type enrichment of differential mCpHs sites. (C) Venn diagram showing the comparison between differential mCpHs and mCpGs.
Differential mCpHs sites associated with OUD. (A) Circus plot showing hypermethylated and hypomethylated sites. (B) Volcano plot of the differential sites. (C) Genomic location of hypermethylated and (D) hypomethylated differential sites.
From all the OUD-associated differential mCpHs identified, we were able to detect three annotated genes previously identified in the OUD GWAS literature (36), APBB2 (q-value = 0.007), CNIH3 (q-value = 4.92E-12) and candidate genes (37), the OPRK1 (q-value = 0.020). One of the top associations mapped to the dopamine neurotrophic factor, CDNF (q-value = 2.96E-25). We also identified genes previously found in epigenome-wide association studies (EWAS) of opioid dependence (6), including RERE (q-value = 4.79E-53), and CFAP77 (q-value = 3.91E-25).
Functional annotation
Functional annotation was performed using enrichR (28) and STRING (29) for all annotated genes and for a subset excluding intergenic sites. When we evaluated all annotated genes using enrichR, we observed an enrichment for glutamate receptor signaling pathway (p-value = 5.67E-05), regulation of neuron differentiation (p-value = 5.27E-04), and dopaminergic neuron differentiation (p-value = 1.24E-03) (Figure 3A). We also detected enrichment of differential mCpHs in Relationship between inflammation, COX-2 and EGFR (p-value = 1.10E-03; Figure 3B), Wnt Signaling pathway (p-value = 5.20E-04; Figure 3B), and dopaminergic neurogenesis (p-value = 3.97E-03; Figure 3B). Using STRING, we observed enrichment for the cholinergic synapse (FDR = 0.047; Figure 3C), glutamatergic synapse (FDR = 0.021; Figure 3C), cell-cell communication (FDR = 0.034; Figure 3D); and transmission across chemical synapses (FDR = 0.036; Figure 3D). Additional functional annotation results are included in Supplementary Table 2.
Functional enrichment. Enrichment analysis was conducted using EnrichR (A,B) and STRING software (C,D) for GO biological process, wikiPathways, KEGG, and RCTM.
Although we observed changes in the enrichment analysis when evaluating the subset excluding the intergenic sites, we also observed some overlapping enrichment, including glutamate receptor signaling pathway (p-value = 4.30E-05), regulation of neuron differentiation (p-value = 0.0016), dopaminergic neuron differentiation (p-value = 0.0024), Relationship between inflammation, COX-2 and EGFR (p-value = 0.0079), Wnt Signaling pathway (p-value = 0.00028), dopaminergic neurogenesis (p-value = 0.0049), cholinergic synapse (FDR = 0.022), glutamatergic synapse (FDR = 0.0084) (Supplementary Table 3). Comparing all annotated genes and the subset, we observed an overlap of enrichment terms of 66 (Allen Brain Atlas), 67 (COVID-19 related genes), 214 (GO biological terms), 20 (GO cellular component), 54 (Go molecular function), 44 (KEGG), 29 (MAGMA drugs and diseases), and 54 (WikiPathways) using enrichR, and 9 (Tissues), 1 (RCTM), 1 (PFAM), 19 (KEGG), 3 (Interpro), 3 (Diseases) using STRING.
Protein-protein interaction analysis
The PPI analysis of all annotated genes from the OUD-associated differential mCpHs was performed using STRING (29) (Figure 4). We chose nine significantly enriched pathways, including nervous system development (FDR = 4.76e-20), multicellular organism development (FDR = 1.04e-18), neurogenesis (FDR = 1.04e-18), system development (FDR = 5.26e-18), alternative splicing (FDR = 2.18e-14), regulation of multicellular organismal process (FDR = 8.62e-12), regulation of molecular function (FDR = 2.79e-09), regulation of transcription by RNA polymerase II (FDR = 1.81e-09), and epilepsy (FDR = 0.0388). Three interactions were associated with cholinergic synapses (BAK1-BCL2L1-BCL2L11-BCL2*, NR2F2-ESR1-SRC-EGFR-PTPN11-PIK3R1*-PIK3CD*-PIK3CA*-PIK3 R2*, PRKCA*-PFKFB2), and two with glutamatergic synapses (GRIN2A*-GRIN1*, PRKCA*-PFKFB2), where * indicates the genes related with synaptic pathways.
Protein-protein interaction network analysis. We used STRING selecting experiments and co-expression with the highest confidence value (0.9).
Genome-wide association studies enrichment of differential mCpHs
Genome-wide association studies enrichment analysis (Table 2 and Supplementary Table 5) of all annotated genes from the OUD-associated differential mCpHs identified enrichment for two opioid gene sets—Methadone dose in opioid dependence (p-value = 3.00E-03; SORCS1, GRK5, E2F7, SPRY2, TRIB2, and BCL11A), and Opioid dependence (p-value = 6.18E-03; CNIH3, TACC2). In addition, we observed enrichment for Cannabis use (p-value = 8.83E-09), Cognitive decline rate in late mild cognitive impairment (p-value = 4.40E-07), Neuroticism (p-value = 6.03E-05), and Smoking initiation (p-value = 8.11E-05).
Genome-wide association studies enrichment analysis of annotated genes from the OUD-associated differential mCpHs.
GeneSet
N
n
P-value
Adjusted p
Cannabis use
10
7
8.83E-09
1.99E-06
Schizophrenia
827
60
3.52E-07
5.82E-05
Cognitive decline rate in late mild cognitive impairment
122
18
4.40E-07
6.65E-05
Cognitive ability, years of educational attainment or schizophrenia (pleiotropy)
197
23
8.86E-07
1.24E-04
Intelligence (MTAG)
313
30
1.51E-06
1.83E-04
Autism spectrum disorder, attention deficit-hyperactivity disorder, bipolar disorder, major depressive disorder, and schizophrenia (combined)
46
10
4.69E-06
4.06E-04
General cognitive ability
256
24
2.34E-05
1.34E-03
Systolic blood pressure x smoking status (current vs. non-current) interaction (2 df test)
111
14
5.01E-05
2.27E-03
Neuroticism
99
13
6.03E-05
2.61E-03
Smoking initiation
51
9
8.11E-05
3.22E-03
Diastolic blood pressure x smoking status (current vs. non-current) interaction (2 df test)
107
13
1.35E-04
4.72E-03
Cognitive performance
84
11
2.24E-04
6.36E-03
Bitter alcoholic beverage consumption
91
11
4.54E-04
1.14E-02
Feeling tense
19
5
4.72E-04
1.14E-02
Systolic blood pressure x smoking status (ever vs. never) interaction (2 df test)
95
11
6.57E-04
1.47E-02
Major depressive disorder
210
18
6.90E-04
1.49E-02
Neuroticism
149
14
1.09E-03
2.04E-02
General factor of neuroticism
104
11
1.40E-03
2.37E-02
Bipolar I disorder
90
10
1.56E-03
2.54E-02
Smoking status (ever vs. never smokers)
190
16
1.60E-03
2.57E-02
Methadone dose in opioid dependence
40
6
3.00E-03
3.97E-02
Alcohol dependence
54
7
3.24E-03
4.14E-02
Bipolar disorder
656
38
3.39E-03
4.21E-02
Cognitive function
85
9
3.67E-03
4.37E-02
Feeling fed-up
30
5
4.16E-03
4.72E-02
Cognitive function
85
11
1.80E-03
1.94E-02
Experiencing mood swings
40
7
2.15E-03
2.25E-02
Bipolar I disorder
90
11
2.86E-03
2.76E-02
Nicotine dependence
24
5
4.23E-03
3.49E-02
Cognitive ability
46
7
4.85E-03
3.76E-02
Feeling hurt
25
5
5.09E-03
3.80E-02
Bipolar disorder with mood-incongruent psychosis
16
4
5.26E-03
3.88E-02
Feeling worry
48
7
6.16E-03
4.42E-02
Opioid dependence
3
2
6.18E-03
4.42E-02
Pulse pressure x alcohol consumption (light vs. heavy) interaction (2 df test)
17
4
6.63E-03
4.70E-02
HDL cholesterol levels x alcohol consumption (regular vs. non-regular drinkers) interaction (2 df)
129
13
6.73E-03
4.72E-02
Significant enrichment GeneSet involved in drug use, psychiatric disorders, and cognitive function. N means the number of genes in the GeneSet and n represents the number of genes detected in the input.
Gene-drug interaction analysis
Drug interaction analysis showed 3,420 interactions (Supplementary Table 6) between the annotated genes with differential mCpHs and drugs. We evaluated the interactions with 15 opioids: Apomorphine, codeine, diacetylmorphine, hydrocodone, methadone, morphine, oxycodone, oxymorphone, tramadol, methadone, propoxyphene, fentanyl, hydromorphone, heroin, levorphanol, meptazinol, naloxone, pethidine, buprenorphine, and pentazocine. As a result, we found 12 annotated genes in our differential mCpHs that interacts with opioids (Figure 5), ATXN2, CBFB, CDH2, GAD1, GRIN1, KCNH2, MAPT, OPRK1, RAD52, RARA, GRP, and TAOK3. We also identified genes with known interactions with pergolide mesylate (KCNH2), lofexidine (ADRA2B), and tizanidine (ADRA2B). Finally, we evaluated drug-gene interactions for genes associated with two specific pathways detected in the functional annotation analysis: glutamatergic synapses (KEGG, Supplementary Figure 2), and cholinergic synapses (KEGG, Supplementary Figure 3). We further evaluated the location of differential mCpH in the genes with described interaction with opioids, and observed that ATXN2, CBFB, CDH2, GRIN1, KCNH2, OPRK1, RARA, and GRP (one mCpH in each gene), MAPT and RAD52 (two mCpHs each) were located inside the gene, GAD1 and GRP (one mCpH) were located upstream, and TAOK3 was annotated in an intergenic region.
Gene-drug interactions with opioids. The figure shows all genes that interact with 15 selected opioids. OPRK1 interacts with most opioid-related drugs.
Discussion
Here, we examined the role of non-CpG methylation in OUD, which to date had not been investigated. Evaluating genome-wide differential mCpHs associated with OUD in human postmortem OFC tissue, revealed 2,351 differential mCpHs enriched for glutamatergic synapses and cholinergic synapses. Gene-drug interaction analysis showed ten annotated genes with differential mCpHs interacting with opioid-related drugs, including lofexidine and tizanidine, both currently used for OUD treatment. These results highlight neuronal mCpH as an important regulatory mechanism in the pathophysiology of OUD.
Interestingly, several of the differential mCpHs mapped to genes previously described in OUD GWAS studies: OPRK1, APBB2, and CNIH3. OPRK1 encodes a kappa opioid receptor and is known to be associated with morphine (38) and alcohol dependence (39). In our analysis, this gene was hypomethylated and showed interactions with several opiates. A recent study reported a role of OPRK1 genetic variants in pain modulation (40). Further, a previous candidate gene study found associations of OPRK1 variants with methadone dosage in a Chinese population (n = 801) (41). APBB2 overexpression is associated with β-amyloid (42) and has been associated with Alzheimer’s disease (43). GWAS studies have suggested an association of APBB2 (p-value = 1.8 × 10–9) with OUD in African Americans (AA) of two independent samples of the SAGE cohort (n = 3,318 and n = 1,311) (44). More recently, variants at APBB2 were also associated amphetamine use in Taiwanese individuals with OUD receiving methadone treatment (n = 344) (45). CNIH3 is a glutamate receptor auxiliary protein that regulates trafficking properties in glutamate receptors. This gene has been associated with daily opioid use in an Australian cohort (n = 1,167), a finding replicated in the Yale Penn cohort (ndaily = 643; nlifetime = 157) (46).
We also identified mCpHs in other genes enriched for opioid-related gene sets: SORCS1, GRK5, E2F7, SPRY2, TRIB2, BCL11A, CNIH3, and TACC2. SORCS1 is a regulator of the synaptic trafficking required for glutamatergic synaptic transmission in mice brain (47) being associated with Alzheimer’s disease risk in a candidate gene study (48). GRK5 is a kinase that phosphorylates G-protein-coupled receptors, including dopamine and opioid receptors. An up-regulation of GRK5 has been described during initial morphine treatment in male Sprague–Dawlay rats and suggested to play a role in the regulation of opioid signaling (49). Genetic variants at the GRK5 gene have been associated with the regulation in methadone dosage in heroin dependence (50). SPRY2 is an intracellular inhibitor of growth factors involved in stress sensitivity (51). TRIB2 is a member of the tribble family; risk alleles mapping to this gene have been also associated with methadone dose in an African Americans population (52). BCL11A is a zinc-finger protein, in which its deletion is associated with higher levels of fetal hemoglobin (HbF). Although we did not observe expression changes of HbF, we did observe a differential gene expression in the HBB associated with OUD in our previous study (15).
Here, we also identified an OUD-associated differential mCpH site near the dopamine neurotrophic factor CDNF, involved in dopaminergic neuron differentiation and shown to have neuroprotective and neurorestorative effects in these neurons (53). Further, two genes with differential mCpHs in OUD, RERE, and CFAP77, were previously identified in an epigenome-wide association study from our group evaluating differential 5mC associated with opioid dependence in whole blood samples from European–American women (n = 220) (6). RERE is highly expressed in the brain and plays a role in brain development, retinoic acid signaling, and gene regulation (54). CFAP77 is involved in cell projection. Despite the unclear role of these two genes in OUD, the replication of these epigenetic markers across tissues suggests a promising important role in the pathophysiology of OUD.
Functional annotations of the differential mCpHs using KEGG highlighted two enriched pathways—glutamatergic and cholinergic synapses. The glutamatergic system has been implicated in opioid-induced hyperalgesia, where patients undergoing chronic opioid therapy endure acute pain and tolerance (55). In addition, studies have suggested that morphine exposure may decrease glutamate transporter expression (56). There is also evidence that glutamate transporters interact with opioid receptors (56). The link between glutamatergic signaling and opioid addiction could be also due to its role in learning and memory formation, where glutamatergic circuits may be recruited to create and maintain memories of reward, aversion, and reinforcement processes (57). The cholinergic system, implicated in several brain functions, including sensory processing, sleep, memory, learning, and attention (58, 59), has been suggested as a promising target for OUD treatment (60). Furthermore, opioids modulate cholinesterase activity, an enzyme involved in regulating cholinergic response (61). Moreover, we observed enrichment of the Wnt signaling pathway, previously involved in opioid-related withdrawal in mice (62). Our results show that dysregulated mCpHs sites target genes involved in both cholinergic and glutamatergic systems as well as Wnt signaling, all previously implicated in opioid-related traits.
We also observed several OUD-associated differential mCpHs genes involved in alternative splicing as predicted with the Protein-Protein Interaction analysis. A recent study assessing the interaction of alternative splicing in opioid-induced hyperalgesia in the nucleus accumbens (NAc) and trigeminal ganglia of mice identified four genes with more than four nodes connections in the NAc: Grin1, Src, Dnm1, and Pxn (63), three of which were detected in our analysis with differential mCpHs, GRIN1, SRC, and DNM1. The GRIN1 gene was also identified with known interactions with opiates. A study evaluating mRNA alternative splicing in the human postmortem dlPFC, NAc, and midbrain associated with OUD identified eight differential spliced genes in morphine addiction (64), one of which was identified in our analysis: GABBR1. In addition, they observed five spliced genes in the dlPFC, NAc, and ventral midbrain (64), including two genes with differential mCpHs found in our work: HERC1 and BIN1. These results support an important gene regulatory role of OUD-linked differential mCpHs.
Genome-wide association studies enrichment of annotated genes showing differential mCpHs showed associations with opioid dependence as well as other substance use such as tobacco, alcohol, and cannabis. Tobacco smoking is considered a major risk factor for OUD development (65). Further, studies have shown co-morbidities between OUD and other substance use disorders, including cannabis (66), nicotine (67), and alcohol (66, 68). We also observed enrichment for cognitive decline and impairment. Cognitive deficit is one of the withdrawal symptoms of OUD (69). Our findings suggest that epigenetic mechanisms may underlie the co-occurrence of SUDs. However, studies in larger and better-powered samples may help disentangle potential shared epigenetic mechanisms in SUDs comorbidities.
Gene-drug interaction analysis showed interactions of annotated differential mCpH genes with opioid-related drugs, including lofexidine (helps with physical symptoms of opioid withdrawal), tizanidine (used to treat withdrawal symptoms in heroin users), and pergolide mesylate (used to treat cocaine dependence). The genes interacting with opioids included GRIN1 and OPRK1, described above. GAD1 is a glutamic acid decarboxylase related to glutamate-dependent acid resistance; candidate gene studies reported genetic variants associated with heroin addiction (70). KCNH2, a Potassium Voltage-Gated Channel Subfamily H Member 2, component of potassium channels, was also detected in our study associated with pergolide mesylate, a drug used to treat cocaine dependence (71). Further, KCNH2 polymorphisms were shown as a potential risk factor associated with QT prolongation (interval between the ventricle depolarization and repolarization) under low doses of methadone treatment in a French cohort (72). We also observed the interaction between ADRA2B and tizanidine, lofexidine, and dopamine. ADRA2B is involved in the thermoregulatory effects induced by the opiate meperidine in mice (73). Interestingly, both tizanidine (74) and lofexidine has been used to treat symptoms related to opioid withdrawal (75).
Comparing the non-CpG 5mC results with CpG 5mC described in Rompala et al. (15), we observed more OUD-associated differential changes in mCpHs than mCpGs. Further, mCpHs enrichment showed a higher number of significantly enriched pathways, including KEGG signaling pathways, glutamatergic and cholinergic synapses. Drug interaction analysis also showed a higher number of annotated genes with differential mCpHs than mCpGs associated with opioids, suggesting that 5mC at non-CpG sites in OFC neurons may be more functionally relevant than 5mC at CpG sites in the context of OUD. Our results highlight the importance of studying 5mC at non-CpG sites. We also evaluated the differential gene expression in association with OUD in OFC previously described in Rompala et al. (15) on the same cohort as in the current study; however, no significant OUD-associated differential expression was observed of the genes identified here with differential mCpH (p-value range = 0.091–0.987; Supplementary Table 7).
We further evaluated whether the identified differential mCpHs would be targeted for activation of gene expression through the prediction of transcriptional factor binding site (TFBS) and enrichment of H3K27Ac markers (Table 3). We observed that seven of the differential mCpHs identified here that was located upstream of the annotated gene showed predicted TFBS. A total of five mCpHs were located in exon, of which four have predicted TFBS. A total of eight were detected inside introns, of which five showed predicted TFBS sites. The only identified differential mCpH in the intergenic region was also detected with predicted TFBS. Further, 14 of these mCpH sites were observed with H3K27Ac marker, suggesting a potential role in gene regulation (Table 3).
Detailed information about transcription factor binding site (TBFS) and H3K27Ac mark.
The table includes information about TFBS and H3K27Ac mark for the sites discussed.
This study has limitations relevant to interpreting the data. We studied a small sample of subjects that included only 12 individuals diagnosed with OUD, of which we only evaluated differences in methylation levels of males, most of European ancestry. Furthermore, the potential effects of comorbidities often associated with OUD and other SUDs were not fully addressed due to lack of power due to the limited sample size. Future work aims to increase the sample size and evaluate the effects of ancestry, sex, as well as the influence of comorbidities in a better-powered, more diverse cohort.
Conclusion
In conclusion, this is the first examination of non-CpG methylation in OUD, identifying 2,351 differential mCpHs that included well-known opioid-related genes, such as the OPRK1 gene. We identified cholinergic and glutamatergic synapses among the functional pathways associated with the differential mCpHs. Drug interaction analysis also showed several genes interacting with opioids, including the OPRK1, GRIN1, and KCNH2, as well as gene interactions with drugs used for the treatment of opioid-related withdrawal symptoms. Our findings reveal mCpHs as a crucial regulatory mechanism in OUD and may help shed light on future therapeutic and preventive strategies for treating this disorder.
Members of the Traumatic Stress Brain Research Group
Victor E. Alvarez, David Benedek, Alicia Che, Dianne A. Cruz, David A. Davis, Matthew J. Girgenti, Ellen Hoffman, Paul E. Holtzheimer, Bertrand R. Huber, Alfred Kaye, John H. Krystal, Adam T. Labadorf, Terence M. Keane, Mark W. Logue, Ann McKee, Brian Marx, Deborah Mash, Mark W. Miller, Crystal Noller, JM-O, William K. Scott, Paula Schnurr, Thor Stein, Robert Ursano, Douglas E. Williamson, Erika J. Wolf, and Keith A. Young.
Data availability statement
The original contributions presented in this study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.
Ethics statement
The studies involving human participants were reviewed and approved by the IRB Committee from the Department of Veteran’s Affairs. The patients/participants provided their written informed consent to participate in this study.
Author contributions
SN was responsible for data analysis, interpretation, and manuscript writing. GR and DN-R contributed to the data analyses. YH contributed to data generation and study design. TSBRG contributed to data collection. JM-O contributed to data generation and interpretation, study design, and supervision of manuscript writing. All authors contributed to the revision of the manuscript and approved its final version.
Funding
This study was supported by the U.S. Department of Veterans Affairs via 1IK2CX002095-01A1 (JM-O) and NIDA R21DA050160 (JM-O).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpsyt.2022.1078894/full#supplementary-material
Abbreviations
mCpH, methylation in non-CpG site; 5mC, methylation in CpG site; 5hmC, hydroxymethylation in CpG site.
https://www.dgidb.org
ReferencesLydenJBinswangerI. The united states opioid epidemic.Semin Perinatol. (2019) 43:123–31.LancetT.A time of crisis for the opioid epidemic in the USA.Lancet (London, England). (2021) 398:277.DeakJZhouHGalimbertiMLeveyDWendtFSanchez-RoigeSGenome-wide association study in individuals of European and African ancestry and multi-trait analysis of opioid use disorder identifies 19 independent genome-wide significant risk loci.Mol Psychiatry. (2022) 27:3970–9. 10.1038/s41380-022-01709-135879402ZhouHRentschCChengZKemberRNunezYShervaRAssociation of OPRM1 functional coding variant with opioid use disorder: a genome-wide association study.JAMA Psychiatry. (2020) 77:1072–80.JangHShinWLeeJDoJ. CpG and non-CpG methylation in epigenetic gene regulation and brain function.Genes. (2017) 8:6.Montalvo-OrtizJChengZKranzlerHZhangHGelernterJ. Genomewide study of epigenetic biomarkers of opioid dependence in European- American women.Sci Rep. (2019) 9:4660.Marie-ClaireCCrettolSCagnardNBlochVMoulySLaplancheJVariability of response to methadone: genome-wide DNA methylation analysis in two independent cohorts.Epigenomics. (2016) 8:181–95. 10.2217/epi.15.11026792095BorrelliKWachmanEBeierleJTaglauerEJainMBryantCEffect of prenatal opioid exposure on the human placental methylome.Biomedicines. (2022) 10:5. 10.3390/biomedicines1005115035625888KozlenkovAJaffeATimashpolskyAApontesPRudchenkoSBarbuMDNA methylation profiling of human prefrontal cortex neurons in heroin users shows significant difference between genomic contexts of hyper- and hypomethylation and a younger epigenetic age.Genes. (2017) 8:6. 10.3390/genes806015228556790HuYSalmeronBKrasnovaIGuHLuHBonciACompulsive drug use is associated with imbalance of orbitofrontal- and prelimbic-striatal circuits in punishment-resistant individuals.Proc Natl Acad Sci USA. (2019) 116:9066–71.SchoenbaumGChangCLucantonioFTakahashiY. Thinking outside the box: orbitofrontal cortex, imagination, and how we can treat addiction.Neuropsychopharmacology. (2016) 41:2966–76. 10.1038/npp.2016.14727510424SchoenbaumGRoeschMStalnakerT. Orbitofrontal cortex, decision-making and drug addiction.Trends Neurosci. (2006) 29:116–24.ShuCSosnowskiDTaoRDeep-SoboslayAKleinmanJHydeTEpigenome-wide study of brain DNA methylation following acute opioid intoxication.Drug Alcohol Depend. (2021) 221:108658. 10.1016/j.drugalcdep.2021.10865833667780LiuADaiYMendezEHuRFriesGNajeraKGenome-wide correlation of DNA methylation and gene expression in postmortem brain tissues of opioid use disorder patients.Int J Neuropsychopharmacol. (2021) 24:879–91. 10.1093/ijnp/pyab04334214162RompalaGNagamatsuSMartínez-MagañaJWangJGirgentiMKrystalJProfiling neuronal methylome and hydroxymethylome of opioid use disorder in the human orbitofrontal cortex.medRxiv[preprint]. (2022). 10.1101/2022.09.09.22279769BoothMMarsicoGBachmanMBeraldiDBalasubramanianS. Quantitative sequencing of 5-formylcytosine in DNA at single-base resolution.Nat Chem. (2014) 6:435–40.LutzPGrossJDhirSMaussionGYangJBramoulleAEpigenetic regulation of the kappa opioid receptor by child abuse.Biol Psychiatry. (2018) 84:751–61. 10.1016/j.biopsych.2017.07.01228886759ZhangHHermanAKranzlerHAntonRSimenAGelernterJ. Hypermethylation of OPRM1 promoter region in European Americans with alcohol dependence.J Hum Genet. (2012) 57:670–5. 10.1038/jhg.2012.9822914673SunLZhaoJGuXLiangLWuSMoKNerve injury-induced epigenetic silencing of opioid receptors controlled by DNMT3a in primary afferent neurons.Pain. (2017) 158:1153–65.Perzel MandellKEaglesNWiltonRPriceASemickSCollado-TorresLGenome-wide sequencing-based identification of methylation quantitative trait loci and their role in schizophrenia risk.Nat Commun. (2021) 12:5251. 10.1038/s41467-021-25517-334475392EllisSGuptaSMoesAWestAArkingD. Exaggerated CpH methylation in the autism-affected brain.Mol Autism. (2017) 8:6. 10.1186/s13229-017-0119-y28316770FriedmanMHuberBBradyCUrsanoRBenedekDKowallNVA’s national PTSD brain bank: a national resource for research.Curr Psychiatry Rep. (2017) 19:73. 10.1007/s11920-017-0822-628840457KruegerFAndrewsS. Bismark: a flexible aligner and methylation caller for bisulfite-seq applications.Bioinformatics. (2011) 27:1571–2. 10.1093/bioinformatics/btr16721493656AkalinAKormakssonMLiSGarrett-BakelmanFFigueroaMMelnickAmethylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles.Genome Biol. (2012) 13:R87. 10.1186/gb-2012-13-10-r8723034086AkalinAFrankeVVlahovicekKMasonCSchubelerD. Genomation: a toolkit to summarize, annotate and visualize genomic intervals.Bioinformatics. (2015) 31:1127–9. 10.1093/bioinformatics/btu77525417204SmedleyDHaiderSBallesterBHollandRLondonDThorissonGBioMart–biological queries made easy.BMC Genomics. (2009) 10:22. 10.1186/1471-2164-10-2219144180KentWSugnetCFureyTRoskinKPringleTZahlerAThe human genome browser at UCSC.Genome Res. (2002) 12:996–1006.ChenETanCKouYDuanQWangZMeirellesGEnrichr: interactive and collaborative HTML5 gene list enrichment analysis tool.BMC Bioinform. (2013) 14:128. 10.1186/1471-2105-14-12823586463SzklarczykDGableALyonDJungeAWyderSHuerta-CepasJSTRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets.Nucleic Acids Res. (2019) 47:D607–13. 10.1093/nar/gky113130476243McKenzieAWangMHaubergMFullardJKozlenkovAKeenanABrain cell type specific gene expression and co-expression network architectures.Sci Rep. (2018) 8:8868.WatanabeKTaskesenEvan BochovenAPosthumaD. Functional mapping and annotation of genetic associations with FUMA.Nat Commun. (2017) 8:1826.HuangDShermanBTanQCollinsJAlvordWRoayaeiJThe DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists.Genome Biol. (2007) 8:R183. 10.1186/gb-2007-8-9-r18317784955Castro-MondragonJRiudavets-PuigRRauluseviciuteILemmaRTurchiLBlanc-MathieuRJASPAR 2022: the 9th release of the open-access database of transcription factor binding profiles.Nucleic Acids Res. (2022) 50:D165–73. 10.1093/nar/gkab111334850907FreshourSKiwalaSCottoKCoffmanAMcMichaelJSongJIntegration of the drug-gene interaction database (DGIdb 4.0) with open crowdsource efforts.Nucleic Acids Res. (2021) 49:D1144–51. 10.1093/nar/gkaa108433237278ShannonPMarkielAOzierOBaligaNWangJRamageDCytoscape: a software environment for integrated models of biomolecular interaction networks.Genome Res. (2003) 13:2498–504. 10.1101/gr.123930314597658CristRReinerBBerrettiniW. A review of opioid addiction genetics.Curr Opin Psychol. (2019) 27:31–5.YuferovVButelmanERandesiMOttJKreekM. Analyses of polymorphisms of intron 2 of OPRK1 (kappa-opioid receptor gene) in association with opioid and cocaine dependence diagnoses in an African-American population.Neurosci Lett. (2022) 768:136364. 10.1016/j.neulet.2021.13636434843875GerraGLeonardiCCorteseED’AmoreALucchiniAStrepparolaGHuman kappa opioid receptor gene (OPRK1) polymorphism is associated with opiate addiction.Am J Med Genet B Neuropsychiatr Genet. (2007) 144B:771–5.ZhangHKranzlerHYangBLuoXGelernterJ. The OPRD1 and OPRK1 loci in alcohol or drug dependence: OPRD1 variation modulates substance dependence risk.Mol Psychiatry. (2008) 13:531–43. 10.1038/sj.mp.400203517622222HoKWallaceMStaudRFillingimR. OPRM1, OPRK1, and COMT genetic polymorphisms associated with opioid effects on experimental pain: a randomized, double-blind, placebo-controlled study.Pharmacogenomics J. (2020) 20:471–81. 10.1038/s41397-019-0131-z31806881ZhangQShiMTangHZhongHLuX. Opioid receptor 1 single nucleotide polymorphisms were associated with the methadone dosage.Genet Test Mol Biomarkers. (2020) 24:17–23.ChangYTescoGJeongWLindsleyLEckmanEEckmanCGeneration of the beta-amyloid peptide and the amyloid precursor protein C-terminal fragment gamma are potentiated by FE65L1.J Biol Chem. (2003) 278:51100–7. 10.1074/jbc.M30956120014527950LiYHollingworthPMoorePFoyCArcherNPowellJGenetic association of the APP binding protein 2 gene (APBB2) with late onset Alzheimer disease.Hum Mutat. (2005) 25:270–7. 10.1002/humu.2013815714520GelernterJKranzlerHShervaRKoestererRAlmasyLZhaoHGenome-wide association study of opioid dependence: multiple associations mapped to calcium and potassium pathways.Biol Psychiatry. (2014) 76:66–74. 10.1016/j.biopsych.2013.08.03424143882LiuCFangCLiuTKuoHLiuSWangSAPBB2 is associated with amphetamine use and plasma beta-amyloids in patients receiving methadone maintenance treatment.Prog Neuropsychopharmacol Biol Psychiatry. (2018) 83:92–8. 10.1016/j.pnpbp.2018.01.00829330135NelsonEAgrawalAHeathABogdanRShervaRZhangBEvidence of CNIH3 involvement in opioid dependence.Mol Psychiatry. (2016) 21:608–14.SavasJRibeiroLWierdaKWrightRDeNardo-WilkeLRiceHThe sorting receptor SorCS1 regulates trafficking of neurexin and AMPA receptors.Neuron. (2015) 87:764–80. 10.1016/j.neuron.2015.08.00726291160ReitzCTokuhiroSClarkLConradCVonsattelJHazratiLSORCS1 alters amyloid precursor protein processing and variants may increase Alzheimer’s disease risk.Ann Neurol. (2011) 69:47–64.FanXZhangJZhangXYueWMaL. Acute and chronic morphine treatments and morphine withdrawal differentially regulate GRK2 and GRK5 gene expression in rat brain.Neuropharmacology. (2002) 43:809–16. 10.1016/s0028-3908(02)00147-8WangSChungRKuoHLiuTFangCLiuSGRK5 is associated with the regulation of methadone dosage in heroin dependence.Int J Neuropsychopharmacol. (2018) 21:910–7. 10.1093/ijnp/pyy06630060048DowALinTChartoffEPotterDMcPhieDVan’t VeerASprouty2 in the dorsal hippocampus regulates neurogenesis and stress responsiveness in rats.PLoS One. (2015) 10:e0120693. 10.1371/journal.pone.012069325822989ChawarCHillmerASangerSD’EliaAPanesarBGuanLA systematic review of GWAS identified SNPs associated with outcomes of medications for opioid use disorder.Addict Sci Clin Pract. (2021) 16:70. 10.1186/s13722-021-00278-y34838141LindholmPVoutilainenMLaurenJPeranenJLeppanenVAndressooJNovel neurotrophic factor CDNF protects and rescues midbrain dopamine neurons in vivo.Nature. (2007) 448:73–7. 10.1038/nature0595717611540FregeauBKimBHernandez-GarciaAJordanVChoMSchnurRDe novo mutations of RERE cause a genetic syndrome with features that overlap those associated with proximal 1p36 deletions.Am J Hum Genet. (2016) 98:963–70. 10.1016/j.ajhg.2016.03.00227087320ZhouHChenSChenHPanH. Opioid-induced long-term potentiation in the spinal cord is a presynaptic event.J Neurosci. (2010) 30:4460–6. 10.1523/JNEUROSCI.5857-09.201020335482AlasmariFSariDAlhaddadHAl-RejaieSSariY. Interactive role of acid sensing ion channels and glutamatergic system in opioid dependence.Neurosci Biobehav Rev. (2022) 135:104581. 10.1016/j.neubiorev.2022.10458135181397PetersJDe VriesT. Glutamate mechanisms underlying opiate memories.Cold Spring Harb Perspect Med. (2012) 2:a012088.LiXYuBSunQZhangYRenMZhangXGeneration of a whole-brain atlas for the cholinergic system and mesoscopic projectome analysis of basal forebrain cholinergic neurons.Proc Natl Acad Sci USA. (2018) 115:415–20. 10.1073/pnas.170360111529259118HampelHMesulamMCuelloAFarlowMGiacobiniEGrossbergGThe cholinergic system in the pathophysiology and treatment of Alzheimer’s disease.Brain. (2018) 141:1917–33.JensenKDeVitoEYipSCarrollKSofuogluM. The cholinergic system as a treatment target for opioid use disorder.CNS Drugs. (2018) 32:981–96.AramjooHRiahi-ZanjaniBFarkhondehTForouzanfarFSadeghiM. Modulatory effect of opioid administration on the activity of cholinesterase enzyme: a systematic review of mice/rat models.Environ Sci Pollut Res Int. (2021) 28:52675–88. 10.1007/s11356-021-16044-134453251WuMLiZLiangLMaPCuiDChenPWnt signaling contributes to withdrawal symptoms from opioid receptor activation induced by morphine exposure or chronic inflammation.Pain. (2020) 161:532–44.ZhangPPerezOSoutheyBSweedlerJPradhanARodriguez-ZasS. Alternative splicing mechanisms underlying opioid-induced hyperalgesia.Genes. (2021) 12:10.HuggettSIkedaAMcGearyJKaunKPalmerR. Opioid use disorder and alternative mRNA splicing in reward circuitry.Genes. (2022) 13:6. 10.3390/genes1306104535741807RajabiADehghaniMShojaeiAFarjamMMotevalianS. Association between tobacco smoking and opioid use: a meta-analysis.Addict Behav. (2019) 92:225–35.BahorikASatreDKline-SimonAWeisnerCCampbellC. Alcohol, cannabis, and opioid use disorders, and disease burden in an integrated health care system.J Addict Med. (2017) 11:3–9. 10.1097/ADM.000000000000026027610582LichensteinSZakiniaeizYYipSGarrisonK. Mechanisms and clinical features of co-occurring opioid and nicotine use.Curr Addict Rep. (2019) 6:114–25. 10.1007/s40429-019-00245-332864292MintzCPresnallNXuKHartzSSahrmannJBierutLAn examination between treatment type and treatment retention in persons with opioid and co-occurring alcohol use disorders.Drug Alcohol Depend. (2021) 226:108886. 10.1016/j.drugalcdep.2021.10888634245997RapeliPKivisaariRAuttiTKahkonenSPuuskariVJokelaOCognitive function during early abstinence from opioid dependence: a comparison to age, gender, and verbal intelligence matched controls.BMC Psychiatry. (2006) 6:9. 10.1186/1471-244X-6-916504127ShiYLiYZhangJXiaoYYanPZhuY. GAD1 but not GAD2 polymorphisms are associated with heroin addiction phenotypes.Neurosci Lett. (2020) 717:134704.MalcolmRHuttoBPhillipsJBallengerJ. Pergolide mesylate treatment of cocaine withdrawal.J Clin Psychiatry. (1991) 52:39–40.HajjAKsoudaKPeoc’hKCurisEMessaliADeveauxLKCNH2 polymorphism and methadone dosage interact to enhance QT duration.Drug Alcohol Depend. (2014) 141:34–8. 10.1016/j.drugalcdep.2014.04.02724875677HockerJParisAScholzJTonnerPNielsenMBeinB. Differential effects of alpha 2-adrenoceptors in the modulation of the thermoregulatory response in mice induced by meperidine.Anesthesiology. (2008) 109:95–100. 10.1097/ALN.0b013e31817c02fc18580178PinelliATrivulzioSSpeziaR. Effects of tizanidine administration on precipitated opioid withdrawal signs in rats.Drug Alcohol Depend. (1998) 50:81–8. 10.1016/s0376-8716(98)00010-6SrivastavaAMarianiJLevinF. New directions in the treatment of opioid withdrawal.Lancet. (2020) 395:1938–48.