Opioid Medication Use and Blood Dna Methylation: Epigenome-Wide Association Meta-Analysis

Abstract

Aim: To identify differential methylation related to prescribed opioid use. Methods: This study examined whether blood DNA methylation, measured using Illumina arrays, differs by recent opioid medication use in four population-based cohorts. We meta-analyzed results (282 users; 10,560 nonusers) using inverse-variance weighting. Results: Differential methylation (false discovery rate <0.05) was observed at six CpGs annotated to the following genes: KIAA0226, CPLX2, TDRP, RNF38, TTC23 and GPR179. Integrative epigenomic analyses linked implicated loci to regulatory elements in blood and/or brain. Additionally, 74 CpGs were differentially methylated in males or females. Methylation at significant CpGs correlated with gene expression in blood and/or brain. Conclusion: This study identified DNA methylation related to opioid medication use in general populations. The results could inform the development of blood methylation biomarkers of opioid use.

Keywords::

• blood
• epigenetics
• methylome
• prescription opioids

Opioid use and dependence are major public health concerns. Opioid medications are prescribed to treat moderate-to-severe pain because they bind to opioid receptors to reduce pain intensity [1]. Approximately 70% of drug overdose-related deaths in the USA involve opioids, including prescription opioids [2]. Improved understanding of biological responses to prescription opioids could reduce these adverse consequences.

DNA methylation is an epigenetic mechanism that can influence gene expression without changing the underlying genetic sequence. Epigenome-wide association studies (EWAS) have shown that DNA methylation in blood can be altered by substance use such as cigarette smoking [3] and alcohol consumption [4]. The extensive methylation signals for these substances have led to the development of biomarkers of habitual use [4,5].

Some small to moderately sized studies (n ≤ 370) have reported methylation alterations related to opioid dependence, defined by the Diagnostic and Statistical Manual of Mental Disorders, fourth edition, criteria [6–8], or to opioid use [9,10]. Candidate gene studies of the OPRM1 gene reported hypermethylation related to opioid dependence or use [6,7,9–11]. An EWAS identified three CpG probes meeting genome-wide significance related to opioid dependence [8]. To date, there have been no large-scale meta-analyses of genome-wide blood DNA methylation and opioid medication use.

Sex differences in pain and analgesia have been reported [12,13]. A recent large study found that men had a greater risk than women for opioid overdose after their initial prescription [14]. In one study of opioid prescriptions to treat noncancer-related pain, men were more likely than women to escalate into higher-dose treatments and die of opioid-related causes [15]. Two of the five candidate gene methylation studies focusing on the OPRM1 gene were limited to men (n < 120) [6,11], and the sole EWAS reporting three CpGs of genome-wide significance was conducted only in women (n = 220) [8]. No studies have reported well-validated DNA methylation alterations related to opioid use using large-population studies or exploring whether associations between opioid use and DNA methylation differ by sex.

This study focused on the use of opioid medications and excluded illicit opioid use. Because any intake of opioid medications could provide an individual a chance to become a chronic user or develop a misuse disorder, generating an accurate biomarker to help identify individuals exposed to opioid medications is essential in opioid- and/or health effects-related research.

We conducted a coordinated, epigenome-wide association meta-analysis across five datasets (10,842 individuals comprising 282 opioid users and 10,560 nonusers), including four of European ancestry and one of African ancestry, from cohorts in the USA and Scotland to identify CpGs differentially methylated in relation to recent opioid medication use. Additionally, we explored whether opioid-related differential methylation differs by sex. To examine the functional implications of identified CpGs, we evaluated associations with nearby gene expression in blood (n = 3684) and brain (n = 515) and overlaps with histone marks and transcription factor (TF) binding motifs. This study of blood DNA methylation in general population studies could contribute to biomarker development and an improved understanding of the epigenomic mechanisms of response to opioid medications.

Materials & methods

Study population

This meta-analysis included 10,842 adults (≥30 years old) from five datasets, within the framework of the Cohorts for Heart and Aging Research in Genomic Epidemiology consortium. Participating studies included the Agricultural Lung Health Study (ALHS), Framingham Heart Study (FHS), Generation Scotland (GS Set1 and Set2), and Genetic Epidemiology Network of Arteriopathy (GENOA). Study participants were mainly of European ancestry (91% of the 10,842 adults). African-ancestry individuals were from GENOA. Each study was approved by its institutional review board. All study participants provided written informed consent.

Recent opioid medication use

The harmonized definition across studies categorized any recent use of opioid medications as exposed (users). In brief, ALHS participants, enrolled between 2009 and 2013, were asked at home visits to provide all prescription and nonprescription medications used in the last 7 days and field technicians recorded the medication names, which were coded using the Slone Drug Dictionary [16], which links drug names to active ingredients and drug classes based on Pharmacologic-Therapeutic Classification from the American Hospital Formulary Service. Recent opioid use was classified based on any use of full opioid agonists in the past 7 days versus no use. In the FHS, participants were asked to bring all prescription and nonprescription medications taken regularly in the past month in their original containers to the study visit. Medication names were converted into Anatomical Therapeutic Chemical codes. Recent opioid use (users vs nonusers) was based on any use of full opioid agonists in the past month. In GS, participants reported prescription and nonprescription medications they regularly took. Any use of opioid medications was employed to define recent opioid use. In GENOA, participants were asked “During the last month, have you used any medication that was prescribed or recommended by a physician?”; if they responded ”yes”, they provided the names of the medications. The Medi-Span Therapeutic Classification was used to identify pharmacologic classes of prescriptions. Participants reporting any use of opioid medications were categorized as recent opioid users. In each study, participants who did not report any use of opioids were categorized as unexposed (nonusers).

DNA methylation

