Alcohol and nicotine codependence-associated DNA methylation changes in promoter regions of addiction-related genes

Altered DNA methylation in addiction-related genes may modify the susceptibility to alcohol or drug dependence (AD or ND). We profiled peripheral blood DNA methylation levels of 384 CpGs in promoter regions of 82 addiction-related genes in 256 African Americans (AAs) (117 cases with AD-ND codependence and 139 controls) and 196 European Americans (103 cases with AD-ND codependence and 93 controls) using Illumina’s GoldenGate DNA methylation array assays. AD-ND codependence-associated DNA methylation changes were analyzed using linear mixed-effects models with consideration of batch effects and covariates age, sex, and ancestry proportions. Seventy CpGs (in 41 genes) showed nominally significant associations (P < 0.05) with AD-ND codependence in both AAs and EAs. One CpG (HTR2B cg27531267) was hypomethylated in AA cases (P = 7.2 × 10−5), while 17 CpGs in 16 genes (including HTR2B cg27531267) were hypermethylated in EA cases (5.6 × 10−9 ≤ P ≤ 9.5 × 10−5). Nevertheless, 13 single nucleotide polymorphisms (SNPs) nearby HTR2B cg27531267 and the interaction of these SNPs and cg27531267 did not show significant effects on AD-ND codependence in either AAs or EAs. Our study demonstrated that DNA methylation changes in addiction-related genes could be potential biomarkers for AD-ND co-dependence. Future studies need to explore whether DNA methylation alterations influence the risk of AD-ND codependence or the other way around.


No significant interactive effects of HTR2B cg27531267 and nearby SNPs on AD-ND codependence.
Considering that genetic variants may either have a direct or an indirect (e.g., via altering DNA methylation patterns) influence on disease risk, we further analyzed the effects of 13 SNPs around HTR2B cg27531267 (50 kb up-or downstream of cg27531267) and the interaction of these SNPs with HTR2B cg27531267 on the susceptibility to AD-ND codependence. All 13 SNPs were located in a tight linkage disequilibrium (LD) block with average R2 of 0.48 in EAs, while the LD pattern was similar but less tight in AAs with average R2 of 0.28 ( Figure S1). As shown in Table 2, none of the 13 SNPs was significantly associated with AD-ND codependence in either AAs or EAs. Moreover, the interactive effect of HTR2B cg27531267 and nearby SNPs on AD-ND codependence was not significant in either AAs or EAs. Additionally, the methylation level of HTR2B cg27531267 was not significantly correlated with the genotypes of 13 nearby SNPs in either AAs or EAs (P correlation > 0.05).
Potential function of AD-ND codependence-associated CpGs. The  CpGs. Among the 17 differentially methylated CpGs (Table 1), 11 were predicted to be located in the core binding site of one or more transcription factors (Table 3). The UCSC Genome Browser was used to query DNase hypersensitivity sites (DHSs) and H3K27Ac marks in DNA sequences harboring differentially methylated CpGs. Among the 17 differentially methylated CpGs, 14 were located in DHSs and five were situated in DNA sequences that are associated with the H3K27Ac histone mark (Table 3). Although the most significant HTR2B cg27531267 was not predicted to be located in TFBSs, it was mapped to DHSs and DNA sequences associated with histone protein modification mark H3K27Ac (Fig. 2).

