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,4,5,6,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,4,5,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 (rg) with other substance use disorders (SUDs) (e.g. alcohol use disorder [AUD]; rg = 0.73) and psychiatric disorders (e.g. attention-deficit hyperactivity disorder; [rg = 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,11,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 rg between genetically-related traits (e.g., rg > 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

Data and participants

The meta-analysis includes summary statistics across seven cohorts examining OUD case vs. OUD control status in AFR and EUR ancestry individuals. We included both published and unpublished OUD GWAS. Previously published GWAS include data from Yale-Penn [3, 6, 7], PGC-SUD [6], and the Partners Biobank [15]. For MVP Releases 1 and 2 (the data releases used in the present analysis), a previous GWAS of OUD cases vs. opioid-exposed controls was reported [7]. MVP data included in the current meta-analysis use a different control definition (unscreened controls) to align better with the control definitions available in most other included samples. GWAS summary data for FinnGen [16] was accessed via a publicly available repository (https://r5.finngen.fi/). GWAS of OUD from iPSYCH [17], BioVU [18], and newly-available data from Yale-Penn subjects (Yale-Penn 3), previously unpublished, were performed by analysts at their respective study sites (Supplemental Materials). We had a total AFR N = 84,877 (Ncase = 5435 Neffective = 20,032), a total EUR N = 554,186 (Ncase = 15,251; Neffective = 56,994), and an overall N = 639,063 (Ncase =20,686; Neffective = 77,026). Other than Yale-Penn, this study involved de-identified data. The work was approved as appropriate by the Central Veterans Affairs (VA) institutional review board (IRB) and site-specific IRBs, including Yale University School of Medicine and VA Connecticut, and was conducted in accordance with all relevant ethical regulations. Cohort-specific summaries of AFR and EUR OUD subjects are presented in Table 1. Specific OUD diagnostic codes are provided in Supplementary Table 1. Additional phenotyping considerations are described in Supplemental Materials.

Table 1 Overview of samples included in GWAS meta-analysis of OUD cases vs. OUD controls.

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, I2) 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 Neffective). 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).

Fig. 1: Summary of OUD GWAS, meta-analysis, and MTAG study design.
figure 1

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

Data from the 1000 Genome Project (1000 G) phase 3 [20] was used for LD reference. Variants were mapped to the nearest gene based upon physical position (<10 kb 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 × 1008) OPRM1 variants in low LD (r2 < 0.1).

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

EUR OUD GWAS summary statistics were used to estimate SNP-heritability (h2SNP), and to characterize OUD genetic correlations (rg) using LD score regression (LDSC) [22]. LDSC analyses were restricted to HapMap3 variants [23]. Effective sample-size was used in all LDSC-based analytic steps. Genetic correlations were estimated with 54 other traits including SUDs, substance use, psychiatric traits, chronic pain, sociodemographic factors, and additional traits of interest (data sources are described in Supplemental Tables). Bonferroni-corrected significance was p = 9.26 × 10−04 (0.05/54). LDSC analyses were not performed in AFR and cross-ancestry meta-analyses because of the inability to use an LD reference panel for recently-admixed populations such as African-Americans or for analyses integrating datasets from diverse ancestry groups [22].

Multi-trait analysis of GWAS summary statistics (MTAG)

Based on LDSC estimates of genetic correlations with OUD, a joint-analysis that included the EUR OUD GWAS and GWAS summary statistics for AUD [11] and CanUD [12] was conducted using MTAG [13]. MTAG enhances statistical power by leveraging the genetic correlation between traits to generate trait-specific estimates for each SNP. The AUD GWAS summary statistics used in the present analysis are from a broader GWAS of problematic alcohol use [11]. MTAG used study-specific effective sample sizes for the respective GWAS. Study-specific effect sizes were transformed to Z-scores to be on a uniform scale across the three included GWAS. Included genetic variants were filtered using default MTAG parameters [13]. Briefly, variants were restricted to those common to all three of the GWAS, with a minor allele frequency (MAF) > 0.01, and present in at least two-thirds of the 90th percentile of the study-specific SNP sample sizes. These MTAG parameters guard against heterogeneity in the distribution of common vs. rare variant effects, ensuring that SNP effects generated from relatively small subsets of the contributing discovery GWAS do not bias the effect estimates across traits [13].

Phenome-wide Association Study (PheWAS)

To examine phenome-wide relationships for OUD and the OUD-MTAG analysis, and to compare their relationships with other clinically-relevant outcomes, we performed phenome-wide association studies (PheWAS) in BioVU [18], a cohort of >66,000 genotyped patients, with phenotypic data currently available for 1338 clinical outcomes from electronic health records [18]. Additional details on the BioVU cohort are provided in Supplemental Materials. Polygenic risk scores (PRS) for OUD and OUD-MTAG were computed using PRS-CS [24], excluding the subset of BioVU participants included in the meta-analysis. The respective PRS were then included in individual logistic regression models regressed on 1291 clinical outcomes with case counts ≥100, covarying for sex, age, and the first 10 genetic principal components. Statistical significance for the PheWAS was defined as p = 3.87 × 1005 (0.05/1291).

Polygenic risk score analysis

PRS are described in Supplemental Materials. Briefly, a leave-one-out PRS analysis was performed by excluding the EUR and AFR Yale-Penn 3 (YP3) cohorts from the respective OUD GWAS and OUD-MTAG analyses allowing for YP3 OUD cases and controls to be used as ancestry-specific PRS target samples.

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 × 1010) mapped to FURIN on chromosome 15, one of two GWS SNPs in the FURIN gene (rs17514846; r2 = 0.91). The second strongest association was with the OPRM1 functional variant (rs1799971; p = 4.92 × 1009). An additional OPRM1 variant was also identified (rs79704991; p = 1.11 × 1008; r2 = 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-pconditioned = 1.66 × 1006; rs79704991-pconditioned = 3.71 × 1006); 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).