We measured genome-wide methylation in blood, collected at the time of the medication questionnaires, using Illumina methylation arrays. FHS used the Illumina HumanMethylation450 BeadChip (450K array; Illumina, Inc., CA, USA) and the other studies used the newer MethylationEPIC BeadChip (EPIC/850K array; Illumina, Inc.). Each study performed preprocessing and quality control filtering, including the removal of methylation values based on detection P, exclusion of poor-quality CpGs and samples, probe-type bias normalization and correction for batch effects. Additionally, we excluded CpGs previously reported as potentially problematic [17], including “ch” probes, probes having an SNP in the extension base or color channel issues, and cross-reactive probes. Autosomal CpGs that passed these quality control steps were included in association analyses.

Differential methylation in relation to opioid use

Each study assessed associations of recent opioid medication use (opioid users vs nonusers) as the predictor with methylation, ranging between 0 (unmethylated) and 1 (methylated) as the response, using robust linear regression to account for potential heteroskedasticity and lessen the impacts of influential outliers on association results. To further reduce potential impacts of extreme outliers, we used winsorization (0.25% on both sides) [18]. Each study adjusted for potential confounders, including age, sex, BMI, smoking (never/former/current as dummy variables), smoking pack-years, ancestry principal components based on genome-wide genetic variants and estimated cell-type proportions (monocyte, B cell, natural killer cell, CD8⁺ T cell, CD4⁺ T cell and granulocyte) [19]. Additionally, the studies adjusted for analytic batch, study site and any study-specific selection factor or accounted for family relationships where appropriate. We meta-analyzed study-specific results using inverse-variance weighted fixed-effects models [20,21], including 778,492 CpGs analyzed in two or more datasets.

Additionally, we explored whether differential methylation related to opioid use differed by sex using studies having more than 15 opioid users in both sex groups, which limited the analysis to ALHS and FHS. Because large-sample statistics were used, we set a threshold of 15 in the exposed group within each analysis based on prior experience in consortium epigenomic meta-analyses. We analyzed the data separately in males and females using robust linear regression, adding a product term for opioid use × sex to a model with both main effects in all participants and with the same covariates as above. We meta-analyzed the results using inverse-variance weighting [20,21], including the 394,524 CpGs overlapping between EPIC/850K array (ALHS) and 450K array (FHS).

To account for multiple testing, we set a genome-wide significance threshold of a Benjamini–Hochberg false discovery rate (FDR) of 0.05 [22]. Unless otherwise noted, CpGs of genome-wide significance refer to those with FDR < 0.05.

We also evaluated differentially methylated regions (DMRs) in relation to opioid use using the comb-p method [23] implemented in the ENmix R package [24]. EWAS meta-analysis results were used as the input. Comb-p combines irregularly spaced adjacent p-values to identify regions containing consecutive low p-values [23]. Regional p-values were corrected for multiple testing using Šidák (default option) [25]. We set a threshold of Šidák P < 0.05 for significant genome-wide DMRs. The DMR analysis complements individual CpG level analysis, as the method aggregates signals across multiple nearby CpGs and a DMR of statistical significance does not necessarily contain a CpG having genome-wide significance.

Functional & biological implications

Using the eFORGE integrative epigenomic approach [26,27], we explored whether our top 1000 CpGs (default option) were enriched in regulatory elements from the Roadmap Epigenomics Consortium [28] across over 20 tissue types. To evaluate overlaps between the significant CpGs and TF binding motifs in the brain, we used eFORGE-TF [26]. Additionally, we examined enrichment of TF motifs in implicated loci (±250 bp) using MEME Suite [29].

To further evaluate the biological context of the methylation findings, we tested enrichment of pathways using robust rank aggregation [30] and preranked gene set enrichment analysis [31], which account for differing numbers of CpGs per gene on the methylation arrays. We used the methylRRA function in the methylGSA [32] R package and reported pathways enriched at nominal significance (p < 0.05). Pathways were from the Kyoto Encyclopedia of Genes and Genomes gene sets via the org.Hs.eg.db [33] R package.

Opioid-related CpGs & nearby gene expression: cis expression quantitative trait methylation analysis

To evaluate functional impacts of genome-wide significant CpGs, we assessed whether methylation at identified CpGs was associated with expression of nearby (±250 kb) genes: cis-expression quantitative trait methylation (eQTM). Linear regression was used to examine associations between methylation (predictor) and the expression of nearby genes (response). To identify cis-eQTMs in blood, we used paired DNA methylation (450K) and gene expression (Affymetrix Human Exon 1.0 ST GeneChip Platform) from FHS [34] (n = 3684). For cis-eQTMs in the brain, we used paired DNA methylation (450K) and gene expression (Illumina HiSeq) from ROS/MAP [35] (n = 515). We set a threshold of FDR < 0.05 for statistical significance.

Opioid-related CpGs & nearby genetic variants: cis methylation quantitative trait loci analysis

To examine whether genome-wide significant CpGs in all participants showed associations with nearby (±500 kb) genetic variants, we used the GoDMC database (N = 32,851 individuals) [36] of methylation quantitative trait loci (meQTL).

Look-up of previously reported opioid-related methylation in our data

Because few data exist on genome-wide DNA methylation and opioid use, we attempted to validate findings from a study using saliva samples [10] in our analysis of blood. We looked up CpGs reported (p < 0.0001) from a recent EWAS of prescription opioid use [10]. We also examined whether CpGs previously related to opioid dependence [8] were related to opioid use in our data. Additionally, we interrogated CpGs annotated to the OPRM1 gene in our data because candidate gene studies reported hypermethylation related to opioid addiction and use [7,10].

Additional analyses

We considered potential residual confounding by substances known for their impact on genome-wide blood DNA methylation. Although we had adjusted for smoking status (never/past/current) and lifetime intensity (pack-years), to evaluate possible residual confounding by smoking, we examined whether CpGs identified for opioid use (FDR < 0.05) in our data overlapped with CpGs associated with current smoking in a previous large EWAS meta-analysis [3]. For alcohol use, we were not able to harmonize and adjust for alcohol consumption due to minimal information on alcohol use in the ALHS, one of the larger studies. However, we investigated whether CpGs identified (FDR < 0.05) in our meta-analysis overlapped with CpGs related to heavy alcohol consumption in a previous EWAS meta-analysis [4]. Because we included one African-ancestry study (GENOA) in our meta-analysis of all participants, we conducted a leave-one-out meta-analysis to verify that our findings (FDR < 0.05) were not driven by this single study of non-European ancestry.

