Abstract
Alterations in DNA methylation have been associated with alcohol exposure and proposed to contribute to continued alcohol use; however, the molecular mechanisms involved remain obscure. We investigated the escalating effects of alcohol use on DNA methylation, gene expression and predicted neural effects in the nucleus accumbens of rhesus macaques that self-administered 4% alcohol for over 12 months. Using an exploratory approach to identify CpG-rich regions, followed by bisulfite sequencing, the methylation levels of 2.7 million CpGs were compared between seven low-binge drinkers and nine heavy–very heavy drinking subjects. We identified 17 significant differential methylation regions (DMRs), including 14 with methylation levels that were correlated with average daily alcohol consumption. The size of the DMRs ranged from 29 to 158 bp (mean=63.7), included 4–19 CpGs per DMR (mean=8.06) and spanned a range of average methylation values from 5 to 34%. Eight of the DMRs mapped to genes implicated in modulating synaptic plasticity. Six of the synaptic genes have not previously been linked to alcohol use. Validation studies of these eight DMRs using bisulfite amplicon sequencing and an expanded set of 30 subjects confirmed the significant alcohol-dose-associated methylation of the DMRs. Expression analysis of three of the DMR-associated genes, LRP5, GPR39 and JAKMIP1, revealed significant correlations between DMR methylation and whole-gene or alternative transcript expression, supporting a functional role in regulating gene expression. Together, these studies suggest that alcohol-associated synaptic remodeling may be regulated and coordinated at the level of DNA methylation.
Introduction
Chronic and excessive alcohol use can lead to alcohol dependence, a relapsing and remitting condition that ultimately costs lives and disrupts families. In an effort to understand the neuroadaptive changes associated with dependent or compulsive drinking, investigations have focused on brain regions that process motivated behaviors and on cellular mechanisms underlying learning and memory. Accumulating evidence implicates the nucleus accumbens (NAc) core (NAcc) in the control of motivated behaviors by discrete cues.1 Thus, the NAcc can be viewed as a relay station selecting and integrating the most relevant environmental stimuli among competing limbic and cortical afferents to drive behavioral output,2 such as alcohol seeking. In the NAc, chronic alcohol use has been linked to changes in dendritic structures3 and neurotransmitter signaling4 thought to contribute to alcohol tolerance, craving and withdrawal.5 Thus, elucidating the molecular mechanisms that link alcohol use and these neural adaptations remains a challenge for fully understanding and treating alcohol dependence.
DNA methylation is an epigenetic mechanism that mediates the effects of the environment into altered chromatin structure, gene regulation and expression.6 Modified DNA methylation at individual loci has been linked to alcohol dependence,7, 8 and global CpG methylation has been reported to be higher in alcoholic populations.9, 10, 11 A recent genome-wide study identified differential DNA methylation in the prefrontal cortex from human post-mortem alcoholic and non-alcoholic subjects.11 A subset of the differentially methylated CpGs mapped to differentially expressed genes, suggesting that DNA methylation may be contributing to transcriptional regulation in alcoholics, although further studies are needed to verify this relationship. Nonetheless, the dependence of self-reported alcohol use in human studies severely limits the evaluation of alcohol dose effects on neural DNA methylation or associated changes in gene expression. In addition, unknown nicotine or drug use, or comorbid psychiatric and medical conditions have the potential to contribute confounding effects.
The nonhuman primate alcohol self-administration model provides an outstanding opportunity to discover both epigenetic and expression effects associated with alcohol use, overcoming the logistical constraints of human studies. In this model, macaques have 22 h per day access to 4% alcohol and water for a period of 12 months, enabling the collection of precise information of all alcohol consumed.12 In addition, macaques exhibit a natural broad distribution of drinking preferences similar to that of humans, enabling the study of varied alcohol use histories without relying on self-reported data. Furthermore, macaques born and raised in captivity typically have complete medical histories, and the immediate collection of tissues at the time of death provides optimal resources for epigenomic and transcriptomic studies. Finally, owing to the high conservation of genomic sequences, and similar neuroanatomy, the macaque offers high translational relevance for the study of alcohol-associated neuroadaptations.
In the present study, we leveraged the features of the macaque alcohol self-administration model to identify long-term alcohol dose effects on DNA methylation and associated gene expression in the NAcc. We combined the targeted selection of CpG-rich genomic regions, bisulfite sequencing and a statistical clustering approach to identify significant differential methylated regions (DMRs) among subjects that differed in their categorical alcohol consumption levels (low-binge (L/BD) and heavy–very heavy drinkers (H/VHD)). Both DNA methylation validation studies and associated gene expression analysis underscore the significant correlation between average daily alcohol dose, DNA methylation and gene expression. Eight of the genes identified map to synaptic genes, including LRP5, GPR39 and JAKMIP1, encoding proteins modulating the balance between excitatory and inhibitory signaling.
Materials and methods
Subjects
Male rhesus macaques (n=30, Macaca mulatta), that were late adolescents, young adults and middle aged adults (4.3–4.9, 5.4–6.6 and 7.1–10.2 years at the start of the 12-month open access period, respectively; n=8, n=11 and n=11 per age group), were included in this study. All of the monkeys were born and reared at the Oregon National Primate Research Center (ONPRC) with their mothers until 2–3 years of age, and they were initially selected to minimize relatedness; the average kinship coefficient of all subjects was 0.003. Monkeys were individually housed, and all subjects underwent the same experimental conditions. Briefly, monkeys were allowed for visual, auditory and olfactory sensory contact with each other in a colony room with 12:12-h light–dark cycle with lights on at 0700 hours. All of the animal procedures used in this study were approved by the ONPRC IACUC and were performed in accordance with the NIH and the National Resource Council’s Guide for the Care and Use of Laboratory Animals.
Drinking procedure
Voluntary and long-term ethanol self-administration was induced using schedule-induced polydipsia as previously described.13 Briefly, the monkeys were trained daily to use the operant panel and induced to drink 0.0, 0.5, 1.0 and 1.5 g kg−1 ethanol (4%) in 30 day epochs. During the following 12 months, subjects had open access (22 h per day) to 4% alcohol and water ad libitum. The alcohol intake data were collected and recorded in an automated manner by computer. Accumulative data on the alcohol self-administration consumption patterns in these macaques during the 12 months of open access have identified four stable and distinguishable categories, named low, binge, heavy and very heavy drinking. A combination of blood ethanol concentration (BEC), average g kg−1 per day consumption and percentage of days over a certain threshold have been identified as accurate parameters in distinguishing these four stable drinking patterns.12 Accordingly, subjects consuming >2 g kg−1 for more than 55% of the days were defined as binge drinkers (BDs). Those subjects consuming >3 g kg−1 for 20% of the days were classified as heavy drinkers (HDs), whereas the VHDs had more than 10% of the days with >4 g kg−1. Lastly, low drinkers (LDs) were those subjects that spent less than 55% consuming more than 2 g kg−1. Importantly, LD drinkers and BDs never or occasionally (respectively) reach BECs above 80 mg dl−1, the baseline measure of human intoxication.12 In contrast, HDs and VHDs routinely measure BECs above 80 mg dl−1. In the present study, the four categories were further combined into two groups based on propensitiy of the subjects to be intoxicated (BEC>80 mg dl−1). Thus, LDs and BDs were combined into L/BDs, whereas HDs and VHDs into H/VHDs.
This study included 17 L/BDs and 13 H/VHDs. Sixteen subjects (seven L/BDs and nine H/VHDs) were included in the genome-wide bisulfite sequencing study; the sample size was expanded to 30 subjects for the bisulfite amplicon sequencing (BSAS) validation study. On the basis of sample availability, 19 of the same animals were included in the RNA analysis (10 L/BDs and 9 H/VHDs).
Tissue collection and genomic DNA isolation
After the 12-month open access period, and without imposed abstinence, a previously described, detailed necropsy protocol14 was used to systematically collect tissues from all subjects. Briefly, monkeys were anesthetized with ketamine (10 mg kg−1), maintained on isoflurane and perfused with ice-cold oxygenated monkey perfusion solution (containing (in mM) 124 NaCl, 23 NaHCO3, 3 NaH2PO4, 5 KCl, 2 MgSO4, 10 D-glucose, 2 CaCl2). Brains were quickly removed and sectioned along the coronal plane using a brain matrix.15The block containing the NAcc was initially selected by each individual’s magnetic resonance imaging and verified using visible landmarks. In macaques, the NAcc is ~2 mm × 2 mm and extends ~3 mm rostral/caudal.16 The core is differentiated from the shell based on visible landmarks. Using the curvature of the internal capsule, the area just ventral to its end is the NAcc. From the frozen 4 mm coronal brain block maintained on dry ice, a small circular dissection of ~1 mm3 was made, taking care to not collect white matter from the tract (dorsal to the core). This relatively small dissection avoids the NAc shell and yields enough tissue for nucleic acid isolation. Genomic DNA and RNA were extracted from the NAcc using the All Prep DNA/RNA/microRNA Universal Kit (Qiagen Sciences, Germantown, MD, USA) following the manufacturer’s recommendations.
High-throughput DNA methylation analysis
Three micrograms of genomic DNA were sheared using a Bioruptor UCD200 (Diagenode, Denville, NJ, USA), generating fragments ~180 bp. The SureSelect XT Human Methyl-Seq library preparation (Agilent Technologies, Santa Clara, CA, USA) was used following the manufacturer’s instructions. The libraries were then bisulfite-treated using EZ DNA Methylation-Gold (Zymo Research, Irvine, CA, USA), and quantified using a 2100 Bioanalyzer (Agilent Technologies). DNA libraries were sequenced on an Illumina HiSeq2500 at the OHSU Massively Parallel Shared Sequence Resource.
CpG methylation rate analysis
The quality of the bisulfite-converted sequencing reads was assessed with FastQC. Reads were trimmed and aligned to the rhesus macaque reference genome (MacaM17), and then the bisulfite conversion rates were evaluated, insuring all libraries were >98% converted, and CpG methylation was evaluated using Bismark.18 The methylation rates were calculated as the ratio of methylated reads over the total number of reads. Methylation rates for CpGs with fewer than 10 reads were excluded from further analysis. The remaining CpGs (2.7 million) had an average of 60 × read coverage. All sequence reads were submitted to the Sequence Read Archive at NCBI under project accession number PRJNA294610. An overview of these results is described in Supplementary Table 1. The differential analysis of the CpG methylation levels is described in Statistical analysis below.
Bisulfite amplicon sequencing
Candidate DMR methylation levels were validated using targeted BSAS. Primers were designed within 200 bp of each DMR, using the Bisulfite Primer Seeker tool from Zymo Research (Supplementary Table 2). Each gDNA (250 ng) was bisulfite-converted using EZ DNA Methylation-Gold (Zymo Research) and 12.5 ng of bisulfite-converted DNA was used for each polymerase chain reaction (PCR). Library construction, analysis of the percent methylation at each CpG in each amplicon and PCR allele bias correction were performed as previously described.19
High-throughput real-time PCR
RNA extracted from the same NAcc tissues was used for quantitative reverse transcriptase-PCR (RT-PCR) analysis. The Fluidigm Reverse Transcription Master Mix (Fluidigm, San Francisco, CA, USA) was used to reverse-transcribe 100 ng of each RNA sample following the manufacturer’s instructions. Next, the complementary DNA was pre-amplified, and unincorporated primers were removed following the manufacturer’s instructions. The reactions were diluted (10 × ) with 43 μl of TE buffer (TEKnova, Hollister, CA, USA).
qPCR was performed using the BioMark HD System and the 96.96 GE Dynamic Arrays (Fluidigm) in triplicate assays. The Fluidigm sample premix and the assay premix were prepared following the manufacturer’s instructions. The samples and reagents were mixed using the Nanoflex IFC controller (Fluidigm). Thermal qPCR conditions were as follows: 95 °C for 60 s, 35 cycles of 95 °C for 5 s and 60 °C for 20 s. Data were processed by automatic threshold for each assay, with derivative baseline correction using the BioMark Real-Time PCR Analysis Software 3.1.2 (Fluidigm). The quality threshold was set at the default of 0.65.
The primer sequences are listed in Supplementary Table 3.
The mRNA expression levels were normalized as previously described,19 except for using the geometric mean of three constitutive genes: B-Actin, Tubulin1 and Phosphoglycerate kinase (PGK1 (ref. 20)). We also confirmed that different levels of alcohol use did not affect their expression in the NAcc (data not shown).
Statistical analysis
The exploratory nature of the present study limited our ability to estimate a priori the sample size needed to detect the effects of a broad range of alcohol on the methylome. Thus, we used a group-comparison (L/BD versus H/VHD) that was previously established, based on the specific drinking parameters during 12 months of alcohol consumption, to maximize the opportunity to detect significant effects of alcohol dose on DNA methylation levels.
One of the limitations of the analysis of genome-wide bisulfite sequencing data is that the significance of the individual CpG sites may be dampened after multiple-testing correction on potentially millions of sites. Thus, we employed the comb-p method that combines P-values in sliding windows and accounts for spatial correlations across the genome.21 In detail, the CpG methylation rates of L/BD (n=7) and H/VHD (n=9) subjects were first subjected to the Wilcoxon two-sample Independent test using the wilcox.test function in the R programming language.22 The comb-p method uses a sliding window correction where each Wilcoxon P-value is adjusted by applying the Stouffer–Liptak–Kechris (slk) method23, 24, 25 of neighboring P-values as weighted according to the observed autocorrelation (ACF) at the appropriate lag. Briefly, comb-p first calculates the ACF at varying distance lags, and then the ACF is used to perform the slk correction where each P-value is adjusted according to adjacent P-values as weighted according to the ACF. Thus, a given P-value will be pulled lower if its neighbors also have low P-values and likely remain insignificant if the neighboring P-values are also high. Next, a q-value score based on the Benjamini–Hochberg false discovery rate correction is calculated. The peak-finding algorithm is used to find enrichment regions. Once the regions are identified, a P-value for each region can be assigned using the Stouffer–Liptak correction. Then, the false discovery rate q-value is used to define the extent of the region, whereas the slk-corrected P-value and a one-step Sidak multiple-testing correction26 is used to define the significance of the region. Parameters for Comb-p were DIST=300, STEP=60 and THRESHOLD=0.05.
The Shapiro–Wilk test (appropriate for small sample sizes) was used to assess the normality of the average methylation and single CpG methylation rates for the eight amplicons analyzed by BSAS, and also the mRNA expression for LRP5, GPR39 and JAKMIP1. The variables analyzed followed a normal distribution. We used the 1.5xIQ (interquartile range) method to identify outliers, which were subsequently excluded from the corresponding analysis. The sample size for each analysis is specified in the Results section as well as in the figure legends.
Before applying Independent t-test to compare the BSAS average methylation rates and mRNA expression between L/BDs and H/VHDs, the Levene’s test was used to test homogeneous variance assumption for the t-test (pARHGEF7-DMR=0.099, pCDH5-DMR=0.028, pJAKMIP1-DMR=0.050, pKIRREL3-DMR=0.676, pGPR39-DMR=0.177, pNTM-DMR=0.064, pLRP5-DMR=0.309, pNBEA-DMR=0.302, pJAKMIP1-exon 1B-mRNA=0.021, pJAKMIP1-exon 1A-mRNA=0.921, pGPR39-mRNA=0.145, pLRP5-mRNA=0.207). When the variance was heterogeneous (CDH5-DMR, JAKMIP1-exon 1B-mRNA), the Welch–Satterthwaite method was used for estimating the s.e.
Pearson correlation analysis was used to explore associations between mRNA expression, DMR average methylation and average ethanol (g kg−1 per day) consumption.
All statistical analyses were carried out using IBM SPSS Statistics (Armonk, NY, USA), R and comb-p software, with values α<0.05.
Results
Description of the subjects and their ethanol-drinking patterns
Self-reported data indicate that human subjects with alcohol use disorders consume a wide range of alcohol, from 0.7 to >4 g kg−1 per day.11, 27, 28, 29 Similarly, the macaques in this study consumed between an average of 0.5 and 4.4 g kg−1 per day during the 12 months of open access to ethanol and water (Supplementary Figure 1b). The subjects were classified as H/VHD if they had more than 20% of days with 3 g kg−1 or more ethanol consumption, and L/BD if they did not exceed the 20% threshold.12 The mean daily alcohol use was 1.9 and 3.2 g kg−1 per day for L/BD and H/VHD, respectively. In addition, we note that the L/BD group very rarely or occasionally reached BECs of 80 mg dl−1, a measure of alcohol intoxication, whereas the H/VHD members regularly exceeded the 80 mg dl−1 threshold. We previously confirmed that the amount of ethanol consumed did not reflect general differences in thirst and were not associated with the age of the subjects.19
Alcohol-associated differentially methylated CpGs in the NAcc
Using an enrichment approach targeting GENCODE promoters,30 CpG islands (CGI), shores and shelves, and DNase I-hypersensitive sites located in or near RefGenes, followed by bisulfite sequencing, we analyzed roughly 3.1 million CpGs per subject. Of these, 2.7 million CpGs per individual met our quality requirements (Supplementary Table 1). In all, 74 032 CpGs (2.7%) had significant methylation differences between L/BD and H/VHD subjects (Wilcoxon test). The significant CpGs showed a greater proportion of hypermethylated CpGs in H/VHDs (51 500/74 032=70%; Supzplementary Figure 2).
Alcohol-associated DMRs in the NAcc
Although the methylation of a single CpG can influence gene expression,31 it is challenging to pinpoint functional methylation changes among the thousands of significant CpGs detected. An alternative approach is to identify contiguous differentially methylated CpGs, which minimizes false discovery and enhances the prospect of identifying functional effects.21, 32 In this study, we applied the comb-p method,21 which identifies regionally correlated P-values, applies a false discovery rate correction to define the extent of the region and a one-step Sidak multiple-testing correction26 to define the significance of the DMR. This approach identified a set of 17 DMRs distinguishing L/BDs and H/VHDs (Table 1). The size of the DMRs ranged from 29 to 158 bp (mean=63.7), included 4–19 CpGs per DMR (mean=8.06) and had an average CpG density of 0.13. Neighboring CpGs within each DMR exhibited concordant DNA methylation differences. Whereas alcohol consumption was generally associated with higher DMR methylation (Table 1), 2 of the 17 DMRs were hypomethylated among H/VHDs (DMRs linked to NTM and LRP5; Table 1). The global CpG methylation within these two DMRs was negatively correlated with the daily average amount of alcohol consumed (seven L/BD and eight H/VHD, one subject had no data for these two DMRs; r=−0.681, P=0.007; Figure 1b). Alternatively, the global average CpG methylation within the remaining 15 DMRs (seven L/BD and nine H/VHD) was positively associated with the daily average amount of ethanol consumed by the 16 subjects included in this study (r=0.752, P=0.0004; Figure 1a). Overall, 14 of the DMRs had average CpG methylation levels that were significantly correlated with average g kg−1 per day ethanol (Table 1), suggesting that DNA methylation of these regions is dynamically modified by the amount of alcohol consumed.