G C A T genes T A C G G C A T Article DNA 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 AlexeyKozlenkov1,2,AndrewE.Jaffe3,4,5,AlisaTimashpolsky1,PashaApontes1, SergeiRudchenko6,MihaelaBarbu6,WilliamByne1,2,YasminL.Hurd2,SteveHorvath7,8and StellaDracheva1,2,* 1 JamesJ.PetersVAMedicalCenter,Bronx,NY10468,USA;[email protected](A.K.); [email protected](A.T.);[email protected](P.A.);[email protected](W.B.) 2 TheFriedmanBrainInstituteandDepartmentofPsychiatry,IcahnSchoolofMedicineatMountSinai, NewYork,NY10029,USA;[email protected] 3 LieberInstituteforBrainDevelopment,JohnsHopkinsMedicalCampus,Baltimore,MD21205,USA; [email protected] 4 DepartmentofBiostatisticsandDepartmentofMentalHealth,JohnsHopkinsBloombergSchoolofPublic Health,Baltimore,MD21205,USA 5 CenterforComputationalBiology,JohnsHopkinsUniversity,Baltimore,MD,21205,USA 6 HospitalforSpecialSurgery,NewYork,NY10021,USA;[email protected](S.R.); [email protected](M.B.) 7 DepartmentofHumanGenetics,DavidGeffenSchoolofMedicine,UniversityofCaliforniaLosAngeles, LosAngeles,CA90095,USA;[email protected] 8 DepartmentofBiostatistics,SchoolofPublicHealth,UniversityofCaliforniaLosAngeles,LosAngeles, CA90095,USA * Correspondence:[email protected];Tel.:+1-718-584-9000(ext.6085) AcademicEditor:DennisR.Grayson Received:14February2017;Accepted:25May2017;Published:30May2017 Abstract: WeemployedIllumina450KInfiniummicroarraystoprofileDNAmethylation(DNAm) inneuronalnucleiseparatedbyfluorescence-activatedsortingfromthepostmortemorbitofrontal cortex(OFC)ofheroinuserswhodiedfromheroinoverdose(N=37),suicidecompleters(N=22) with no evidence of heroin use and from control subjects who did not abuse illicit drugs and died of non-suicide causes (N = 28). We identified 1298 differentially methylated CpG sites (DMSs)betweenheroinusersandcontrols,and454DMSsbetweensuicidecompletersandcontrols (p<0.001). DMSsandcorrespondinggenes(DMGs)inheroinusersshowedsignificantdifferences inthepreferentialcontextofhyperandhypoDM.HyperDMSswereenrichedingenebodiesand exonsbutdepletedinpromoters,whereashypoDMSswereenrichedinpromotersandenhancers. Inaddition,hyperDMGsshowedpreferenceforgenesexpressedspecificallybyglutamatergicas opposedtoGABAergicneuronsandenrichmentforaxonogenesis-andsynaptic-relatedgeneontology categories,whereashypoDMGswereenrichedfortranscriptionfactoractivity-andgeneexpression regulation-relatedterms. Finally,wefoundthattheDNAm-based“epigeneticage”ofneuronsfrom heroinuserswasyoungerthanthatincontrols. Suicide-relatedresultsweremoredifficulttointerpret. Collectively,thesefindingssuggestthattheobservedDNAmdifferencescouldrepresentfunctionally significantmarksofheroin-associatedplasticityintheOFC. Keywords: DNAmethylation;drugaddiction;heroin;suicide;brain;neurons;human Genes2017,8,152;doi:10.3390/genes8060152 www.mdpi.com/journal/genes Genes2017,8,152 2of18 1. Introduction OpioidoverdoseisnowthesecondleadingcauseofaccidentaldeathamongadultsintheU.S.[1,2]. Prescription opioids are believed to have served as a gateway to the use of heroin which is more readilyavailableandlessexpensivethanprescriptionmedications[3–5]. Theopioidepidemichas raisedattentiontotherelativelylimitedknowledgeaboutthepathophysiologyunderlyingheroinuse disorder,particularlyfrominsightsgainedthroughinterrogationofthehumanbrain. Aswithotherdrugsofabuse,thesusceptibilitytoopioidaddictionisknowntobeinfluenced roughlyequallybygeneticandenvironmentalfactors[6–8],suggestinganimportantroleforepigenetic regulation. Epigeneticmechanismsmediatelong-termchangesingeneexpressionwithoutchanges in DNA sequence [9,10]. Recent studies in animal models provide robust evidence that repeated exposure to drugs of abuse induces changes in gene expression through alterations in all major modes of epigenetic regulation (e.g., cytosine DNA methylation (DNAm), histone modifications, and non-coding RNAs), and in several instances, the contribution of such epigenetic changes to addiction-relatedbehavioralabnormalitiesinanimalshasbeendirectlydemonstrated[11]. Specific examplesincludechangesinthelevelofthetranscription-activatinghistonemodificationH3K27ac indifferentbrainregionsinseveralexperimentalmodelsofdrugaddiction[12,13],alterationsofthe repressivehistonemarksH3K9me2/3inthenucleusaccumbensduringchroniccocaineoropioid addiction[14,15],andaplethoraofepigeneticmechanismswhichhaverecentlybeenimplicatedinthe downregulationofBDNFexpressionintheventraltegmentalareaafterchronicopioidexposure[16]. DNAmisastablemostlyrepressiveepigeneticmodificationthatisextremelyimportantforboth theestablishmentofcell-type-specificphenotypesinthenervoussystem[17]andthemediationof environmentallyinducedchangesintheadultbrainincludingmemoryformation,stressresponses, anddepression[18–20]. ChangesinlevelsofDNAmatspecificgeneloci,aswellaschangesinlevels ofmodifyingproteins,havealsobeenobservedafterexposuretodrugsofabuse,witharecentstudy supportingtheroleofDNAhydroxymethylation(5hmC)inchroniccocaineaddiction[21,22].Whereas suchinformationfromanimalstudiesisvaluable,increasingthelimitedknowledgeabouttheDNAm landscapeinthebrainsofhumanswithsubstanceusedisordersiscriticalforthedevelopmentofnovel moreeffectivetreatments. SeveralfeaturesofDNAmwereconsideredinthedesignofthepresentstudy. First,recentstudies indicatethatthecompositionanddynamicsofDNAmarenotonlydistinctinthebraincomparedto othertissues[23,24],butalsodiffersignificantlyacrossdifferentregionsofthebrain[25,26].Studies ofcognitiveprocessesthataccompanydrugaddictionduringthelastdecadesuggestthataddiction involves neuroplasticity mechanisms similar to traditional models of learning and memory [27] and that these mechanisms underlie the role of the medial (m) prefrontal cortex (PFC) in drug self-administration and the long-lasting propensity to relapse [28–30]. Moreover, human imaging studies [31–33] and rodent studies that employed a reinstatement model of drug relapse [34–36] stronglysuggestinvolvementoftheventralaspectsofthemPFCinaddiction. Inparticular, these studiessuggestthatdiminishedoutputfromtheventralmPFCcontributestodrugseekingbehavior byimpairingtheabilitytoactivelyinhibitbehavioralresponsestodrug-conditionedstimuli. Forthese reasons,inthepresentstudyweexaminedDNAminautopsyspecimensfromthemedialorbitalfrontal cortex(mOFC)—a ventralsubregion ofthemPFC—which participatesinregulatinggoal-directed behavioranddecision-makingandhasbeenimplicatedindrugaddictionbyhumanimagingand animalstudies[37,38]. In addition, a number of studies have reported DNAm alterations in the human brain that are associated with addiction (specifically with alcohol use disorder, e.g., [39,40]). These studies were performed using bulk (cellularly heterogeneous) tissues. However, several recent reports haveclearlydemonstratedrobustdifferencesinDNAmandhistonemodificationpatternsbetween neuronalandglialcellpopulationsinhumanandrodentbrains[24,41,42]. Thesecelltype-specific epigeneticlandscapesmightultimatelydeterminetheselectivevulnerabilitytoneurodevelopmentalor environmentalinsultsthatcouldculminateindrugaddiction. Thus,inordertoincreasethelikelihood Genes2017,8,152 3of18 ofidentifyingcelltype-specificsignaturesofheroinaddiction-associatedepigeneticvariations,our studiesemployedneuronalcellsthatwereseparatedfromgliausingfluorescenceactivatednuclei sorting(FANS). Recently,ahighlyaccuratemulti-tissueepigeneticbiomarkeroftissueage(“epigeneticclock”, also known as “Horvath clock”) based on DNAm levels has been introduced [43]. This approach allows one to estimate the “epigenetic age” of any organ, tissue, or cell type including sorted neurons. Mathematically it is defined as a weighted average across 353 CpG sites. The resulting ageestimate(inunitsofyears)isreferredtoas“DNAmethylationage”(DNAmage)or“epigenetic age”. Recent studies support the idea that epigenetic age estimates can serve as biomarkers of biologicalage. Forinstance,theepigeneticageofbloodhasbeenfoundtobepredictiveofall-cause mortality[44–48],frailty[49],lungcancer[50],andcognitiveandphysicalfunctioning[51]. Further, theutilityoftheepigeneticclockmethodusingvarioustissuesandorganshasbeendemonstrated inapplicationssurroundingAlzheimer’sdisease[52],centenarianstatus[47,53],development[54], Downsyndrome[55],HIVinfection[56],Huntington’sdisease[57],obesity[58],lifetimestress[59], andParkinson’sdisease[60]. ItwasalsoshownthatdifferentbrainregionshavedifferentDNAmage, withthecerebellumdisplayingtheyoungestageofalltestedregions[53]. In the present report, we compared DNAm in neuronal populations of the mOFC among individualswhoabusedheroinanddiedofheroinoverdose,suicidecompleterswithoutanyevidence ofheroinuse,andacontrolgroupconsistingofspecimensfrompsychiatricallynormalindividuals whodiedofnon-suicidecausesanddidnotabuseillicitdrugs. Allbrainspecimenscamefromthe same brain collection. Suicide completers were included as a comparison group with a different pathophysiology. The study revealed significant methylation disturbances in heroin abusers that mayrepresentafunctionallyrelevantepigeneticsignatureofheroinaddictioninthehumanbrain. Moreover,whiletheepigeneticageofneuronsfromthesuicideandcontrolindividualsdidnotdiffer fromtheirbiologicalage,theepigeneticageofneuronsfromheroinuserswasyounger. 2. MaterialsandMethods 2.1. SubjectsandTissue Human postmortem brain specimens were obtained from our Brain Collection at the Icahn SchoolofMedicineatMountSinaithathasbeenextensivelyusedformanymolecularandepigenetic studies [42,61–69]. Brains were collected at autopsy at the Department of Forensic Medicine, SemmelweisUniversity(Hungary)orNationalInstituteofForensicMedicine,KarolinskaInstitutet (Sweden). Allmaterialwasobtainedunderapprovedlocalethicalguidelines. Thecohort(N =88) consistsofaHungarian/Swedishpopulationofheroinabuserswhodiedofheroinoverdose(N=37), suicidecompleterswithoutanyevidenceofheroinuse(N =22), andcontrolsubjectswhodidnot abuseillicitdrugsanddiedofnon-suicidecauses(suchascardiacfailure,viralinfection,oranaccident; N=29)(SupplementaryMaterialsTableS1). Oneofthecontrolsubjectswaslaterexcludedfromthe analysesduetoambiguityofsex,whichwasdiscoveredduringthedataprocessing. Thecauseand mannerofdeathandpossiblepsychiatricdiagnosesweredeterminedbyaforensicpathologistafter evaluatingautopsyresults,circumstancesofdeath,datafromextensivetoxicologicaltesting,police reports,familyinterviews,andmedicalrecords. Exclusioncriteriawerepostmorteminterval(PMI) of>24h,HIV-positivestatus,historyofalcoholism,useofillicitdrugs(otherthanheroininheroin abusers)andthepresenceofAxis1psychiatricdisorders. Allspecimensfromheroinusersincluded inthestudywerefromindividualswhodiedfromheroinintoxication. Theseindividualswerenot receivingmethadoneorbuprenorphinetreatment, andhadpositivebloodand/orurinelevelsfor opiatesatthetimeofdeath. Theiraveragetimeofheroinusewherethisinformationwasavailable (N =23individuals)was3.75years(rangingfrom0.5to10years). Thecontrolandsuicidegroups werenegativeforbloodopiatesandhadnohistoryofanydrugaddiction. Allsuicidesubjectsdied fromhanging. Thecontrolsubjectsdiedfrommyocardialinfarction,pulmonaryembolism,electric Genes2017,8,152 4of18 shockorviralinfection. Subjectsfromall3groupsshowednegativetoxicologyforothercommon drugsofabuseandcommontherapeuticagents. EtOHwasdetectedinthebloodand/orurineof 6 heroin, 5 suicide, and 3 control subjects and EtOH levels were not significantly different among groups. Nicotinetoxicologywasnotconducted,buttobaccouseisfrequentinthegeneralpopulation fromwhichallthesubjectswerecollected. Autopsybrainspecimenswerecutcoronallyin1cmslabs,frozen,andkeptat−80◦C.Fromthese specimensweharvestedtheventralextentofthePFCcommonlyreferredtoastheorbitalfrontalcortex (OFC;Brodmannarea11). Specifically,wedissectedtheregionscontainingthemedialorbital(MOrG) andinferiororbital(IOrG)gyrijustanteriortothetransverselyrunningorbitalsulcus. Theregionwas dissectedinasingleblockboundedbytheolfactorysulcusmediallyandtheinferiororbitalsulcus laterallyasdescribedpreviously[42]. 2.2. NucleiSeparationbyFluorescenceActivatedNucleiSorting(FANS) Cellstructureisnotpreservedinfrozenautopsybrainspecimens. However,thenucleiofdifferent cell types remain intact. Antibodies against the RNA-binding protein NeuN, which is expressed exclusively in the neuronal nuclei, have been used to separate neuronal from glial nuclei using FANS (e.g., [24,41,70–73]). In a recent study, we optimized previously published FANS protocols byemployingtheDNA-bindingdye7-AADandanti-NeuNantibodiesdirectlyconjugatedwiththe fluorophore[42]. Thisprotocolwasusedinthepresentstudy. Inshort, foreachspecimen, mOFC tissuewasgroundusingmortarandpestleonliquidnitrogen,resuspendedinice-coldLysisBuffer (0.1% Triton, 0.32 M sucrose, 5 mM CaCl , 3 mM MgCl , 10 mM Tris-HCl), filtered through a cell 2 2 strainer,andcentrifugedfor5minat300g. ThepelletwasresuspendedinBlockingBuffer(1%goat serum,2mMMgCl ,TBS)andincubatedfor45minwithAlexa488-conjugatedanti-NeuNantibodies 2 (EMDMillipore,Billerica,MA,USA)(1:1000dilution). Next,asecondcentrifugationstep(15min, 2800g)throughalayerof1.1Msucrosewasdone,andtheresultingpelletwasresuspendedinPBS. TheDNAdye7-AAD(Sigma-Aldrich, St. Louis, MO,USA)wasaddedtoafinalconcentrationof 2µg/mL,andthesamplewassubjectedtotheFANSprocedureusingVantagewithDiVa(excitation wavelength 488 nm). Finally, the sorted nuclear fractions were precipitated by centrifugation at 4000rpmfor20minat4◦Candstoredfrozenat−80◦CuntilDNAisolation. Thelatterwasperformed usingproteinaseKtreatmentandtworoundsofphenolextractionfollowedbyethanolprecipitation. This protocol allowed us to routinely obtain well-separated NeuN(+) and NeuN(−) nuclear fractions, with the width of separation reaching up to an order of magnitude of the NeuN signal intensity[42]. Inaddition,theDNAcontentofbothfractionswaswell-definedbecauseaggregates, nuclei of dividing cells, and debris were excluded in the process of sorting. Validation of the cell-typespecificityoftheobtainedneuronalandnon-neuronalpopulationswasdonepreviouslyby demonstratingtheenrichmentfortheknownneuronalorglial-specifictranscriptsinRNAsamples extracted from the sorted NeuN(+) and NeuN(−) populations, respectively [42]. We estimated the proportions of neurons and glia in our FANS-separated neuronal nuclear preparations using the algorithm from [74] and NeuN(+) and NeuN(−) reference data from [41] (see [75] for details). For comparison, we also included DNA methylation data for 6 NeuN(−) (glial) specimens from our published study [42]. The results are presented in a Supplementary Materials Figure S1, and demonstratehighpurityofourneuronalnuclearpreparations. 2.3. DNAMethylationMeasurementandAnalysis For each specimen, DNA was extracted from the neuronal fraction and subjected to sodium bisulfitetreatmenttogeneratemethylation-specificbasechangesbeforehybridization. Batcheffect wasminimizedbyperformingthebisulfitetreatmentsimultaneouslyforall88specimens, andby randomizingtheplacementofheroin,suicideandcontrolsamplesacrossthearrays. DNAsampleswerebisulfiteconvertedusingEZDNAmethylationkit(ZymoResearch,Tustin, CA,USA).Specifically,500ngofhighqualitygenomicDNA(A260/260≥1.8;A260/230=2.0–2.2) Genes2017,8,152 5of18 wasdenaturedbyincubationwithNaOH-containingZymoM-Dilutionbufferfor15minat37◦C. Next,thedenaturedDNAwasincubatedwithbisulfite-containingCT-conversionreagentfor16h at50◦Cinathermocycler. Every60minthereactionwasheatedto95◦Cfor30s. All88samples wereprocessedonthesameplate. TheInfiniummethylationassaywascarriedoutasdescribed[76]. In short, 4 µL of bisulfite-converted DNA (~150 ng) was used in the whole-genome amplification reaction. Afteramplification,theDNAwasfragmentedenzymatically,precipitatedandre-suspended inhybridizationbuffer. AllsubsequentstepswereperformedfollowingthestandardInfiniumprotocol (UserGuidepart#15019519A).ThefragmentedDNAwasdispensedontotheHM450KBead-Chip[77], whichwasfollowedbyhybridizationinahybridizationovenfor20h. Afterhybridization,thearray wasprocessedthroughaprimerextensionandanimmunohistochemistrystainingprotocoltoallow detectionofasingle-baseextensionreaction. Finally,theBead-Chipwascoatedandthenimagedonan IlluminaiScan(Illumina,SanDiego,CA,USA).ThesameiScanarrayscannerwasusedforprocessing allsamples. All88samplesthatwereusedinthisstudypassedIlluminaqualitycontrolrequirementsand receivedastatusof“SuccessfulSample”. Thedatawerefurtherevaluatedforqualityusingthe“minfi” R/Bioconductor package [78]. All samples passed by the suggested criteria of median M and U intensitiesgreaterthan10.5. Allbutonesamplepassedatestforconcordancebetweenestimatedand reportedsex. Thissamplewasremovedfromthesubsequentanalyses(seeSupplementaryMaterials Table S1). The data were then pre-processed with stratified quantile normalization described in the minfi paper [78]. Probes with annotated/dbSNP-labeled SNPs in the single base extension or target CpG site were filtered, as were probes on the sex chromosomes, leaving 456,513 probes for subsequentanalysis. Fordifferentialmethylationanalysis,weemployedlimmaR/Bioconductorpackagetoperform linearregressionandmoderatedt-statistics(withempiricalBayes)adjustingforage,sexandtissuepH andthefirstfour“negativecontrol”principalcomponents,whichtypicallycapturebatchandslide effectsasdescribedpreviously[75].Differentiallymethylated(DM)CpGsites(DMSs)wereassignedto genesusingthedistanceToNearestfunctionintheGenomicRangesRpackage(http://bioconductor.org/ packages/release/bioc/html/GenomicRanges.html)andtheKnownGenedatasetfromtheUniversity ofCalifornia,SantaCruiz)(UCSC)(http://genome.ucsc.edu). 2.4. EstimationofDNAMethylationAge DNAmage(alsoreferredtoasepigeneticage)wascalculatedfromtheneuronalsamplesprofiled withtheIlluminaInfinium450Kplatformasdescribedin[43]. Briefly,theepigeneticclockisdefined asamultivariatelinearmodelforpredictingagebasedontheDNAmlevelsof353CpGs. TheseCpGs andtheirweights(coefficientvalues)werechosenusingseveralindependentdatasetsbyregressing chronologicalageonCpGs. Predictedage,referredtoasDNAmage,correlateswithchronological ageinsortedcelltypes(CD4+Tcells,monocytes,Bcells,glialcells,neurons),tissues,andorgans, including: wholeblood,brain,breast,kidney,liver,lung,andsaliva[43]. Inourstudy,theepigenetic clock method was implemented using publicly available R software scripts and in a web-based calculator. Themeasureofepigeneticageaccelerationwasdefinedasarawresidualresultingfrom regressingDNAmageonchronologicalage. Bydefinition,epigeneticageaccelerationdoesnotcorrelat withchronologicalage(r=0). 2.5. GeneOntologyFunctionalAnnotationAnalysis Geneontology(GO)analysiswasperformedwiththeonlinesoftwaretoolWebGestalt(www. webgestalt.org,[79]). EnrichedGOtermswithadjustedp-values<0.01(Benjamini-Hochbergmultiple test adjustment) were considered statistically significant. In addition, we applied the cutoff of N ≥ 10 genes for the minimal number of genes associated with a specific GO term. For the GO analysis,weincludedthegeneswhichhadbothhyperDMSsandhypoDMSsintobothhyperDMgene (G)andhypoDMGlists. Genes2017,8,152 6of18 2.6. GeneExpressionComparisonbetweenGlu-andGABANeuronsintheHumanPFC WeemployedourrecentlydevelopedFANS-basedprotocol[69]toisolatetwotypesofneuronal nuclei: (1) the MGE-derived GABA neurons, and (2) glutamatergic neurons with a small (~10%) admixtureofnon-MGEderivedGABAneurons(denoted“Gluneurons”). PFCtissuesamplesfrom 3controlsubjectswereusedfortheexperiments. Afternuclearisolation,weemployedanoptimized versionofourRNAisolationprotocol(see[69]). Specifically,40,000nucleiweresorteddirectlyinto 150uLofExtractionBufferfromthePicoPureRNAisolationkit(ThermoFisherScientific,Waltham, MAUSA).RNAisolationwasthenperformedaccordingtothekitprotocol,withtheinclusionofthe on-columnDNAsetreatmentstep,andtheRNAwaselutedin15uLElutionBuffer. RNA-seqlibraries werethenpreparedfrom10ngRNAusingtheSMARTerStrandedPicoRNA-seqlibrarypreparation kit(Clontech,MountainView,CA,USA),andsequencedonaHiSeq2500Illuminasequencer(Illumina, SanDiego,CA,USA)usingpaired-end50bpprotocol.FASTQfilesweretrimmedtoremovelowquality readsandadapters(usingScytheandSicklesoftwarepackages(https://github.com/vsbuffalo/scythe, https://github.com/najoshi/sickle)),furthertrimmedby3bpfromthe5(cid:48) endofread1assuggested inthelibrarykitprotocol,andthenmappedtohumanhg19genomewithSTARsoftware[80]. EdgeRR package [81] was used to perform differential expression (DE) analysis between GABA and Glu neurons,withDEcriteriaofFDR<0.01,abs(FoldChange)>2,andacutoffof0.1countspermillion (CPM).SupplementaryMaterialsTableS2presentstheresultinglistsoftheDEgenes. 3. Results 3.1. mOFCNeuronsofHeroinUsersShowYoungerEpigeneticAgeThanNeuronsofNon-AddictedIndividuals Ithasbeensuggestedthatepigeneticageacceleration(thatmeasuresdeviationsbetweenDNAm ageandchronologicalage)capturesaspectsofthebiologicalageofthebraintissue[43]. Weestimated theepigeneticage(alsoknownasDNAmage)ofeachneuronalmOFCsamplebyaveragingtheDNAm levelsof353CpGsprofiledwiththeIlluminaInfinium450Kassay. Asexpected,neuronalDNAm agewashighlycorrelatedwithchronologicalageofsubjectsatthetimeofdeathacrossallsamples (correlationr=0.85,p=1.2×10−25;Figure1a). Wedefinedameasureofepigeneticageacceleration asresidualresultingfromregressingDNAmageonchronologicalage. Thus,apositivevalueofage accelerationindicatesthattheDNAmageofasampleishigherthanexpectedbasedonchronological age, whereas a negative value indicates that the DNAm age is younger than expected based on chronologicalage. Whereastherewerenodifferencesintheneuronalageaccelerationbetweenthe control(N=28)andsuicidesubjects(N=22),theheroinsubjects(N=37)showedayoungerDNAmvs. chronologicalage(p=0.022;Figure1b). Asignificantlyyoungerneuronalepigeneticageinheroin abuserswasobservedwhenwecomparedtheirageaccelerationtothatofallsubjectswhodidnot abuseheroin(controlsandsuicides,N=50;p=0.0082;Figure1c). Thedifferenceinageacceleration wasalsoobservedwhenonlyyoungadults(subjectswithchronologicalage<40)wereconsidered intheanalysis(N=32andN=36fornon-abusersandabusers,respectively;p=0.0082;Figure1d). Becauseinourcohorttheaverageageofheroinindividualsisyoungerthancontrolindividuals,the latteranalysisprovidedadditionalvalidationofourfindingsusinggroupsthatwerebetterbalanced byage. Genes2017,8,152 7of18 Genes 2017, 8, 152 7 of 18 Figure 1. Neurons of heroin users show younger epigenetic age than neurons of non-addicted Figure 1. Neurons of heroin users show younger epigenetic age than neurons of non-addicted individuals. (a) Correlation of neuronal DNAm age with chronological age of all subjects at time of individuals. (a)CorrelationofneuronalDNAmagewithchronologicalageofallsubjectsattime death. Red circles, heroin subjects; black circles, all other subjects; (b) Analysis of epigenetic age ofdeath. Redcircles,heroinsubjects;blackcircles,allothersubjects;(b)Analysisofepigeneticage acceleration in the neurons among heroin users, suicide subjects, and controls. The heroin subjects accelerationintheneuronsamongheroinusers,suicidesubjects,andcontrols. Theheroinsubjects have a younger DNAm vs. chronological age (p = 0.022); (c) Analysis of epigenetic age acceleration in haveayoungerDNAmvs.chronologicalage(p=0.022);(c)Analysisofepigeneticageaccelerationin the neurons comparing heroin users and all subjects who did not abuse heroin. A significantly theneuronscomparingheroinusersandallsubjectswhodidnotabuseheroin.Asignificantlyyounger younger neuronal DNAm age in heroin users was observed (p = 0.0082); (d) Same as in (c), but for neuronalDNAmageinheroinuserswasobserved(p=0.0082);(d)Sameasin(c),butforsubjects subjects with chronological age <40. The titles of the bar plots report the results of a non-parametric withchronologicalage<40. Thetitlesofthebarplotsreporttheresultsofanon-parametricgroup group comparison test (Kruskal Wallis test). comparisontest(KruskalWallistest). 3.2. Differential DNA Methylation Analysis 3.2. DifferentialDNAMethylationAnalysis We compared samples from heroin users, suicide completers and control individuals to identify Wecomparedsamplesfromheroinusers,suicidecompletersandcontrolindividualstoidentify differences in DNAm, adjusting for age, sex and tissue pH as covariates. After correcting for multiple differencesinDNAm,adjustingforage,sexandtissuepHascovariates. Aftercorrectingformultiple comparisons, we did not identify any differentially methylated (DM) CpGs sites between control and comparisons,wedidnotidentifyanydifferentiallymethylated(DM)CpGssitesbetweencontroland suicide specimens using FDR < 0.1, whereas there were 3 sites that were DM between heroin users suicidespecimensusingFDR<0.1,whereastherewere3sitesthatwereDMbetweenheroinusersand and controls. We then used a more liberal threshold of significance (nominal p-value < 0.001). This controls. Wethenusedamoreliberalthresholdofsignificance(nominalp-value<0.001). Thisresulted resulted in 454 DM sites (DMSs) between suicide and control subjects (DMsui sites) and 1311 sites in454DMsites(DMSs)betweensuicideandcontrolsubjects(DMsuisites)and1311sitesthatwere that were DM between heroin users and controls (DMher sites) (Supplementary Materials Table S3), DMbetweenheroinusersandcontrols(DMhersites)(SupplementaryMaterialsTableS3),andweused and we used these DM sites to focus on gene-set-level analyses. Among the DMher sites, there were theseDMsitestofocusongene-set-levelanalyses. AmongtheDMhersites,therewere12DMSsthat 12 DMSs that belonged to the CpH sequence context; these were excluded from the subsequent belongedtotheCpHsequencecontext;thesewereexcludedfromthesubsequentanalysis,retaining analysis, retaining 1298 DMher CpG sites. The absolute differences in average DNAm levels (beta 1298DMherCpGsites. TheabsolutedifferencesinaverageDNAmlevels(betavalues)weregenerally values) were generally small. Nevertheless, the majority of DMSs (1247 out of 1298 for DMher, and small. Nevertheless,themajorityofDMSs(1247outof1298forDMher,and431outof454forDMsui 431 out of 454 for DMsui sites) showed beta value differences >1%, and the smallest beta value sites)showedbetavaluedifferences>1%,andthesmallestbetavaluedifferenceamongtheremaining difference among the remaining DMSs was 0.62%. Because of the relatively small number of DMSs DMSswas0.62%. BecauseoftherelativelysmallnumberofDMSswithlow(<1%)differences,we with low (<1%) differences, we retained all DMSs for the subsequent analyses. The DMher and DMsui retained all DMSs for the subsequent analyses. The DMher and DMsui sites represented largely sites represented largely non-overlapping populations, with only 29 CpG probes showing significant non-overlappingpopulations,withonly29CpGprobesshowingsignificantDNAmchangesinboth DNAm changes in both heroin users and suicide subjects compared with controls (Supplementary heroinusersandsuicidesubjectscomparedwithcontrols(SupplementaryMaterialsTableS4). Materials Table S4). In our cohort, the average age of heroin individuals was younger than control individuals. In our cohort, the average age of heroin individuals was younger than control individuals. In In addition to including age as a covariate in our initial analysis, we also performed a secondary addition to including age as a covariate in our initial analysis, we also performed a secondary analysis analysisinordertobetteraccountfortheroleofageonourresults. Specifically,weleftoutN=10 in order to better account for the role of age on our results. Specifically, we left out N = 10 oldest oldestcontrolindividualstobalancethegroupsbyage,analyzingN=18controlsandN=37heroin control individuals to balance the groups by age, analyzing N = 18 controls and N = 37 heroin individuals. In this new cohort there were no significant difference in age between the heroin and control subjects (p = 0.09). We detected no noticeable effect of age on the regression estimates for the Genes2017,8,152 8of18 individuals. Inthisnewcohorttherewerenosignificantdifferenceinagebetweentheheroinand controlsubjects(p=0.09). Wedetectednonoticeableeffectofageontheregressionestimatesforthe heroinvs. controlcomparison(SupplementaryMaterialsFigureS2),suggestingthatwehadproperly adjustedforageeffectsintheinitialdifferentialmethylationanalysis. Theinformationabouttheyearsofheroinusewasavailableonlyforasubsetofheroinsubjects. Inordertotestifthedurationofheroinabuseinfluencedthedifferentialmethylationdetectedinour study,weleftouttheheroinindividuals(N=14)whichdidnothavetheyears-of-useinformation,and usingthisfilteredcohort,comparedtheresultsofthedifferentialmethylationanalysesemployingtwo differentlinearmodels: (1)heroin/non-heroinstatus,andthecovariatesfromtheinitialanalysis,and (2)Log10(YearsofHeroinUse+1)andthecovariatesfromtheinitialanalysis. Wefoundnosignificant differencesbetweentheresultsfromthesetwomodels(SupplementaryMaterialsFigureS3). Wethendefineddifferentiallymethylatedgenes(DMGs)basedonthepositionoftheDMSs(see MaterialsandMethods). Therewereatotalof1159uniqueDMGsassociatedwithDMhersites(DMher genes)and407genesassociatedwithDMsuisites(DMgenes)(SupplementaryMaterialsTableS5); 87geneswerecommonforbothDMherandDMsuigenelists. AmongtheDMhergenes,660were associatedwithatleastonehyperDMSand531wereassociatedwithatleastonehypoDMS,with 32DMhergenesassociatedwithbothahyper-andahypoDMS.Therewere243hyperDMsuigenes and170hypoDMsuigenes,with6DMsuigenesoverlappingbetweenthesetwolists(Supplementary MaterialsTableS5). TheDMGsthatwerepresentinboththehyperDMGandhypoDMGlistswere excludedfromtheanalysisforneuronsubtypespecificity(seebelow). 3.3. EnrichmentofDifferentiallyMethylatedSitesinGenomicFeaturesandRegulatoryElements Wefirstaskedifheroin-orsuicide-associatedhyper-orhypoDMSswereenriched(ordepleted) withingenicfeatures(i.e.,promoters,exons,introns,entiregenebodies,intergenicregions),putative regulatoryelements(i.e.,predictedenhancers)orCpGisland-relatedfeatures(i.e.,CpGislands,shores, orregionsoutsideofislandsandshores). Wealsotestedifsuchenrichmentscouldbespecifictothe directionoftheDNAmchanges(hyper-vs. hypoDM). We detected clear differences in the patterns of enrichment or depletion between hyper- and hypoDMSswithinmanyoftheassessedfeaturesintheneuronsfromheroinusers(Table1,statistical significancebyFisher’sexacttest). Inparticular,withingenicfeatures,hyperDMhersiteswerestrongly enrichedinexonsandgenebodies,moderatelyenrichedinintrons,butsignificantlydepletedfrom promoters. Incontrast,hypoDMhersiteswereenrichedinpromoters,butdidnotshowsignificant enrichmentordepletioninothergenicfeatures. Likewise, we detected differences between hyper- and hypoDMher sites in their pattern of enrichmentordepletionwithinputativetranscriptionalenhancers. Thegenomiccoordinatesofthese enhancerswerebasedontheChIP-seqpeaksfortheH3K27achistonemarkinthehumanPFCtissue (dataobtainedfromtheNIHRoadmapEpigenomicsMappingConsortium,REMC[82]).Weperformed this analysis separately for the H3K27ac peaks distal from transcription start site (TSS) (>1000 bp distance),whichareindicativeofputativeactiveenhancers,andfortheTSS-proximalH3K27acpeaks, whichareindicativeofactivepromoters. ThehypoDMhersiteswerestronglyenrichedwithindistal peaks(putativeenhancers),whereasthehyperDMhersitesshowedneitherenrichmentnordepletionin putativeenhancers. Inlinewiththeirenrichmentinpromoters,thehypoDMhersiteswereenrichedin proximalH3K27acpeaks,whereasthehyperDMhersiteswerestronglydepletedfrombothpromoters andproximalH3K27acpeaks. Bothhyper-andhypoDMhersitesweredepletedfromthegenomicareasoutsideofCpGislands or shores (the latter is defined as 2000 bp-width-areas at ≤2000 bp from the nearest CpG island). Inaddition,hyperDMSswereslightlyenrichedinCpGislands,whereashypoDMSswereenriched withinshores. Whenthedistributionofthesuicide-associatedDMSswithindifferentgenicfeatureswasanalyzed, onlyfewsignificantresultsweredetected. HyperDMsuisitesweredepletedwithinpromotersaswell Genes2017,8,152 9of18 aswithinproximalH3K27acpeaks,andweremoderatelyenrichedinexons,intronsandgenebodies (Table1). NoenrichmentordepletionwithinanygenicfeaturewasfoundforhypoDMsuisites,which mightbeexplainedbytherelativelysmallnumberofDMswithinthisgroup. Table1.Enrichmentordepletionofhyper-andhypomethylatedsites(DMSs)withingenomicfeatures andputativeenhancers.StatisticalsignificancewascalculatedbyFisher’sexacttest.Promoterswere definedasregionswithin−1000bparoundthetranscriptionstartsite. Significant(p-value<0.01) enrichments or depletions are shown in italic bold. Moderately significant (p-value < 0.05) enrichments/depletionsareshowninitalic.H3K27acChIP-seqdatasetwasobtainedfromtheRoadmap EpigenomicsMappingConsortiumdatabase(REMC,http://www.roadmapepigenomics.org). HeroinDMSs SuicideDMSs GenomicFeatures/ HyperDMSs HypoDMSs HyperDMSs HypoDMSs PutativeEnhancers Odds Odds Odds Odds p-Value p-Value p-Value p-Value Ratio Ratio Ratio Ratio Exons 1.41 1.6×10−4 0.85 0.19 1.36 0.045 0.85 0.54 Introns 1.21 0.013 0.97 0.75 1.31 0.030 1.15 0.38 Genicfeatures GeneBody 1.26 0.0025 1.06 0.52 1.29 0.044 0.94 0.70 Intergenic 1.11 0.20 0.86 0.15 0.95 0.78 1.14 0.43 Promoters 0.47 4.8×10−15 1.32 0.0034 0.53 6.1×10−5 0.83 0.34 H3K27ac DistalH3K27ac 1.18 0.07017 1.52 3.0×10−5 0.94 0.76 1.11 0.57 peaks(enhancers) ChIP-seqpeaks (humanPFC) ProximalH3K27ac 0.33 2.0×10−15 1.52 1.2×10−4 0.67 0.048 0.60 0.052 peaks(promoters) CpG CpGislands 1.18 0.03 1.02 0.85 0.91 0.52 0.86 0.42 island-related Shores 1.16 0.075 1.22 0.044 0.86 0.42 0.85 0.29 features Otherregions 0.76 3.5×10−4 0.85 0.06 1.22 0.12 1.19 0.26 3.4. GeneOntologyAnalysisofDifferentiallyMethylatedGenes ToobtaininsightintothepossiblebiologicalimplicationsofthedetectedDNAmchangesinthe OFCofheroinusersandsuicidecompleters,weperformedafunctionalannotationandenrichment analysisoftheDMGsusingtheonlinesoftwaretoolWebGestalt(www.webgestalt.org[79]). Thistool integrates information from a number of public resources to assess the enrichment of gene sets in variousfunctionalcategories,andpresentstheresultsasalistaswellasatree-likediagramdisplaying thehierarchicalrelationshipoftheenrichedcategories. Wefocusedouranalysisonthegeneontology (GO) module of WebGestalt, which comprises 3 large subsets of GO terms: “Biological Process”, “MolecularFunction”and“CellularComponent”. Notably,whenhyperDMhergenes(N=660)andhypoDMhergenes(N=531)wereanalyzed,we obtainedessentiallynon-overlappingsetsofenrichedGOterms(SupplementaryMaterialsTableS6and FigureS4A,B).Amongthe“Cellular_Component”GOsubclass,theGOlocalizationtermsassociated withaxonsorsynapticcompartments(“axon”,“synapticmembrane”)wereenrichedwiththehyperDM butnotwiththehypoDMgeneset. Incontrast,severalGOcategoriesassociatedwithgeneexpression regulationandtranscriptionfactoractivity(“transcriptionfactorbinding”,“sequence-specificDNA bindingRNApolymeraseIItranscriptionfactoractivity”andothers)wereenrichedinthehypoDM butnotinthehyperDMgenesets. Inaddition,“Biological_Process”GOterms“transmissionofnerve impulse”,“axonogenesis”and“cell–cellsignaling”werespecifictothehyperDMgenes. Ontheother hand,thehypoDMgeneswereenrichedforthe“regulationofneurondifferentiation”GOterm,which includedseveraltranscriptionfactorsandgrowthfactors(e.g.,NGF). Whensuicide-associatedDMGswereanalyzed,nosignificantenrichmentswerefoundforthe hypoDMgeneset,whichwasonlycomprisedof170genes. ThehyperDMsuigeneswereenrichedfor the“synaptictransmission”GOterm,aswellasforseveralotherbroadGOtermsnotdirectlyrelated tothenervoussystemfunction(SupplementaryMaterialsFigureS4C). Genes2017,8,152 10of18 3.5. NeuronalSubtype-SpecificityofHyper-orHypomethylatedGenes Our analysis of DM between heroin users or suicide subjects and controls was performed in FANS-separated neuronal nuclei. This experimental design could not distinguish between the contributions of specific neuronal subpopulations that exist in human PFC (e.g., excitatory glutamatergicorinhibitoryGABAergicneurons)totheobserveddifferencesinDNAm. Toinvestigate neuronal subtype-specificity of the hyper- and hypoDM genes, we compared the DMG sets with genes that were differentially expressed (DE) between Glu projection neurons and MGE-derived GABAinterneuronsinthePFC.MGE-derivedGABAinterneuronscomprise~60–70%ofallcortical GABAneurons, andexpressparvalbumin(PVALB)orsomatostatin(SST)[83]. Thispopulationof GABAneuronsconstitutesanessentialcomponentofinhibitoryneuronalnetworksinthehumanand mammaliancortexandhasbeenimplicatedinmanyneurologicalandpsychiatricdiseases(including schizophrenia,majordepression,autism,andepilepsy)[84–87]. Therefore,itwasimportantforthe interpretationoftheresultsofourstudytoclassifytheDMGsasGABA-orGlu-specific. Tothisend, weemployedtheFANS-basednucleiisolationprotocolthatwasrecentlydevelopedinourlab[69]and thatallowedustoseparateandobtainneuronalnucleifromtheseneuronalpopulations. PFCtissues obtainedfrom3controlsubjectswereanalyzed. The DE analysis identified a number of genes that were mostly expressed in Glu or GABA neurons(Glu-specificandGABA-specificDEgenes(DEGs),N =1339and817genes,respectively), Next,weoverlappedtheDEGandDMgenesetsandassessedtheenrichmentusingthehypergeometric distribution test. A significant difference in patterns of enrichment was observed between hyper- and hypoDM genes (Table 2). For both DMher and DMsui genes, the hyperDMGs were strongly enrichedforGlu-specificDEgenes,butnotforGABA-specificgenes(Table2). Incontrast,hypoDMher and hypoDMsui genes showed no enrichment for Glu-specific DEGs. We, however, detected that hypoDMsuigenesweresignificantlyenrichedforGABA-specificDEGs(Table2). Table2. OverlapbetweenheroinorsuicideDMGsandGlu-orGABA-specificDEGs. Theoverlaps werecalculatedseparatelyforhyper-andhypoDMGs. DMGsassociatedwithbothhyperDM-and hypoDMsiteswerenotincluded.Statisticalsignificancewascalculatedwiththehypergeometrictest. Significantoverlaps(p-value<0.01)areshowninbold. Onlygenespresentinboththe450Karray annotationandtheDEanalysisdatasetwereincludedintheanalysis(13,544genes). Glu-DEGenes N=1074Genes GABA-DEGenes N=701Genes DM/DEGenes Overlapping p-Value Overlapping p-Value hyperDMGs N=537genes heroin 70 2.2×10−5 27 0.59 hypoDMGs N=405genes heroin 40 0.087 31 0.019 hyperDMGs N=201genes suicide 34 1.9×10−5 11 0.47 hypoDMGs N=133genes suicide 14 0.17 15 3.6×10−3 4. Discussion InthisstudyweemployedtheIllumina450KInfiniumDNAmethylationarraytoprobethe differencesinDNAmintheOFCofheroinabuserswhodiedofheroinoverdose,suicidecompleters, andcontrolsubjectswhodiedofnon-suicidecauses. Therewasnoevidenceofanyillicitdrugusefor thelattertwogroups. WeanalyzeddifferentialDNAm(1)betweensuicidecompletersandcontrols and (2) between heroin users and controls. Applying a conservative nominal p-value (<0.001) we detected454CpGsitesthatwereDMbetweensuicideandcontrolsubjects,and1298CpGDMsites betweenheroinusersandcontrols. Wethenassessedgenomicfeatures,neuron-subtypespecificityas wellasfunctionalcategoriesandpathwaysthatwereassociatedwiththeseDNAmdifferences.