Results

This meta-analysis included 10,842 participants from five datasets (Table 1). On average, participants were 58 years old. About 3% (282/10,842) reported recent use (up to 30 days) of opioid medications, and this proportion did not differ by sex (males: 83/3,108; females: 117/3,348). We visualized meta-analysis results in Miami plots (Figure 1) [37] and quantile-quantile plots (Supplementary Figure 1) with λ ranging between 1.1 and 1.3, suggesting minimal inflation.

Table 1. Characteristics of participating studies.

	ALHS			FHS			GS Set1^‡	GS Set2^‡	GENOA^‡
	All	Males	Females	All	Males	Females	All	All	All
Ancestry	European			European			European	European	African
Country	USA			USA			UK	UK	USA
Characteristics, mean ± standard deviation or N (%)
N total	2286	1173	1113	4170	1935	2235	1864	1566	956
N exposed^† (%)	130 (6%)	53 (5%)	77 (7%)	70 (2%)	30 (2%)	40 (2%)	36 (2%)	20 (1%)	26 (3%)
Age, years	63 ± 11	63 ± 11	62 ± 11	59 ± 13	59 ± 13	59 ± 13	51 ± 11	54 ± 10	58 ± 10
BMI, kg/m^2	30 ± 6	30 ± 5	30 ± 7	28 ± 5	29 ± 5	27 ± 6	27 ± 5	27 ± 5	31 ± 7
Smoking status
– Never	1523 (67%)	680 (58%)	843 (76%)	1790 (43%)	742 (38%)	1048 (47%)	935 (50%)	813 (52%)	573 (60%)
– Former	667 (29%)	429 (37%)	238 (21%)	1940 (47%)	952 (49%)	988 (44%)	585 (31%)	527 (34%)	223 (23%)
– Current	96 (4%)	64 (5%)	32 (3%)	440 (11%)	241 (12%)	199 (9%)	344 (18%)	226 (14%)	160 (17%)
Pack-years among ever smokers	18 ± 21	21 ± 23	12 ± 17	12 ± 19	12 ± 20	12 ± 19	18 ± 17	16 ± 17	21 ± 18
Methylation array	EPIC/850K			450K			EPIC/850K	EPIC/850K	EPIC/850K
N CpGs analyzed	817,235	817,235	817,235	473,864	473,864	473,864	841,753	773,860	838,270
λ inflation factor	1.23	1.11	1.16	0.92	1.27	0.79	0.87	0.97	1.12

† N exposed to opioid medications.

‡ Excluded from meta-analyses stratified by sex due to individuals exposed to opioids being ≤15 in either males or females. GS Set1 included 697 males and 1167 females, and GS Set2 contained 671 males and 895 females. Of GENOA participants, 276 were males and 680 were females.

Figure 1. In a Miami plot, each dot represents -log10(P) of a single CpG.

Each plot has two panels: upper for association results with positive regression coefficients and lower for association results with negative regression coefficients, with -log10 (p) on the Y axis and 22 chromosomes on the X axis. Horizontal lines depict p cutoffs for statistical significance after multiple-testing correction: Bonferroni and Benjamini–Hochberg false discovery rate. CpGs having uncorrected p > 0.05 were not displayed. (A) Epigenome-wide association study (EWAS) meta-analysis results in all participants. (B) EWAS meta-analysis results in males only. (C) EWAS meta-analysis results in females only.

We identified six CpGs differentially methylated in relation to recent opioid medication use (FDR < 0.05; Table 2). Of the six CpGs, four showed higher methylation in opioid users compared with nonusers. The CpGs annotated to six different genes: RUBCN (KIAA0226), CPLX2, TDRP, RNF38, TTC23 and GPR179. We did not observe significant heterogeneity across studies in the genome-wide significant findings (Supplementary Figure 2). Leave-one-out meta-analysis results confirmed that signals at the six CpGs were not driven by GENOA, a study of African ancestry (Supplementary Table 1).

Table 2. Six CpGs differentially methylated (false discovery rate < 0.05) in relation to opioid medication use (N_exposed = 282; N_unexposed = 10,560).

Chromosomal position	Probe	Effect^†	Standard error^‡	p-value	Direction^§	Gene name^¶	Mean methylation^#
3: 197459812	cg07255038	0.0057	0.0011	3.21 × 10^-7	++-++	KIAA0226	0.8602
5:175223982	cg07295964	-0.0016	0.0003	8.23 × 10^-8	---+-	CPLX2	0.0388
8:528687	cg23406159	0.0061	0.0012	1.46 × 10^-7	+++++	TDRP	0.7650
9:36479098	cg09409129^††	0.0016	0.0003	3.65 × 10^-7	+?--+	RNF38	0.9187
15:99789777	cg13401703	-0.0109	0.0021	2.01 × 10^-7	-+--+	TTC23	0.2523
17: 36499760	cg23186467^††	0.0059	0.0012	3.63 × 10^-7	+?+++	GPR179	0.7501

† Regression coefficient from statistical models.

‡ Standard error of regression coefficient.

§ Direction of associations. Order of studies is ALHS, FHS, GS Set1, GS Set2 and GENOA.

¶ Gene name listed on the Illumina annotation. When no gene names were reported on the Illumina annotation, the authors used Zhou et al. [16] and Homer v4.9.1 [50] to annotate genes.

# Weighted average methylation across participating studies.

†† CpG present only on the EPIC array; 212 exposed/6,460 unexposed.