Discussion
Prior evidence supports the involvement of many of the genes included in the present study in neurobiologic processes underlying drug reward and addiction. However, with a few exceptions (such as those coding for nicotinic acetylcholine receptors or alcohol metabolizing enzymes), most of these genes have not been previously found to harbor genetic variants that are associated at a genome-wide level with AD and/or ND. There are other mechanisms by which these loci could exert more substantial effects on AD and/or ND: namely, DNA methylation changes in addiction-related genes could confer vulnerability to AD and/or ND. In the present study, we identified methylation alterations in the promoter regions of a number of addiction-related genes in African Americans (AAs) and European Americans (EAs) with AD-ND codependence. We found that differentially methylated genes are involved in several critical pathways for AD and/or ND. A post hoc power analysis demonstrated that 196 EAs samples provided 80% statistical power for an effect size greater than 0.004, while 256 AA samples provided 80% statistical power to detect an effect size greater than 0.003, indicating that our sample size provided adequate statistical power to detect methylation changes. As summarized in Table 1, four of the differentially methylated genes were serotonergic (HTR2B, HTR2C, HTR3A, and SLC6A4), three were dopaminergic (DRD2, DRD4, and SLC6A3), two were GABAergic (GABRB1 and GABRB2), one was glutamatergic (GAD1), and three (OPRD1, OPRK1, and PENK) were opioidergic. Three other genes (RGS17, NCAM1, and MBD3), which do not belong to the above-listed neurotransmitter systems, also showed promoter DNA methylation changes in subjects with AD-ND codependence. RGS17 encodes a regulator of G-protein signaling (RGS) 58 , which can inactivate the G protein and rapidly switch off G-protein-coupled receptor signaling pathways 59 . Our candidate gene studies have demonstrated that variation in RGS17 is associated with a variety of different substance dependence disorders 60 . NCAM1 encodes a neural cell adhesion protein, a member of the immunoglobulin superfamily 61 . It is located in the DRD2-ANKK1-TTC12-NCAM1 gene cluster region, and SNPs or haplotypic variants in this region were found to be associated with dependence on alcohol, nicotine, or other drugs of abuse [62][63][64] . MBD3 encodes methyl-CpG binding domain protein 3, a nuclear protein that is potentially involved in chromatin remodeling and histone modifications 65,66 . Because genetic association studies did not reveal a strong effect (i.e., with genome-wide significance) of variants in the above genes on AD and/or ND, we postulate that either inherent DNA methylation of these genes results in AD and/or ND or long-term alcohol misuse or smoking leads to AD and/or ND through epigenetic modifications. It should be noted that these epigenetic   changes may not contribute to the risk of AD and ND simultaneously or to the same extent. Significant CpGs identified in the present study are associated with AD-ND codependence, which was considered to be a new phenotype, i.e., the co-occurring risk of AD and ND. Our findings suggest that the impact of DNA methylation of addiction-related genes on AD-ND codependence risk may be larger than that of genetic variants carried by these genes. Considering the possible correlation of methylation levels of CpGs and genotypes of nearby SNPs as reported in our previous studies 67 , we further investigated the effect of 13 SNPs within ± 50 kb of HTR2B cg27531267, as well as the interaction of cg27531267 with these SNPs on AD-ND codependence risk. Although HTR2B cg27531267 was differentially methylated in both AAs and EAs (albeit in opposite directions) with AD-ND codependence, no interactive effect of HTR2B cg27531267 with proximal SNPs on AD-ND codependence was observed in either population. Additionally, none of the 13 nearby SNPs was significantly associated with AD-ND codependence ( Table 2). This finding is consistent with the results from our previous GWAS research that variants in HTR2B were not significantly associated with either AD 26 or ND 37 .
The present study provides further evidence that DNA methylation changes within the regulatory (promoter) regions of addiction-related genes are associated with AD and/or ND, presumably because they may result in altered gene transcription. The rationale is that promoter DNA methylation may directly interfere with the binding of transcription factors (TFs) to the regulatory regions. Among the top 17 CpGs located in promoter regions of 16 genes (Table 1), 11 were predicted to be located in the core binding site of one or more TFs (Table 3). It is well known that chromatin structure mediates the interaction of TFs and DNA 68,69 . Our findings suggest that promoter DNA methylation may also modulate TF-DNA interactions and subsequently influence gene transcription. Moreover, methylated CpG falling within DNase hypersensitivity sites (DHSs) may impede the association of TFs to DNA, thus inhibiting the accessibility of chromatin. Because we did not extract RNA from blood samples (which were used only for genomic DNA extraction), we were unable to perform RT-qPCR to confirm gene expression changes caused by differentially methylated CpGs. Additionally, methylation of CpGs mapped to DNA sequences that are associated with histone marks (e.g., H3K27Ac mark, which is often found near active regulatory elements) may change histone protein epigenetic status, thus influencing the compact structure of chromatins. We also noticed that some of the 17 differentially methylated CpGs (including the most significant HTR2B cg27531267) were located in DHSs or DNA sequences that are associated with histone mark H3K27Ac. Taken together, subjects with altered DNA methylation in promoter regions of addiction-related genes may have  an increased risk of AD and/or ND. Note that, further studies are needed to replicate the above findings in independent samples and extend the findings to other racial/ethnic groups. The major limitation of the present study is that we only investigated the association of CpGs in promoter regions of a number of preselected addiction-related genes. For a more complete understanding of the epigenetic mechanism of AD-ND codependence, it will be necessary to use a high-resolution DNA methylation array (such as the Illumina MethylationEPIC BeadChip) assays or whole genome bisulfite sequencing (WGBS) to identify DNA methylomic changes. Another limitation is that we did not consider the relative proportion of different types of blood cells in the DNA methylation data analysis. Inter-individual differences in DNA methylation levels due to different blood cell composition may confound the findings. In future, when we have the high-density DNA methylation data (such as those generated by the Illumina MethylationEPIC BeadChip), we could use the method developed by Jaffe and Irizarry 70 to estimate the relative proportions of CD4 + and CD8 + T-cells, natural killer cells, monocytes, granulocytes, and B-cells in blood samples and then incorporate the cell proportion estimation into the data analysis. Although DNA methylation patterns in the peripheral blood may not reflect those in the brain, peripheral blood samples are easier to collect than brain tissue samples, and blood DNA methylation changes in regulatory regions of genes could be accessible biomarkers.
We are aware that there are more significant CpGs in EAs than AAs, even though the EA sample was smaller than the AA sample. Additionally, it is unknown why HTR2B cg27531267 showed an opposite methylation direction in AAs (hypomethylation) and EAs (hypermethylation) with AD-ND codependence. One possible explanation is that the DNA methylation status of CpGs can be influenced by genetic variation. Similar to many disease-associated SNPs that are specific to a certain population, CpG methylation levels may also be population-specific, leading to inconsistent results between EAs and AAs. For example, the top CpG cg27531267 showed a mean methylation level of 0.045 in AA controls, which was significantly higher than that of EA controls (β = 0.035, P = 1.8 × 10 −8 ). Among the 13 SNPs that are proximal to HTR2B cg27531267 and included in CpG-SNP interaction analysis, two were either not existent or rare in EAs (Table 2).
In summary, the present study examined DNA methylation alterations in promoter regions of addiction-related genes among individuals with AD-ND codependence. We identified both specific and shared DNA methylation changes in the two populations. The overlap of the differentially methylated promoter CpGs and TFBSs or DHSs and the location of differentially methylated CpGs in DNA sequences that are associated with specific histone marks (e.g., H3K27Ac) imply that promoter CpG methylation may modulate gene transcription and influence an individual's susceptibility to AD-ND codependence. Considering the reversibility of DNA methylation, the findings from the present study could provide the basis for effective pharmacotherapies for AD-ND co-dependence that target specific epigenetic marks in promoter regions of addiction-related genes.  126). The mean age of AA cases and AA controls was 42 ± 8 and 37 ± 14 years, respectively (P < 0.05). The mean age of EA cases and EA controls was 43 ± 12 and 37 ± 16 years, respectively (P > 0.05). Both cases and controls were chosen from a large sample of subjects recruited for studies of the genetics of substance dependence. Subjects were interviewed using an electronic version of the Semi-Structured Assessment for Drug Dependence and Alcoholism (SSADDA) 71 . Lifetime diagnoses for AD and ND codependence were made according to the criteria of the Diagnostic and Statistical Manual of Mental Disorders, 4th edition (DSM-IV) [American Psychiatric Association, 1994]. None of the subjects had a lifetime major psychotic disorder such as schizophrenia and bipolar disorder. Additionally, no control subjects (139 AAs and 93 EAs) were affected with alcohol or drug abuse or dependence. Subjects gave informed consent as approved by the institutional review board at each clinical site, and certificates of confidentiality were obtained from the National Institute on Drug Abuse and the National Institute on Alcohol Abuse and Alcoholism. All methods were carried out in accordance with relevant guidelines and regulations. All experimental protocols were approved by the Human Investigation Committee of the above three institutes.