Fig. 2: OUD and OUD-MTAG manhattan plots.
figure 2

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

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 × 1008), an intronic variant. Rs1799971 was also GWS in the cross-ancestry meta-analysis (p = 4.91 × 1008), and is not in strong LD with rs9478500 (EUR r2 = 0.03; AFR r2 = 0.002; ALL r2 = 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 Neffective) 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 × 1008).

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

Gene-based analysis

Gene-based analyses are described in Supplemental Materials. Both FURIN (p = 3.09 × 1007) and OPRM1 (p = 3.59 × 1007) were significant in EUR gene-based analysis (Supplementary Fig 4). Additional results are reported in Supplemental Table 35 and Supplementary Fig 5.

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

The liability scale SNP-heritability (h2SNP) 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.26 x 10−04) 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 (rg = 1.02; p = 2.38 x 10214) [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 (rg = 0.82; p = 1.14 x 1047) [12] and AUD (rg = 0.77; p = 6.36 x 1078) [11]. Modest genetic correlations were found for measures of substance use (e.g., the quantity/frequency alcohol use measure AUDIT-C) (rg = 0.14; p = 8.15 x 1003) [10].

Fig. 3: OUD genetic correlation results.
figure 3

EUR OUD GWAS genetic correlations (rg) 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 (rg = 0.52; p = 2.89 x 1018) and PTSD (rg = 0.52; p = 3.87 x 1019), lower back pain (rg = 0.61; p = 1.22 x 1009), inability to work due to being sick or disabled (rg = 0.57; p = 1.31 x 1020), and scores on the Townsend Deprivation Index (rg = 0.56; p = 1.13 x 1025). OUD was negatively genetically correlated with measures of sexual activity (age of first sexual intercourse [rg = −0.64; p = 4.43 x 1076]), indices of educational attainment (age of school completion [rg = −0.54; p = 9.41x1028]) and cognitive performance (rg = −0.38; p = 1.54 x 1020), and levels of past month “Heavy Do It Yourself” physical activity (rg = −0.38; p = 7.37 x 1013), 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 (rg=0.82; p = 1.14 x 1047) and AUD (rg = 0.77; p = 6.36 x 1078) in EUR. The OUD-MTAG analysis resulted in an increase in effective sample size from the original EUR OUD GWAS Neffective = 56,994 (GWAS mean x2 = 1.18) to an equivalent sample size of N = 128,748 (GWAS mean x2 = 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 1215 along with the OUD loci.

Table 3 Genome-wide significant (p ≤ 5.00 × 1008) Lead SNP associations in OUD-MTAG analysis (of 441 total GWS SNPs).

The top OUD-MTAG association was with rs11229119 (p = 7.03 x 1011) on chromosome 11 mapping to both TMX2 and CTNND1. The second strongest was with NICN1 (rs77648866; p = 1.82 x 1010) on chromosome 3. Additional GWS associations included FOXP2 (rs1989903; p = 2.47 x 1010), PDE4B (rs7519259; p = 2.68 x 1010), SLC39A8 (rs13135092; p = 3.60 x 1010), NCAM1 (rs1940701; p = 9.63 x 1010), RABEPK (rs864882; p = 1.24 x 1009), PLCL2 (rs55855024; p = 7.89 x 1009), and FTO (rs7188250; p = 3.63 x 1008). One of the FURIN variants identified in the EUR OUD GWAS was also GWS in the OUD-MTAG (rs17514846; p = 2.30 x 1008). The top OPRM1 association was with rs1799971 (p = 1.39 x 1006). 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.14 x 10−06) genes (Supplementary Table 7; Supplementary Fig 8).

The OUD-MTAG GWAS was significantly genetically correlated (Bonferroni p ≤ 9.26 x 10−04) with 46 traits including the largest previously published GWAS of OUD to date at rg = 0.98 (p = 1.22 x 1077) [7]. All estimates of genetic correlation for the OUD-MTAG analysis can be found in Supplementary Table 8.

Phenome-wide Association Study (PheWAS)

The top PheWAS association for OUD was with substance addiction and disorders (OR = 1.53; p = 2.12 x 1069). Additional top OUD associations included tobacco use disorder (OR = 1.26; p = 3.38 x 1056), chronic pain (OR = 1.25; p = 2.32 x 1028), alcohol-related disorders (OR = 1.35; p = 1.04 x 1023), mood (OR = 1.13; p = 1.27 x 1022) and anxiety (OR = 1.14; p = 1.00 × 1021) disorders, viral hepatitis C (OR = 1.33; p = 3.04 x 1020), and suicidal ideation or attempt (OR = 1.49; p = 2.17 x 1019), amongst others (Supplementary Fig 9; Supplementary Table 9).

Similar patterns of association were found for the OUD-MTAG PheWAS. The top associations were with tobacco use disorder (OR = 1.30; p = 1.37 x 1068), substance addiction and disorders (OR = 1.42; p = 1.15 x 1046), and alcohol-related disorders (OR = 1.42; p = 7.36 × 1031). OUD-MTAG also demonstrated significant associations with mood (OR = 1.12; p = 7.76 x 1021) and anxiety (OR = 1.13; p = 1.45 x 1020) disorders, chronic pain (OR = 1.20; p = 2.51 x 1020), viral hepatitis C (OR = 1.32; p = 1.73 x 1018), and suicidal ideation or attempt (OR = 1.47; p = 4.52 x 1018) (Supplementary Fig 9; Supplementary Table 10).

OUD polygenic risk score (PRS) analysis

In Yale-Penn 3 (YP3) EUR (N = 1959, 440 OUD cases), the EUR OUD PRS was a significant predictor of OUD (beta = 0.45; s.e. = 0.058; p = 2.9 x 1013) with the PRS specifically accounting for 2.41% of OUD variance (Supplementary Table 11). The OUD-MTAG PRS was a stronger predictor of OUD (beta = 0.61;s.e. = 0.066; p = 2.0 x 1016) explaining 3.81% of OUD variance.

In YP3 AFR (N = 1017, 171 OUD cases), both AFR-derived and EUR-derived OUD PRS were generated. The EUR-derived PRS (beta = 0.19;s.e. = 0.0945; p = 0.042; R2 = 0.29%) explained a small proportion of OUD variance in YP3 AFR but the AFR-derived PRS (beta = 0.073;s.e. = 0.049;p = 0.136; R2 = 0.11%) was not a significant predictor.

Discussion

We present a large genetic study of OUD, with an overall sample size of 639,063 (EUR = 554,186; AFR = 84,877) individuals (Ncases = 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 (r2 = 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 (r2 = 0.02). Two OPRM1 variants were also GWS in the cross-ancestry OUD GWAS (rs1799971 and rs9478500; r2 = 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; r2 = 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 (h2SNP) of 12.75% (Z = 11.28) compared to the previous largest OUD GWAS (h2SNP = 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 (rg = 0.82) and AUD (rg = 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 (rg = 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].