In sex-stratified meta-analyses, we found 74 CpGs related to opioid use (FDR <0.05): 68 CpGs in males (top 30 CpGs in Table 3; all 68 CpGs in Supplementary Table 2) and six CpGs in females (Table 4); one additional CpG, cg07295964 (CPLX2), significant in females but not males, was also identified in the meta-analysis of all participants. The directions of associations were skewed toward positive in both sexes. Although none of the sex-specific findings showed genome-wide significant interactions between opioid use and sex, 51% showed at least nominal evidence of interaction (P_interaction < 0.05; Supplementary Table 2 & Tables 3 & 4), much more than expected by chance (P_enrichment < 2.2 × 10^-16).

Table 3. Top 30 differentially methylated CpGs (false discovery rate < 0.05) in relation to opioid medication use in males (N_exposed = 83; N_unexposed = 3025), sorted by chromosomal location.

Chromosomal position	Probe	Effect^†	Standard error^‡	p-value	Direction^§	Gene name^¶	Mean methylation^#	Pinteraction^††
1:32083031	cg11573219	0.0287	0.0052	3.76 × 10^-8	++	HCRTR1	0.3369	0.0001
1:42385662	cg16145216	0.0328	0.0065	4.78 × 10^-7	++	HIVEP3	0.3683	0.0816
1:156094519	cg17453456	0.0251	0.0045	2.30 × 10^-8	++	LMNA	0.3100	0.0744
1:240655362	cg01176047	0.0147	0.0028	2.07 × 10^-7	++	GREM2	0.8172	0.2044
2:208635723	cg18866599	-0.0013	0.0002	3.27 × 10^-8	-+	FZD5	0.0463	0.0002
3:29494928	cg26400275	0.0057	0.0010	3.75 × 10^-8	++	RBMS3	0.8773	0.0911
3:58508149	cg23652987	0.0091	0.0019	8.95 × 10^-7	++	ACOX2	0.8186	0.9107
4:47841314	cg02955940	0.0082	0.0017	1.24 × 10^-6	++	CORIN	0.8380	0.7118
4:74486177	cg07136054	-0.0010	0.0002	7.42 × 10^-8	-+	RASSF6	0.0398	8.62 × 10^-6
5:68788113	cg16151977	-0.0038	0.0007	2.72 × 10^-7	--	OCLN	0.0537	0.0111
5:178424364	cg14800603	0.0125	0.0025	5.14 × 10^-7	++	GRM6	0.8196	0.2742
5:179246250	cg25991122	-0.0016	0.0003	3.19 × 10^-7	--	SQSTM1	0.0529	0.0318
5:180633772	cg25566730	-0.0004	0.0001	2.16 × 10^-7	--	CTC-338M12.1, TRIM7	0.0319	0.0020
6:6005992	cg19060550	-0.0011	0.0002	4.51 × 10^-7	-+	NRN1	0.0355	0.0036
6:27513092	cg11497372	-0.0015	0.0003	5.52 × 10^-9	--	AL021918.1	0.0723	0.0022
6:28911926	cg15088539	-0.0068	0.0014	6.85 × 10^-7	--	LINC01556	0.0828	0.0017
6:150311552	cg25106322	-0.0010	0.0002	8.41 × 10^-7	-+	RAET1K	0.0389	0.0002
7:2106321	cg12098228	-0.0078	0.0016	6.87 × 10^-7	--	MAD1L1	0.0882	0.5122
8:97505764	cg13096260	-0.0010	0.0002	6.57 × 10^-7	--	SDC2	0.0690	0.0004
9:139871049	cg00563932	0.0199	0.0041	1.09 × 10^-6	++	PTGDS	0.4954	0.0080
10:79212324	cg18659450	0.0056	0.0011	4.84 × 10^-7	++	KCNMA1	0.8915	0.6008
11:33279091	cg04704963	0.0015	0.0003	2.00 × 10^-7	++	HIPK3	0.0291	0.0190
12:52603903	cg09627757	0.0050	0.0011	2.38 × 10^-6	++	LOC283404	0.8854	0.2494
12:63025730	cg08244085	0.0109	0.0022	5.70 × 10^-7	++	MIRLET7I	0.1085	0.0009
12:100967381	cg13069247	-0.0025	0.0005	4.92 × 10^-7	--	GAS2L3	0.0659	0.0095
12:102513805	cg06192289	0.0032	0.0007	1.14 × 10^-6	++	NUP37, C12ORF48	0.0771	0.0259
12:130524707	cg08968978	0.0142	0.0029	1.24 × 10^-6	++	LOC100190940	0.8470	0.6107
16:68024755	cg03384419	0.0027	0.0006	1.66 × 10^-6	++	DPEP2	0.9156	0.0380
16:68624584	cg10486655	-0.0010	0.0002	2.21 × 10^-6	-+	RP11-615I2.1, RP11-615I2.6	0.0442	0.0008
17:20905737	cg11691341	0.0071	0.0015	1.60 × 10^-6	++	USP22	0.8771	0.3097

† Regression coefficient from statistical models.

‡ Standard error of regression coefficient.

§ Direction of associations. Order of studies is ALHS and FHS.

¶ Gene name listed on the Illumina annotation. When no gene names were reported on the Illumina annotation, the authors used Zhou et al. [16] and Homer v4.9.1 [50] to annotate genes.

# Weighted average methylation across participating studies.

†† P for interaction between opioid use and sex.

Table 4. Seven CpGs differentially methylated (false discovery rate <0.05) in relation to opioid medication use in females (N_exposed = 117; N_unexposed = 3231), sorted by chromosomal location.

Chromosomal position	Probe	Effect^†	Standard error^‡	p-value	Direction^§	Gene name^¶	Mean methylation^#	Pinteraction^††
1:14033779	cg10375300	0.0070	0.0014	2.80 × 10^-7	++	PRDM2	0.8890	0.1614
3:3840360	cg12136608	-0.0009	0.0002	2.49 × 10^-7	-+	LRRN1	0.0382	0.0049
5:139485587	cg10872490	0.0035	0.0007	1.29 × 10^-7	+-	LINC01024	0.9049	0.0023
5:175223982	cg07295964^‡‡	-0.0019	0.0004	3.37 × 10^-7	--	CPLX2	0.0375	0.2006
7:157955692	cg02161503	0.0088	0.0017	1.11 × 10^-7	++	PTPRN2	0.8582	0.5626
12:55029660	cg19880901	0.0044	0.0008	1.35 × 10^-7	++	LACRT	0.8754	0.9825
13:114085732	cg03049340	0.0048	0.0009	2.31 × 10^-7	++	ADPRHL1	0.8734	0.0022