Subjects
The self-reported genetic background of all subjects included in the present study was verified using a set of 41 ancestry informative markers (AIMs), including 36 short tandem repeat markers and five SNPs, implemented in the program STRUCTURE 72 . Subjects were defined as EAs if their ancestry proportion scores were less than 0.5; otherwise, they were considered AAs. To minimize the influence of subjects' genetic background on the results of the DNA methylation analysis, subjects' ancestry proportions were considered as a covariate in the differential CpG analysis.
DNA extraction and bisulfite modification. Genomic DNA was extracted from the peripheral blood of cases and controls using the PAXgene Blood DNA Kit (PreAnalytiX, Hombrechtikon, Switzerland). One microgram of genomic DNA was treated with the bisulfite reagent included in the EZ DNA Methylation Kit (Zymo Research, Orange, CA, USA). Unmethylated cytosines were converted to uracils while methylated cytosines remained unchanged 73 . Bisulfite-converted DNA samples were then used in the custom-designed Illumina GoldenGate DNA methylation assay.
After bisulfite conversion of genomic DNA, the remaining methylation assay steps were the same as those previously described 50 . Image processing and intensity data extraction were performed using the Illumina GenomeStudio TM Methylation Module v.1.0 Software. The background normalization algorithm was used to minimize background variation within the array by using built-in negative control signals. The methylation level (defined as β ) of each individual CpG site was estimated as the ratio of intensities between methylated and unmethylated alleles. The β value was calculated as β = [Max(Cy5,0)]/[Max (Cy3,0) + Max(Cy5,0) + α ]. A constant offset α (by default, α = 100) was added to the denominator to adjust β values when both methylated and unmethylated probe intensities were low. The β value ranges from 0 (i.e., completely unmethylated) to 1 (i.e., completely methylated).
Statistical and bioinformatics analyses. The statistical analysis of DNA methylation data was implemented in R (version 3.1.3) (http://www.r-project.org/). To identify differentially methylated CpGs in subjects with AD and ND codependence, we used the linear mixed-effects model handled by the R package CpGassoc, which was designed for analyzing the association between methylation at CpG sites across the genome and a SCieNTiFiC RepoRts | 7:41816 | DOI: 10.1038/srep41816 phenotype of interest 74 . In the linear mixed-effect model built-in CpGassoc, the methylation level of CpGs was the response variable, the status of AD-ND codependence and other covariates (age, sex, and ancestry proportions) were fitted via fixed effects, and batch factor (referring to methylation chips) was fitted as a random effect. A post hoc power analysis by R package SIMR 75 served to calculate the power for a linear mixed model based on the observed data structure. For both EAs and AAs, we evaluated the sample power by varying effect sizes and using 100 simulations.
The interactive effect of differentially methylated CpGs and nearby single nucleotide polymorphisms (SNPs) on AD-ND codependence risk was examined using the generalized linear mixed-effects model by the lme4 package in R (http://CRAN.R-project.org/package= lme4). Although the correlation of CpGs and SNPs can be detected when they are within one million bases apart, the strength of correlation between CpGs and SNPs is dramatically decreased when their distance is increased 67 . Here, we just included those SNPs that are 50 kb away from the CpG site in the SNP-CpG interaction analysis. SNP genotype data were generated via the Illumina HumanOmni1-Quad v1.0 microarray, as described in our previous GWAS on AD 26 and ND 37 . Genotype data were cleaned with a commonly used procedure, that is, SNPs were excluded if they met any of following criteria: 1) minor allele frequency (MAF) was < 0.05; 2) missing genotyping rate was > 10%; or 3) P value of Hardy-Weinberg disequilibrium test was < 1.0 × 10 −3 . The association of SNP genotypes (in an additive model) and CpG methylation levels was analyzed using PLINK 76 . For both EAs and AAs, linkage disequilibrium (LD) plots of SNPs near CpGs were generated by R package LDheatmap 77 .
Additionally, the function of differentially methylated promoter CpGs was predicted by bioinformatics programs. The online program PROMO 78,79 was applied to predict whether differentially methylated promoter CpGs were located in transcription factor binding sites (TFBS) as defined in the TRANSFAC database 80 . The USCS Genome Browser (http://genome.ucsc.edu) was used to query whether differentially methylated promoter CpGs were located in DNase hypersensitivity sites (DHSs) or in chromosomal regions that were associated with H3K27Ac mark (which are often found near active regulatory elements).