Abstract

Introduction

Opioid use disorder (OUD) has a serious negative impact on public health and is a leading cause of preventable death [1]. Although opioid misuse and progression to OUD [2] are influenced by heritable factors, discovery of OUD risk loci has been limited [3–7]. Difficulties in advancing OUD genetic discovery are largely due to lack of adequately powered cohorts of genetically informative samples [8, 9].

Genome-wide association studies (GWAS) examining single nucleotide polymorphism (SNP) effects on OUD have been underpowered [8, 9]. Nevertheless, recent progress in GWAS of OUD include the identification and confirmation of a genome-wide significant (GWS) functional variant (rs1799971) in OPRM1 [7]. Earlier OUD GWAS identified associations with variation in several genes including KCNG2, KCNC1, APBB2, CNIH3, RGMA, and OPRM1 [3–6], but the validity of those associations remains largely untested due to the lack of powerful independent OUD cohorts. OUD GWAS have also demonstrated genetic correlations (r_g) with other substance use disorders (SUDs) (e.g. alcohol use disorder [AUD]; r_g=0.73) and psychiatric disorders (e.g. attention-deficit hyperactivity disorder; [r_g=0.36]) [7].

Large-scale GWAS meta-analyses have advanced discovery of novel loci for SUDs (e.g., AUD, problematic alcohol use (PAU), cannabis use disorder (CanUD) [10–12]. This study applies similar meta-analytic methods for OUD, combining GWAS effects across multiple studies and two ancestral groups.

Multi-trait methods (e.g., MTAG; Multi-trait analysis of GWAS) [13] have the potential to increase power. MTAG capitalizes on the r_g between genetically-related traits (e.g., r_g>0.70) to increase the equivalent sample size. MTAG is an attractive option for boosting power for sets of similar traits like SUDs [11, 14], and holds particular promise for disorders such as OUD for which only limited cases are available for analysis. MTAG can generate estimates of trait-specific effects that leverage information from multiple GWAS summary statistics while accounting for both known and unknown sample overlap across the discovery samples [13]. MTAG can maximize the genetic information available for OUD by leveraging the statistical power of GWAS of non-opioid SUDs.

We conducted a large-scale GWAS meta-analysis of OUD in samples of African (AFR) and European (EUR) ancestry individuals. We maximized the informativeness of the available samples by performing a multi-trait analysis that incorporates SUDs that are highly genetically correlated with OUD.

Methods

Ancestry-specific and cross-ancestry GWAS meta-analysis

GWAS samples were combined using an effective sample-size weighted meta-analysis in METAL[19]. Ancestry-specific and cross-ancestry meta-analyses were performed. Measures of cross-sample heterogeneity (Cochran’s Q, I²) and genomic inflation (λ_GC) were used to evaluate potential bias influenced by heterogeneity between cohorts or by population stratification. Included GWAS summary statistics were limited to variants present in at least 80% of the analysis-specific effective sample size (e.g., 80% of EUR N_effective). The 80% effective sample size inclusion threshold ensured that variant effects present only in smaller cohorts did not disproportionately influence the overall results. This effectively meant that a variant needed to be present in MVP, PGC-SUD, and at least one additional cohort for it to be included (Fig. 1).

Open in a separate window

Fig. 1

Summary of OUD GWAS, meta-analysis, and MTAG study design.

Overview of European-ancestry opioid use disorder (OUD) genome-wide association study and OUD multi-trait analysis.

Data from the 1000 Genome Project (1000G) phase 3 [20] was used for LD reference. Variants were mapped to the nearest gene based upon physical position (<10kb from assigned gene) and further characterized by gene-mapping approaches using expression quantitative trait locus (eQTL) associations and 3D chromatin interactions (Hi-C)(Supplemental Materials). Conditional analyses were conducted using GCTA-COJO [21] to examine the conditional independence of genome-wide significant (GWS; p=5.00×10⁻⁰⁸) OPRM1 variants in low LD (r²<0.1).

Results

Ancestry-specific and cross-ancestry GWAS meta-analyses

In the ancestry-specific analyses, there were three GWS variants (Fig. 2) in EUR (Table 2). The top association (rs11372849; p=9.54×10⁻¹⁰) mapped to FURIN on chromosome 15, one of two GWS SNPs in the FURIN gene (rs17514846; r²=0.91). The second strongest association was with the OPRM1 functional variant (rs1799971; p=4.92×10⁻⁰⁹). An additional OPRM1 variant was also identified (rs79704991; p=1.11×10⁻⁰⁸; r²=0.02)(OPRM1 regional plots—Supplementary Fig. 1). GCTA-COJO [21] was used for conditional analysis of the two GWS OPRM1 variants demonstrating low LD (rs1799971 conditioned on rs79704991 and vice versa). In these analyses, each variant fell below GWS when conditioning on the effect of the other (conditioned rs1799971-p_conditioned=1.66×10⁻⁰⁶; rs79704991-p_conditioned=3.71×10⁻⁰⁶); although, there were no statistically significant differences in effect estimates for the respective OPRM1 variants conditioned vs. unconditioned effects. No GWS variants were identified in the AFR GWAS (Supplementary Fig. 2).

Open in a separate window

Fig. 2

OUD and OUD-MTAG manhattan plots.

Manhattan plots of (A) European-ancestry OUD GWAS results and (B) OUD-MTAG multi-trait GWAS results.

Table 2

Genome-wide significant (p≤5.00E-08) associations in (A) EUR OUD GWAS, and (B) cross-ancestry OUD GWAS.

(A) EUR OUD analysis
Chr	Position	Marker	A1	A2	Gene	EUR MAF	EUR Z	EUR P-value	Direction	AFR MAF	AFR Z	AFR P-value
15	91419432	rs11372849	T	TC	FURIN	0.46	6.12	9.54E-10	+++???+?	-	-	-
6	154360797	rs1799971	A	G	OPRM1	0.13	5.85	4.92E-09	++++++++	0.03	0.83	0.407
15	91416550	rs17514846	A	C	FURIN	0.46	−5.77	7.87E-09	----+--?	0.18	−0.94	0.347
6	154319449	rs79704991	T	G	OPRM1	0.13	5.71	1.11E-08	++++++++	0.08	0.59	0.552
6	154315310	rs12200046	T	C	OPRM1	0.13	5.64	1.70E-08	++++++++	0.08	0.67	0.505
6	154309808	rs10499276	T	C	OPRM1	0.13	5.52	3.38E-08	++++++++	0.08	0.49	0.621
6	154304242	rs34069531	T	C	OPRM1	0.13	5.51	3.67E-08	++++++++	0.08	0.61	0.542
6	154377925	rs3778146	T	C	OPRM1	0.17	−5.5	3.84E-08	---+----	0.08	−1.54	0.125
6	154378223	rs9478500	T	C	OPRM1	0.17	−5.49	4.12E-08	---+----	0.08	−1.76	0.078
6	154379152	rs3823010	A	G	OPRM1	0.17	5.48	4.25E-08	+++-++++	0.08	1.59	0.111
6	154378739	rs3778147	A	G	OPRM1	0.17	5.46	4.83E-08	+++-++++	0.08	1.47	0.141

Open in a separate window

(B) Cross-ancestry OUD analysis
Chr	Position	Marker	A1	A2	Gene	Z	P-value	Direction	EUR MAF	EUR Z	EUR P-value	AFR MAF	AFR Z	AFR P-value
6	154378223	rs9478500	T	C	OPRM1	−5.62	1.95E-08	---+-----+-+	0.17	−5.49	4.12E-08	0.18	−1.76	0.078
6	154379934	rs9285542	T	C	OPRM1	5.56	2.73E-08	+++-+++++-++	0.17	5.43	5.61E-08	0.16	1.74	0.082
6	154379152	rs3823010	A	G	OPRM1	5.53	3.26E-08	+++-++++++--	0.17	5.48	4.25E-08	0.11	1.59	0.111
6	154381012	rs3778148	T	G	OPRM1	5.52	3.34E-08	+++-++++++-+	0.17	5.43	5.73E-08	0.10	1.68	0.094
6	154355100	rs6936615	A	G	OPRM1	−5.52	3.38E-08	---+--?---+-	-	-	-	0.10	−1.72	0.085
6	154377925	rs3778146	T	C	OPRM1	−5.51	3.54E-08	---+------++	0.17	−5.5	3.84E-08	0.12	−1.54	0.125
6	154383658	rs3778150	T	C	OPRM1	−5.49	3.97E-08	---+-----+-+	0.17	-5.38	7.31E-08	0.18	−1.69	0.091
6	154382139	rs3778149	C	G	OPRM1	−5.48	4.20E-08	---+------++	0.17	−5.44	5.42E-08	0.12	−1.58	0.114
6	154382473	rs7772959	A	G	OPRM1	5.48	4.34E-08	+++-++++++--	0.17	5.43	5.66E-08	0.12	1.58	0.114
6	154362254	rs9322445	A	G	OPRM1	−5.47	4.54E-08	---+------+-	0.17	−5.44	5.49E-08	0.10	−1.56	0.120
6	154382367	rs7773995	T	C	OPRM1	5.46	4.78E-08	+++-++++++--	0.17	5.42	5.84E-08	0.12	1.56	0.120
6	154360797	rs1799971	A	G	OPRM1	5.46	4.91E-08	+++++++++-++	0.13	5.85	4.92E-09	0.03	0.83	0.407
Bold=Lead SNP

Open in a separate window

The cross-ancestry OUD GWAS identified two GWS risk variants mapping to OPRM1 (Supplementary Fig. 3) (Table 2). The top association was with rs9478500 (p=1.95×10⁻⁰⁸), an intronic variant. Rs1799971 was also GWS in the cross-ancestry meta-analysis (p=4.91×10⁻⁰⁸), and is not in strong LD with rs9478500 (EUR r²=0.03; AFR r²=0.002; ALL r²=0.04). The top FURIN association in EUR (rs11372849) was uninformative in three of four AFR ancestry cohorts and did not meet the threshold (80% of N_effective) we set to be included in the analysis. The second top FURIN association (rs17514846) fell below GWS in the cross-ancestry GWAS (p=6.00×10⁻⁰⁸).

SNP-heritability and Linkage-Disequilibrium (LD) Score Regression

The liability scale SNP-heritability (h²_SNP) estimate was 12.75% (s.e.=1.1%) in EUR using effective sample-size adjusted prevalence rates and a population prevalence of 0.021 [25]. Genome-wide inflation was mild with respect to sample size and favored OUD polygenicity as indicated by the LDSC inflation factor (λ_GC=1.18), intercept=1.01 (s.e.=0.011), and attenuation ratio=0.05 (s.e.=0.049).

OUD showed statistically significant (p≤9.26x10⁻⁰⁴) genetic correlations with 40 traits including substance use, SUDs, psychiatric traits, pain outcomes, physical health, and sociodemographic characteristics (Fig. 3; Supplementary Table 6). The OUD trait in the current study was strongly genetically correlated with the largest published GWAS of OUD to date (r_g=1.02; p=2.38x10⁻²¹⁴) [7], suggesting that OUD is being captured consistently across the studies, as might be expected given the substantial overlap in OUD cases between the two studies, although the control definitions differed between analyses. OUD was also strongly genetically correlated with other SUDs, including CanUD (r_g=0.82; p=1.14x10⁻⁴⁷) [12] and AUD (r_g=0.77; p=6.36x10⁻⁷⁸) [11]. Modest genetic correlations were found for measures of substance use (e.g., the quantity/frequency alcohol use measure AUDIT-C) (r_g=0.14; p=8.15x10⁻⁰³) [10].

Open in a separate window

Fig. 3

OUD genetic correlation results.

EUR OUD GWAS genetic correlations (r_g) with mental health, pain, physical health, sociodemographic, and substance use traits of interest.

OUD also demonstrated statistically significant genetic correlations with many mental health, pain, physical health, and sociodemographic traits. The strongest positive correlations across the respective domains were with Generalized Anxiety Disorder (r_g=0.52; p=2.89x10⁻¹⁸) and PTSD (r_g=0.52; p=3.87x10⁻¹⁹), lower back pain (r_g=0.61; p=1.22x10⁻⁰⁹), inability to work due to being sick or disabled (r_g=0.57; p=1.31x10⁻²⁰), and scores on the Townsend Deprivation Index (r_g=0.56; p=1.13x10⁻²⁵). OUD was negatively genetically correlated with measures of sexual activity (age of first sexual intercourse [r_g=−0.64; p=4.43x10⁻⁷⁶]), indices of educational attainment (age of school completion [r_g=−0.54; p=9.41x10⁻²⁸]) and cognitive performance (r_g=−0.38; p=1.54x10⁻²⁰), and levels of past month “Heavy Do It Yourself” physical activity (r_g=−0.38; p=7.37x10⁻¹³), among others (Supplementary Table 6).

Multi-trait analysis of European GWAS summary statistics (MTAG)

MTAG was supported by strong genetic correlation for OUD with CanUD (r_g=0.82; p=1.14x10⁻⁴⁷) and AUD (r_g=0.77; p=6.36x10⁻⁷⁸) in EUR. The OUD-MTAG analysis resulted in an increase in effective sample size from the original EUR OUD GWAS N_effective=56,994 (GWAS mean x²=1.18) to an equivalent sample size of N=128,748 (GWAS mean x²=1.40). The increase resulted in the identification of 18 independent GWS OUD-MTAG risk loci (Fig. 2; Table 3), some previously associated at either the variant level, or that reside in genes associated with, psychiatric and substance use outcomes in previous GWAS. Seven of the OUD-MTAG loci mapped to the nearest gene via brain eQTL data and Hi-C interactions; 8 loci were not mapped to the nearest gene via brain eQTL and Hi-C data but were implicated with additional genes in their respective genomic regions. Some loci fell in complex genomic regions with many mapped genes. These OUD-MTAG loci regions are summarized in Supplemental Fig. 6, 7 and Supplementary Tables 12–15 along with the OUD loci.

Table 3

Genome-wide significant (p≤5.00×10⁻⁰⁸) Lead SNP associations in OUD-MTAG analysis (of 441 total GWS SNPs).

Chr	Pos	MarkerName	Allele 1	Allele 2		Gene(s)	Z	P-value	p-value in AUD	p-value in CanUD
11	57535966	rs11229119	T		C	TMX2-CTNND1	−6.50	7.03E-11	2.60E-07	6.68E-05
3	49469449	rs77648866	A		G	NICN1	6.33	1.82E-10	5.90E-07	5.89E-03
7	114137940	rs1989903	A		G	FOXP2	−6.20	2.47E-10	1.23E-03	3.52E-09
1	66434743	rs7519259	A		G	PDE4B	5.73	2.68E-10	9.94E-08	1.36E-05
4	103198082	rs13135092	A		G	SLC39A8	6.51	3.60E-10	4.89E-18	0.88
11	112869404	rs1940701	T		C	NCAM1	-6.10	9.63E-10	4.91E-06	1.67E-03
9	127968109	rs864882	T		C	RABEPK	−5.60	1.24E-09	2.45E-05	0.06
6	19076417	rs9350100	T		C	RP11-254A17.1	5.62	1.76E-09	3.37E-05	2.89E-03
1	91208451	rs2166171	T		C	BARHL2	−5.88	4.80E-09	1.18E-06	2.76E-05
3	16850764	rs55855024	A		C	PLCL2	5.73	7.89E-09	1.50E-06	0.20
7	75622281	rs6467958	T		C	TMEM120A	−5.50	1.06E-08	2.24E-05	4.02E-03
11	113477081	rs11214677	T		C	TMPRSS5	−5.65	1.16E-08	1.23E-04	1.65E-05
1	28989020	rs6667501	A		G	GMEB1	5.55	1.74E-08	3.66E-04	2.21E-03
19	45453763	rs10422888	A		G	CTB-129P6.11	5.55	1.99E-08	5.83E-05	6.22E-03
15	91416550	rs17514846	A		C	FURIN	−5.59	2.30E-08	2.85E-03	0.09
15	47645174	rs73403005	A		G	SEMA6D	−5.69	2.46E-08	6.06E-09	2.38E-03
13	96932868	rs2389631	A		C	HS6ST3	−5.61	2.53E-08	8.72E-05	0.06
16	53834607	rs7188250	T		C	FTO	5.64	3.63E-08	4.41E-12	0.61
									Bold=GWS in non-OUD MTAG sumstats

Open in a separate window

The top OUD-MTAG association was with rs11229119 (p=7.03x10⁻¹¹) on chromosome 11 mapping to both TMX2 and CTNND1. The second strongest was with NICN1 (rs77648866; p=1.82x10⁻¹⁰) on chromosome 3. Additional GWS associations included FOXP2 (rs1989903; p=2.47x10⁻¹⁰), PDE4B (rs7519259; p=2.68x10⁻¹⁰), SLC39A8 (rs13135092; p=3.60x10⁻¹⁰), NCAM1 (rs1940701; p=9.63x10⁻¹⁰), RABEPK (rs864882; p=1.24x10⁻⁰⁹), PLCL2 (rs55855024; p=7.89x10⁻⁰⁹), and FTO (rs7188250; p=3.63x10⁻⁰⁸). One of the FURIN variants identified in the EUR OUD GWAS was also GWS in the OUD-MTAG (rs17514846; p=2.30x10⁻⁰⁸). The top OPRM1 association was with rs1799971 (p=1.39x10⁻⁰⁶). Of the 18 GWS loci, three were GWS in the AUD GWAS and one was GWS in the CanUD GWAS used for MTAG (Table 3).

The OUD-MTAG gene-based analysis resulted in 66 Bonferroni significant (p≤0.05/15,927=3.14x10⁻⁰⁶) genes (Supplementary Table 7; Supplementary Fig 8).

The OUD-MTAG GWAS was significantly genetically correlated (Bonferroni p≤9.26x10⁻⁰⁴) with 46 traits including the largest previously published GWAS of OUD to date at r_g=0.98 (p=1.22x10⁻⁷⁷) [7]. All estimates of genetic correlation for the OUD-MTAG analysis can be found in Supplementary Table 8.

Discussion

We present a large genetic study of OUD, with an overall sample size of 639,063 (EUR=554,186; AFR=84,877) individuals (N_cases=20,686 [EUR=15,251; AFR=5435]). This study is the first to provide evidence of a GWS single-variant GWAS association between FURIN and OUD. We support findings from previous OUD GWAS implicating OPRM1 as a risk locus for OUD [7], including the coding variant rs1799971 and additional OPRM1 associations that remained statistically significant in a cross-ancestry analysis of AFR and EUR populations. We add evidence of gene-based associations with OUD and provide estimates of OUD SNP-heritability and genetic correlations with numerous traits. Further, we apply a multi-trait approach for OUD genetic discovery utilizing the high degree of genetic correlation across SUDs (OUD, AUD, CanUD) to increase power, yielding an equivalent sample size of 128,748 and 18 GWS OUD-MTAG risk loci. PheWAS of OUD and OUD-MTAG demonstrated similar patterns of associations across the phenome, and the OUD-MTAG PRS explained a larger amount of variance in OUD case status (3.81%) compared to the OUD PRS (2.41%), suggesting that the OUD-MTAG and OUD analyses are capturing similar phenomenon.

Compared to other complex psychiatric traits, there are comparatively small samples available for genetic analyses of SUDs, particularly those involving illegal substances (heroin, cocaine) [8, 9]. A strategy that increases statistical power by incorporating other sets of samples—for example, from GWAS of closely-related but non-identical traits such as other SUDs—could help advance our understanding of the genetic architecture of OUD. This study brought much more information to bear on the analysis of OUD risk variation, resulting in the identification of many more loci. These associations included two loci from the EUR OUD GWAS (OPRM1 and FURIN), and 18 loci identified in the OUD-MTAG analysis (also including FURIN). The OUD-MTAG loci did not include OPRM1. The absence from the MTAG analysis of any statistically significant association mapped to OPRM1, a locus that should be highly-specific to OUD, was unexpected.

FURIN was associated with OUD risk in both SNP-based and gene-based analyses. FURIN (Furin, Paired Basic Amino Acid Cleaving Enzyme) encodes the endoprotease furin enzyme that serves a primary role in regulating synaptic neuronal activity, including the synthesis of brain-derived neurotropic factor and regulation of neurotrophin levels in the brain [26]. Genetic variation in FURIN has been associated with multiple psychiatric outcomes including schizophrenia [27, 28] and studies examining genetic and phenotypic features shared between schizophrenia and bipolar disorder [29, 30]. The two top FURIN SNPs associated with OUD are in strong LD (r²=0.91); the second strongest FURIN association, rs17514846, has been significantly associated with multiple cardiovascular and hypertension outcomes [31, 32], and was also GWS in a GWAS of parents’ attained age (current age of parents or parental age at death) [33]. A statistically-significant FURIN gene-level association being driven by rs17514846 was reported between FURIN and opioid addiction [34]. In a targeted follow-up in the FURIN gene region, there was significant association between rs11372849 (also lead SNP in the current study) and opioid addiction. Accumulating evidence linking FURIN and opioid outcomes, including the FURIN GWAS associations reported herein along with evidence of gene-based associations with opioid addiction [34], reflect the high degree of co-morbidity between OUD and psychiatric and physical health traits.

Our findings support previous OUD GWAS implicating OPRM1 genetic variation in opioid addiction and OUD [7, 34] and extend GWS findings for OPRM1 as a risk factor in a cross-ancestry analysis. The top association in the EUR OUD GWAS was with the OPRM1 coding variant (rs1799971) identified in an earlier OUD GWAS, all cases of which are also included in the current study, plus an additional OPRM1 variant (rs79704991) in low-LD with rs1799971 (r²=0.02). Two OPRM1 variants were also GWS in the cross-ancestry OUD GWAS (rs1799971 and rs9478500; r²=0.02). Rs9478500 was previously GWS for opioid addiction in EUR [34]. There is evidence of OPRM1’s complex haplotype structure and the potential for multiple independent OPRM1 risk loci influencing risk for OUD and SUDs dating to 2006 [35]. Conditional analysis of the top OPRM1 variants (rs1799971 and rs79704991; r²=0.02) demonstrated that these variants are not independent as indicated by each variant falling below GWS when conditioned on the other; however, the variant effects remained nominally significant and there were no significant differences in the conditioned vs. unconditioned effect sizes. Future studies of larger cohorts with diverse ancestral backgrounds will be needed to distinguish the effects of OUD risk alleles across the OPRM1 region, including the known-functional rs1799971 variant which may or may not be the variant motivating previous findings.

We found an estimated SNP-heritability (h²SNP) of 12.75% (Z=11.28) compared to the previous largest OUD GWAS (h²SNP=11.30%; Z=6.27) [7]. However, the comparison between these two studies is not direct: the largest previous GWAS [7] used opioid-exposed controls, while we used a broader control definition, including not only individuals who were opioid-exposed, but subjects with no OUD assessment. This was necessitated by the fact that many of the available datasets did not define exposed controls and would have been excluded had we used the exposed control definition. Findings from the current study do not establish whether the control definition impacted the detection of genetic loci or the genetic architecture of OUD.

OUD was positively genetically correlated with other SUDs (CanUD, AUD) and psychiatric conditions (PTSD, depression, schizophrenia), with lower correlations for measures of substance use (as opposed to dependence; e.g., AUDIT-C), suggesting that OUD is more akin to measures of substance dependence than use per se. OUD was genetically correlated with multiple forms of chronic pain (e.g., lower back pain) and indicators of impairment (inability to work, decreased physical activity) and significantly genetically correlated with socioeconomic hardship (Townsend Deprivation Index) and lower levels of educational attainment. These patterns of genetic correlation are consistent with and may reflect high rates of co-occurrence of OUD with SUDs and psychiatric disorders in epidemiologic studies [25, 36]. Beyond epidemiologic estimates, SUDs and psychiatric disorders have also been demonstrated to be risk markers for severe opioid-related outcomes, such as opioid overdose [37]. Lower educational attainment and greater economic hardship have also been associated with higher rates of opioid overdose and opioid overdose-related deaths [38]. These patterns of genetic correlation are consistent with the complex clinical presentation of OUD.

We examined the utility of using MTAG to increase the information available from the limited number of genotyped OUD subjects. The OUD multi-trait analysis was feasible given the high genetic correlations with CanUD (r_g=0.82) and AUD (r_g=0.77) and increased by an order-of-magnitude the number of GWS risk loci detected. While this provides proof-of-concept for this approach given that many of the loci identified via OUD-MTAG have previously been implicated with psychiatric and substance use outcomes, OPRM1 was not GWS in the OUD-MTAG analysis, so increased detection may have come at the cost of specificity for OUD. However, only 4 of the 18 OUD-MTAG GWS associations were GWS in the respective AUD and CanUD GWAS used as MTAG instruments, so the MTAG results did not simply reflect the findings from AUD and CanUD GWAS. The OUD-MTAG PRS also outperformed the OUD PRS in predicting OUD case status in the Yale-Penn holdout analysis (OUD-MTAG PRS: 3.81%, OUD PRS: 2.41%) indicating that the OUD-MTAG is capturing OUD risk, and that the OUD-MTAG PRS is more powerful than the OUD PRS alone as also reflected by the comparative number of GWS loci identified in the respective GWAS.

A PheWAS across 1291 clinical outcomes also demonstrated convergent patterns of association between OUD and OUD-MTAG with common comorbidities (including SUDs, psychiatric traits, chronic pain, viral hepatitis C), supporting that these two analyses capture genetic factors that underlie similar clinical presentations and related impairment. Additionally, summary data from the OUD MTAG analysis including multiple SUDs was highly genetically correlated (r_g=0.98) with OUD [7], so it appears that the OUD-MTAG did capture genetic information relevant to OUD risk, though measuring the risk for OUD through a genetic liability for SUDs more broadly. That is, genetic risk for OUD may be a combination of a broader addiction liability (OUD-MTAG loci) combined with the opioid-specific genetic effects (e.g., OPRM1) that were found in the OUD single-trait analysis that are also influencing risk.

The distinction between substance-specific genetic effects and general SUD liability is of interest. Quantitative genetic studies have demonstrated both substance-specific influences, as well as heritable factors that contribute to SUDs more broadly [39, 40]. Up to 38% of variation in opioid dependence was reportedly accounted for by opioid-specific factors that were not shared with other SUDs [41]. Molecular genetic studies have begun to disentangle common vs. substance-specific genetic influences, reporting evidence to suggest the presence of a common unitary addiction factor that can account for risk across SUDs, in addition to substance-specific influences [42, 43]. Larger-scale OUD studies will be needed to parse genomic influences specific to OUD from those underlying risk for SUDs more broadly, but this will require many more genotyped OUD cases, because it cannot be accomplished via statistical methods alone.

The present study has limitations. Despite including all genotyped OUD subjects available, the OUD-only component of the present study is smaller than GWAS for other substance use behaviors [14] because OUD cases are underrepresented in available datasets. MTAG yielded a much larger sample, but at the apparent cost of a reduction in specificity marked by the non-significance of OPRM1 in the OUD-MTAG. To maximize sample size while maintaining OUD diagnosis to define case status in extant datasets, we used an unscreened control group, which although not optimal, allowed for the inclusion of additional subject cohorts. Important consideration must be given to OUD control definitions [6, 34]. Additionally, inadequate subject numbers limited our ability to identify risk variants in non-EUR populations. This must be addressed by purposeful recruitment of AFR and other non-EUR OUD subjects.

We report novel findings from a large-scale GWAS meta-analysis of OUD and employed multi-trait approaches that advanced discovery. These identified genomic risk factors for the development of OUD and its underlying biology highlight the need to assemble large OUD datasets that include individuals from diverse ancestral backgrounds. To advance our scientific understanding of OUD risk will require study of a range of opioid-related traits (e.g., clinically diagnosed OUD, non-dependent opioid use, and prescription painkiller use) [44].

Supplementary information

Acknowledgements

Author contributions

JDD and JG conceived the initial study design. JDD, HZ, FRW, ASH, ECJ, RLK, HJE, AA, RP, HRJ, and JG provided ongoing feedback and refinement of the analytic plan. JDD, HZ, MG, DFL, FRW, SS, ECJ, YZN, DD, ADB, VMR, MVJ, RLK, ACJ, HJE, AA, RP, HRK, and JG provided contributions to the individual GWAS included in the meta-analysis. JDD performed the GWAS meta-analysis, genetic correlation analysis, multi-trait analysis of GWAS, gene-based analysis, functional characterization analysis, and PRS analysis. SS and MVJ performed the PheWAS analysis. DFL performed the conditional analysis. JDD, HZ, DFL, FRW, SS, ASH, ECJ, RLK, HJE, AA, RP, HRK, and JG contributed to data interpretation. JDD, SS, MVJ, and JG contributed to data visualization. JDD drafted, and JDD and JG refined the manuscript. All authors critically evaluated and revised the manuscript.

Competing interests

HRK is a member of scientific advisory boards for Dicerna Pharmaceuticals, Sophrosyne Pharmaceuticals, and Enthion Pharmaceuticals; a consultant to Sophrosyne Pharmaceuticals; and a member of the American Society of Clinical Psychopharmacology’s Alcohol Clinical Trials Initiative, which during the past three years was supported by Alkermes, Dicerna, Ethypharm, Lundbeck, Mitsubishi, and Otsuka. HRK and JG are named as inventors on PCT patent application #15/878,640 entitled: “Genotype-guided dosing of opioid agonists,” filed January 24, 2018 and issued on January 26, 2021 as U.S. Patent No. 10,900,082. JG and RP were paid for editorial work for the journal “Complex Psychiatry.” All other authors report no biomedical financial interests or potential conflicts of interest.

Footnotes

References