† Regression coefficient from statistical models.

‡ Standard error of regression coefficient.

§ Direction of associations. Order of studies is ALHS and FHS.

¶ Gene name listed on the Illumina annotation. When no gene names were reported on the Illumina annotation, the authors used Zhou et al. [16] and Homer v4.9.1 [50] to annotate genes.

# Weighted average methylation across participating studies.

†† P for interaction between opioid use and sex.

‡‡ CpG showing genome-wide significance in meta-analysis of all participants.

In regional analyses of all participants, 20 regions were differentially methylated in relation to opioid use (Šidák P < 0.05; Supplementary Table 3). One DMR (Chr15:99789621-99790022) included a significant genome-wide CpG (cg13401703). Sex-stratified analyses identified 24 DMRs for males and 15 for females (Supplementary Table 4). One additional region, Chr6:31650734-31651291 (LY6G5C), identified in females but not males showed genome-wide significance in analysis of all participants. Just four DMRs identified in sex-stratified analyses contained a significant genome-wide CpG (Supplementary Table 5). Table 5 displays the top ten DMRs in all participants and by sex.

Table 5. Top ten differentially methylated regions related to opioid medication use (Šidák p < 0.05), sorted by chromosomal location.

Chromosome	Start^† (bp)	End (bp)	Šidák P	N CpGs^‡	Minimum P^§	Annotated gene^¶
All participants
1	243876785	243877024	0.0011	3	5.77 × 10^-5	AKT3
5	134526123	134526269	0.0001	4	2.73 × 10^-6	C5ORF66
6	29648160	29648756	3.86 × 10^-8	20	0.0022	ZFP57
6	31650734	31651291	5.18 × 10^-10	18	0.0015	LY6G5C
6	31683050	31683352	1.78 × 10^-6	6	2.59 × 10^-6	LY6G6E, LY6G6D
15	99789621	99790022	2.33 × 10^-11	8	2.01 × 10^-7	TTC23
16	450836	451040	1.71 × 10^-7	3	2.98 × 10^-5	DECR2
17	76037034	76037364	0.0002	5	0.0031	TNRC6C
17	81060148	81060259	1.99 × 10^-5	3	0.0002	AC144831.1
22	38092718	38092989	0.0010	8	0.0021	TRIOBP
Males only
1	42385580	42385941	1.49 × 10^-8	3	4.78 × 10^-7	HIVEP3
6	28911872	28912166	1.02 × 10^-5	8	6.85 × 10^-7	LINC01556
6	31543539	31543686	4.91 × 10^-6	8	0.0017	TNF
6	84140682	84140997	6.75 × 10^-5	7	0.0004	ME1
7	73720807	73721031	1.22 × 10^-5	3	4.09 × 10^-5	CLIP2
8	22735110	22735478	6.16 × 10^-6	5	0.0019	PEBP4
17	17465161	17465476	1.65 × 10^-5	5	0.0002	PEMT
17	80859885	80860250	7.14 × 10^-7	3	3.38 × 10^-5	TBCD
17	80870528	80870577	0.0002	2	4.87 × 10^-6	TBCD
21	40759533	40759694	7.15 × 10^-7	5	0.0009	WRB
Females only
1	236686563	236686804	1.54 × 10^-7	4	8.32 × 10^-5	LGALS8, LGALS8-AS1
6	11779816	11779941	0.0011	3	5.45 × 10^-6	C6ORF105
6	29521012	29521272	0.0008	11	0.0011	OR2I1P
8	39172019	39172120	0.0004	6	0.0020	ADAM5P
11	2020278	2020560	0.0003	10	0.0008	H19
11	70517293	70517411	0.0005	3	0.0001	SHANK2
17	75143117	75143184	0.0010	2	3.14 × 10^-5	SEC14L1
17	76037034	76037562	3.75 × 10^-9	6	0.0001	TNRC6C
19	19281174	19281474	0.0009	7	0.0012	LOC729991-MEF2B, MEF2B
21	35831870	35832180	0.0009	9	0.0015	KCNE1

† Physical position (base pair), National Center for Biotechnology Information human reference genome assembly Build 37.3.

‡ Number of CpGs in the region.

§ Minimum p-value among uncorrected p-values of CpGs in each region.

¶ Gene name listed on the Illumina annotation. When no gene names were reported on the Illumina annotation, the authors used Zhou et al. [16] and Homer v4.9.1 [50] to annotate genes.

Integrative epigenomic analyses using eFORGE identified enrichment for active transcriptional histone mark H3K36me3, which recruits DNA repair machinery upon damage to maintain genomic stability [38], in several tissue types, including blood and brain (Supplementary Figure 3), supporting potential functional implications. In TF binding motif analyses using eFORGE-TF, we found an overlap with the motif V_ZID_01, corresponding to the transcription factor ZBTB6, in brain (p = 0.00002). An additional TF motif analysis using Analysis of Motif Enrichment [39] identified the motif ZBT18_HUMAN.H11MO.0.C near all six significant CpGs but only 22 of 100 random nonsignificant CpGs (p = 0.01).

Pathway analysis identified enrichment for 19 gene sets, including the chemokine signaling pathway, a mechanism related to pain modulation. Most pathways were related to metabolism and inflammation (Supplementary Table 6).

