DNA Methylation Profiling of Human Prefrontal Cortex Neurons in Heroin Users Shows Significant

18 Pages·2017·0.68 MB·English

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.:[phone](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.