By linking the significant CpGs (FDR < 0.05 in either the all-participants or sex-specific analyses) to transcriptome data, we identified CpGs related to the expression of nearby (±250 kb) genes: cis-eQTMs. In blood (FHS, n = 3684) [34], three CpGs showed genome-wide significant associations with gene expression (FDR < 0.05). At nominal significance (p < 0.05), 29 CpGs showed associations (Supplementary Table 7). In brain (ROS/MAP, n = 515), four of 45 CpGs having one or more nearby transcripts were related to nearby gene expression at FDR < 0.05, and 11 (24%) showed nominally significant (p < 0.05) associations with expression (Supplementary Table 8). Six CpGs were related to the expression of nearby genes in both blood and brain.

In the GoDMC database [36] of paired genetic and Illumina450K methylation data in 32,851 individuals, of the six CpGs significant in all participants, cg13401703 (TTC23) was associated with two nearby (both within 145 kb) genetic variants (meQTLs): chr15:99642518 and chr15:99932797.

Given the paucity of data on epigenome-wide blood DNA methylation and opioid use, we looked up findings from a recent EWAS [8] (N = 220) of opioid dependence. Of the three CpGs meeting genome-wide significance in the previous study, two were marginally related to opioid use in our data (P = 0.07 for cg21381136; P = 0.04 for cg17426237). Of the remaining 499 CpGs reported at a threshold of p < 0.0001 and available in the data, 69 (14%) were related to opioid use at p < 0.05 (P_enrichment = 4.09 × 10^-14; Supplementary Table 9). Because several candidate gene studies reported higher methylation at the OPRM1 gene with opioid dependence and use [6,7,9–11], we interrogated 61 CpGs annotated to the gene available in our data after CpG filtering. Of these, four were at least nominally (p < 0.05) associated with opioid use in our meta-analysis of all participants. Another three met p < 0.05 in sex-stratified meta-analyses. Finding these seven CpGs is slightly more than would be expected by chance (P_enrichment = 0.032; Supplementary Table 10).

Although the meta-analysis results were adjusted for smoking, we assessed possible residual confounding by smoking by examining whether the opioid-related CpGs, in either overall or sex-stratified analyses, were associated with current smoking in a previous large Illumina450K EWAS meta-analysis [3]. Of the 78 CpGs (450K) available for look-up, only 11 were associated with current smoking (Supplementary Table 11). We also addressed potential residual confounding by alcohol consumption, which we did not adjust for in the meta-analyses due to limited alcohol intake data in ALHS, by looking up the significant CpGs in a previous large 450K EWAS meta-analysis of heavy alcohol use [4]. Of the 78 CpGs available for look-up, only nine CpGs were related to heavy alcohol consumption (Supplementary Table 11). Just four CpGs were associated with both smoking and alcohol consumption in these look-ups. None of the four genome-wide significant 450K CpGs in analyses of all participants were related to smoking or alcohol consumption.

Discussion

To our knowledge, this is the first large-scale, epigenome-wide association meta-analysis of blood DNA methylation and recent opioid medication use in general population studies. We identified several CpGs significantly differentially methylated in relation to opioid use in all participants. Additionally, we found CpGs showing differential methylation specific to men or women. Our findings were novel, meaning CpGs and genes identified in this study were not previously reported in studies of DNA methylation and opioids. Notably, our findings included CpGs related to the expression of nearby genes (cis-eQTM) in blood and brain. Integrative epigenomic results linked our methylation findings to enrichment of an active histone mark and an overlap with a transcription factor motif, suggesting functional implications. Pathway analysis results highlighted the chemokine signaling pathway, an emerging therapeutic option for pain management [40,41]. In addition to identifying novel methylation signals in relation to opioid use, our data provide the first validation of some previous opioid-related methylation findings in the general population. Extending knowledge from previous studies that focused on opioid dependence, this comprehensive investigation of blood DNA methylation and opioid medication use in the general population will help elucidate molecular mechanisms and population-level impacts of these commonly prescribed medications.

Notably, the opioid-related CpGs identified in the analyses of all participants implicate genes linked to brain development or brain disorders. Three of the six genes identified, KIAA0226, CPLX2 and GPR179, contain one or more CpGs associated with fetal brain age (n = 179 samples), suggesting a role in fetal brain development (Supplementary Table 12) [42]. A genetic mutation in KIAA0226 (RUBCN) was related to the developmental disorder ataxia [43,44]. Genetic variants in the CPLX2 gene have been related to cognitive function in schizophrenia [45]. Genetic variants in TDRP have been associated with attention deficit hyperactivity disorder [46]. Another of the differentially methylated genes, RNF38, has been shown in animal studies to be regulated by the serotonin system in the brain [47,48], which modulates behaviors and psychological functions. Finally, the gene TTC23, which was highlighted in our differential methylation analysis and in our analysis linking methylation to nearby gene expression, was linked to schizophrenia in a recent EWAS [49].

Integrative epigenomic analyses highlighted enrichment of an active transcriptional histone mark, H3K36me3, in both blood and brain. This suggests that the blood findings could help inform epigenomic mechanisms in brain. Additionally, the TF ZBTB18 of the motif ZBT18_HUMAN.H11MO.0.C, common across all six significant CpGs in the overall analysis, has been shown in mice to play a pivotal role in normal brain development through the time-sensitive repression of several proneurogenic genes [50].

We also performed sex-specific analyses because some previous reports suggested differential effects of opioids by sex and some methylation studies analyzed only one sex. We found CpGs showing significant associations in males or females, suggesting some sex differences in the effects of opioids on methylation. We acknowledge that none of the identified CpGs showed genome-wide significant interactions. However, much larger sample sizes are needed to fully examine interactions, and sex-specific methylation findings may shed light on sex differences in opioid use and health effects.

Data on blood DNA methylation and opioid use are sparse. Several studies have examined DNA methylation in relation to opioid dependence. These studies differed from ours because they compared individuals with opioid dependence with individuals exposed to opioids but not dependent [8], analyzed data with small-to-moderate sample sizes [9] or were restricted to either women [8] or men [11]. Most used candidate gene approaches [6,7,9–11], except for one EWAS [8] of opioid dependence in 220 European-ancestry women. Two candidate gene studies (both n ≤ 123) examined methylation of the OPRM1 gene in relation to opioid use: long-term (for at least 1 year prior to sample selection) opium use using blood samples [11] and short-term prescription opioid use using saliva samples [10]. Compared with studies of individuals with opioid dependence in clinical settings, our study of opioid medication use in the general population and molecular impacts has value for deeper understanding of the population-level impacts of opioid medication use.

When studying the biological effects of opioid use, the primary tissue of interest is brain. However, data on DNA methylation in brain are understandably very limited. When developing biomarkers, blood, an easily accessible tissue, is useful. To help interpret the blood methylation at implicated CpGs, we evaluated correlations between methylation in blood and brain [51]. Of four CpGs differentially methylated in relation to opioid use in all participants and available in the blood and brain methylation correlation tool [51] (450K), methylation in blood at cg13401703 (TTC23) was highly correlated with methylation at four brain regions (Pearson correlation: 0.63 to 0.71; Supplementary Figure 4). The observed high correlation suggests that these findings in blood may indicate methylation alterations in the brain.

This study has some limitations. Although the authors attempted to harmonize the definition of recent opioid medication use, it differed by study, and the exact time course or dosage of ingestion was unknown. This is because this meta-analysis combined data from population-based studies in which inventories of recently used medications were available. In ALHS, we asked about medications in the past 7 days, the shortest period of any of the studies. The reported recent use in our participating studies is likely a mix of short- and long-term use. Longer use would be predicted to have greater impacts on DNA methylation, and the study may have missed these or underestimated them. However, in ALHS, we did not visit participants with abdominal or chest surgery in the past 3 months, which will have reduced short-term postoperative prescriptions. There were no data to investigate whether the reasons for taking opioids differ between sexes, which might contribute to sex-specific findings. The study participants were predominantly older individuals. Indications for time-limited prescription opioid use are likely to change over the life course. The participating studies are population-based studies of free-living people. Blood was collected during visits specifically conducted for the study – not during visits to the hospital. Population-based studies of older people will include people with various age-related conditions. However, we find that if participants are feeling acutely ill on a given day or so chronically disabled that participation is difficult, they will not agree to the study visit. It is standard in population-based studies not to exclude people based on chronic disease or age. However, it is standard to carefully adjust for age in all analyses. We were not able to replicate these meta-analysis findings in an independent population, as the analysis included as many population-based studies as we were able to identify with the available data. Instead, we provide evidence of validation for some previously reported methylation signals related to opioid dependence or use in our data.

This study has several strengths. To our knowledge, this is the largest study of genome-wide DNA methylation and opioid medication use. We carried out sex-specific analyses and included studies of African ancestry and European ancestry. We identified several DMRs. By linking our methylation findings to transcriptome in blood and brain, we found CpGs with potential functional impacts. Integrative epigenomic results link our data to regulatory elements, including histone marks and TF motifs, in blood and/or brain. Because this meta-analysis included individuals participating in population-based studies not enriched for opioid dependence, the results are particularly useful for the development of biomarkers of opioid medication use.

Publication bias is an issue for meta-analyses involving published studies. However, this is a de novo meta-analysis where none of the individual studies had previously published on this association. Therefore, publication bias is not an issue in this study.

In population-based studies, it is expected there will be a low proportion of opioid users and thus that most of the participants are in the unexposed group. The power of the study comes primarily from the number of exposed, but the imbalance with the expected large proportion of unexposed in these studies is not a source of bias.

To better understand health effects of an exposure, objective measures of exposure are useful. Some epidemiologic studies may have asked participants to self-report use of opioid medications but likely have limited information on duration or dosage. Where available, linkage to prescription registries can be used but lack information on what is actually taken. Furthermore, participants may be reluctant to acknowledge the use of opioids, leading to underreporting. Given that many studies have EWAS data but limited data on medication use, a well-validated biomarker of opioid use would increase the usefulness of existing studies to elucidate the full range of the health effects of this exposure. In clinical settings, objective biomarkers could improve the identification of individuals using these medications and help clinicians manage opioid prescriptions and alternative treatment options for chronic pain.

Conclusion

In summary, this large-scale epigenome-wide association study of blood DNA methylation and recent opioid medication use extends the literature by greatly increasing the sample size and including general-population studies with recent genome-wide DNA methylation arrays in two ancestries. This epigenome-wide meta-analysis identifies loci that have not previously been linked to opioid use. The findings could serve as a resource that can be leveraged for the future development of epigenomic biomarkers of opioid use and to extend our understanding of the effects of these drugs.

Summary points

• This is the first epigenome-wide association meta-analysis of opioid medication use and blood DNA methylation.

• The study identified several CpGs differentially methylated in relation to opioid medication use in all participants.

• There were many differential methylation signals specific to males or females.

• Using integrative epigenomic tools, we found the differential methylation signals overlap with regulatory elements in blood and/or brain.

• Methylation levels at significant CpGs were correlated with the expression of nearby genes in blood and/or brain, furthering functional relevance.

Author contributions

M Lee and SJ London conceived and designed the study. M Lee, R Joehanes, DL McCartney and M Kho completed study-specific analyses. M Lee, R Joehanes, A Hüls, S Sikdar and A Burkholder performed follow-up analyses. AB Wyss, C Liu, RM Walker, SLR Kardia, TS Wingo, J Ma, A Campbell, AP Wingo, T Huan, A Keshawarz, DA Bennett, JA Smith, KL Evans, D Levy and SJ London contributed to study design and/or supervised data analysis in each cohort. M Lee wrote the first draft of the manuscript. All authors contributed to interpretation of the results and/or critical revision of the manuscript. All authors approved the final version.

Ethical conduct of research

Each study was approved by its institutional review board. All study participants provided written informed consent.

Acknowledgments

We thank Frank Day of NIEHS and Jianping Jin of Westat, Inc. (Durham, NC, USA) for providing expert assistance.

Supplementary data

To view the supplementary data that accompany this paper please visit the journal website at: www.tandfonline.com/doi/suppl/10.2217/epi-2022-0353

Financial & competing interests disclosure

This work was supported in part by the Intramural Research Program of the National Institutes of Health, National Institute of Environmental Health Science (Z01-ES043012). Infrastructure for the CHARGE Consortium is supported in part by the National Heart, Lung, and Blood Institute grant HL105756. The Agricultural Lung Health Study is supported by the Intramural Research Program of the National Institutes of Health, National Institute of Environmental Health Sciences (Z01-ES102385, Z01-ES049030, Z01-ES043012 and for ABW, contract no. HHSN273201600003I) and the National Cancer Institute (Z01-CP010119). This work was supported in part by American Recovery and Reinvestment Act (ARRA) funds through NIEHS contract N01-ES-55546. The Framingham Heart Study is funded by National Institutes of Health contract N01-HC-25195 and HHSN268201500001I. The laboratory work for this investigation was funded by the Division of Intramural Research, National Heart, Lung, and Blood Institute, National Institutes of Health. The analytical component of this project was funded by the Division of Intramural Research, National Heart, Lung, and Blood Institute, and the Center for Information Technology, National Institutes of Health, Bethesda, MD. Disclaimer: The views expressed in this manuscript are those of the authors and do not necessarily represent the views of the National Heart, Lung, and Blood Institute; the National Institutes of Health; or the U.S. Department of Health and Human Services. Generation Scotland received core support from the Chief Scientist Office of the Scottish Government Health Directorates (CZD/16/6) and the Scottish Funding Council (HR03006) and is currently supported by the Wellcome Trust (216767/Z/19/Z). Genotyping of the GS:SFHS samples was carried out by the Genetics Core Laboratory at the Edinburgh Clinical Research Facility, University of Edinburgh, Scotland, and was funded by the Medical Research Council UK and the Wellcome Trust (Wellcome Trust Strategic Award “STratifying Resilience and Depression Longitudinally” [STRADL] reference 104036/Z/14/Z). The DNA methylation profiling and analysis were supported by Wellcome Investigator Award 220857/Z/20/Z and grant 104036/Z/14/Z (PI: AM McIntosh) and through funding from NARSAD (ref: 27404; awardee: DM Howard) and the Royal College of Physicians of Edinburgh (Sim Fellowship; awardee: HC Whalley). Support for the Genetic Epidemiology Network of Arteriopathy was provided by the National Heart, Lung and Blood Institute (U01 HL054457, RC1 HL100185, R01 HL087660, R01 HL119443 and R01 HL133221). ROS/MAP is supported by P30AG10161, P30AG72975, R01AG15819, R01AG17917, U01AG46152, U01AG61356. ROS/MAP resources can be requested at https://www.radc.rush.edu. A Hüls is supported by the HERCULES Center (NIEHS P30ES019776). The authors declare that they have no relevant conflicts of interest.

No funded writing assistance has been used in the creation of this manuscript.

Data sharing statement

Complete meta-analysis results can be found at Zenodo 10.5281/zenodo.7545108.

Additional information

Funding

This work was supported in part by the Intramural Research Program of the National Institutes of Health, National Institute of Environmental Health Science (Z01-ES043012). Infrastructure for the CHARGE Consortium is supported in part by the National Heart, Lung, and Blood Institute grant HL105756. The Agricultural Lung Health Study is supported by the Intramural Research Program of the National Institutes of Health, National Institute of Environmental Health Sciences (Z01-ES102385, Z01-ES049030, Z01-ES043012 and for ABW, contract no. HHSN273201600003I) and the National Cancer Institute (Z01-CP010119). This work was supported in part by American Recovery and Reinvestment Act (ARRA) funds through NIEHS contract N01-ES-55546. The Framingham Heart Study is funded by National Institutes of Health contract N01-HC-25195 and HHSN268201500001I. The laboratory work for this investigation was funded by the Division of Intramural Research, National Heart, Lung, and Blood Institute, National Institutes of Health. The analytical component of this project was funded by the Division of Intramural Research, National Heart, Lung, and Blood Institute, and the Center for Information Technology, National Institutes of Health, Bethesda, MD. Disclaimer: The views expressed in this manuscript are those of the authors and do not necessarily represent the views of the National Heart, Lung, and Blood Institute; the National Institutes of Health; or the U.S. Department of Health and Human Services. Generation Scotland received core support from the Chief Scientist Office of the Scottish Government Health Directorates (CZD/16/6) and the Scottish Funding Council (HR03006) and is currently supported by the Wellcome Trust (216767/Z/19/Z). Genotyping of the GS:SFHS samples was carried out by the Genetics Core Laboratory at the Edinburgh Clinical Research Facility, University of Edinburgh, Scotland, and was funded by the Medical Research Council UK and the Wellcome Trust (Wellcome Trust Strategic Award “STratifying Resilience and Depression Longitudinally” [STRADL] reference 104036/Z/14/Z). The DNA methylation profiling and analysis were supported by Wellcome Investigator Award 220857/Z/20/Z and grant 104036/Z/14/Z (PI: AM McIntosh) and through funding from NARSAD (ref: 27404; awardee: DM Howard) and the Royal College of Physicians of Edinburgh (Sim Fellowship; awardee: HC Whalley). Support for the Genetic Epidemiology Network of Arteriopathy was provided by the National Heart, Lung and Blood Institute (U01 HL054457, RC1 HL100185, R01 HL087660, R01 HL119443 and R01 HL133221). ROS/MAP is supported by P30AG10161, P30AG72975, R01AG15819, R01AG17917, U01AG46152, U01AG61356. ROS/MAP resources can be requested at http://www.radc.rush.edu. A Hüls is supported by the HERCULES Center (NIEHS P30ES019776). The authors declare that they have no relevant conflicts of interest. No funded writing assistance has been used in the creation of this manuscript.